Add proper CWRES computation using FOCE linearization (Hooker et al. 2007)

Replace the incorrect ad-hoc CWRES normalization with proper FOCE-based
conditional weighted residuals. The Jacobian F = df/deta is computed via
central finite differences and used to derive both CWRES and the FOCE
variance-covariance matrix, replacing the more expensive numDeriv::hessian
computation when residuals=TRUE.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: roninsightrx

DESCRIPTION

@@ -12,6 +12,6 @@ Depends: utils, stats
 Imports: PKPDsim, numDeriv, mvtnorm
 Suggests: testthat (>= 3.0.0)
 Remotes: InsightRX/PKPDsim
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3
 Encoding: UTF-8
 Config/testthat/edition: 3

R/calc_cwres.R

@@ -0,0 +1,234 @@
+#' Calculate CWRES (Conditional Weighted Residuals)
+#'
+#' Computes CWRES following the FOCE approximation method described in
+#' Hooker et al. (2007). The Jacobian of model predictions with respect to
+#' random effects (etas) is computed via central finite differences.
+#'
+#' CWRES = L^{-1} * (y - f(eta_hat) + F * eta_hat)
+#'
+#' where:
+#' - f(eta_hat) = individual predictions (IPRED) in the transformed space
+#' - F = Jacobian df/deta evaluated at eta_hat
+#' - L * L^T = V = F * Omega * F^T + Sigma (FOCE variance)
+#' - Sigma = diagonal residual error variance matrix
+#'
+#' @references Hooker AC, Staatz CE, Karlsson MO. Conditional weighted
+#'   residuals (CWRES): a model diagnostic for the FOCE method.
+#'   Pharm Res. 2007;24(12):2187-2197.
+#'
+#' @param eta_hat numeric vector of estimated etas (MAP estimates)
+#' @param ipred_raw numeric vector of raw (untransformed) individual
+#'   predictions at eta_hat
+#' @param y numeric vector of observed values
+#' @param omega_full full omega covariance matrix
+#' @param error list with `prop` and `add` residual error components
+#' @param obs_type integer vector of observation types
+#' @param transf transformation function (identity or log for LTBS)
+#' @param model PKPDsim model object
+#' @param parameters_population list of population parameters
+#' @param nonfixed character vector of non-fixed parameter names
+#' @param as_eta character vector of parameters estimated directly as eta
+#' @param covariates list of PKPDsim covariates
+#' @param regimen PKPDsim regimen object
+#' @param lagtime lagtime specification
+#' @param t_obs numeric vector of observation times
+#' @param obs_type_sim observation type vector for simulation
+#' @param int_step_size integrator step size
+#' @param iov_bins IOV bin specification
+#' @param A_init initial state vector
+#' @param t_init initialization time
+#' @param steady_state_analytic steady state settings (or NULL)
+#' @param delta perturbation size for finite differences
+#' @param ... additional arguments passed to PKPDsim::sim_ode
+#'
+#' @return list with components:
+#'   \item{cwres}{numeric vector of CWRES values}
+#'   \item{vcov}{FOCE variance-covariance matrix of eta estimates
+#'     (n_eta x n_eta), or NULL if computation failed. Derived from the
+#'     same Jacobian F used for CWRES: vcov = (Omega^-1 + F' Sigma^-1 F)^-1}
+#'
+calc_cwres <- function(
+    eta_hat,
+    ipred_raw,
+    y,
+    omega_full,
+    error,
+    obs_type,
+    transf,
+    model,
+    parameters_population,
+    nonfixed,
+    as_eta,
+    covariates,
+    regimen,
+    lagtime,
+    t_obs,
+    obs_type_sim,
+    int_step_size = 0.1,
+    iov_bins = NULL,
+    A_init = NULL,
+    t_init = 0,
+    steady_state_analytic = NULL,
+    delta = 1e-4,
+    ...
+) {
+  n_obs <- length(y)
+  n_eta <- length(eta_hat)
+
+  if (n_obs == 0 || n_eta == 0) {
+    return(list(cwres = numeric(0), vcov = NULL))
+  }
+
+  # Transformed individual predictions at eta_hat
+  ipred_transf <- transf(ipred_raw)
+
+  # Compute Jacobian F = d(transf(f))/d(eta) via central finite differences
+  F_matrix <- matrix(0, nrow = n_obs, ncol = n_eta)
+  for (j in seq_len(n_eta)) {
+    eta_plus <- eta_hat
+    eta_minus <- eta_hat
+    eta_plus[j] <- eta_hat[j] + delta
+    eta_minus[j] <- eta_hat[j] - delta
+
+    pred_plus <- simulate_with_etas(
+      eta_plus, parameters_population, nonfixed, as_eta,
+      model, covariates, regimen, lagtime, t_obs, obs_type_sim,
+      int_step_size, iov_bins, A_init, t_init,
+      steady_state_analytic, ...
+    )
+    pred_minus <- simulate_with_etas(
+      eta_minus, parameters_population, nonfixed, as_eta,
+      model, covariates, regimen, lagtime, t_obs, obs_type_sim,
+      int_step_size, iov_bins, A_init, t_init,
+      steady_state_analytic, ...
+    )
+
+    F_matrix[, j] <- (transf(pred_plus) - transf(pred_minus)) / (2 * delta)
+  }
+
+  # Residual error variance diagonal (in transformed space)
+  sigma_diag <- error$prop[obs_type]^2 * ipred_transf^2 +
+    error$add[obs_type]^2
+
+  # Omega submatrix for estimated parameters
+  omega_est <- omega_full[seq_len(n_eta), seq_len(n_eta), drop = FALSE]
+
+  # FOCE approximate marginal variance: V = F * Omega * F^T + Sigma
+  V <- F_matrix %*% omega_est %*% t(F_matrix) +
+    diag(sigma_diag, nrow = n_obs)
+
+  # Cholesky decomposition (lower triangular)
+  L <- tryCatch(
+    t(chol(V)),
+    error = function(e) NULL
+  )
+  cwres <- if (is.null(L)) {
+    warning("CWRES computation failed: variance matrix is not positive definite.")
+    rep(NA_real_, n_obs)
+  } else {
+    # CWRES = L^{-1} * (y_transf - ipred_transf + F * eta_hat)
+    # Derivation:
+    #   E_FOCE(y) = f(eta_hat) - F * eta_hat  (since E[eta] = 0)
+    #   y - E_FOCE(y) = y - f(eta_hat) + F * eta_hat
+    y_transf <- transf(y)
+    as.numeric(
+      solve(L, y_transf - ipred_transf + F_matrix %*% eta_hat)
+    )
+  }
+
+  # FOCE Hessian: H = Omega^{-1} + F' * Sigma^{-1} * F
+  # vcov of etas = H^{-1}
+  # This reuses the Jacobian F already computed for CWRES, so no extra
+  # simulations are needed (replaces the numDeriv::hessian computation).
+  omega_est_inv <- tryCatch(solve(omega_est), error = function(e) NULL)
+  vcov <- NULL
+  if (!is.null(omega_est_inv)) {
+    H_foce <- omega_est_inv +
+      t(F_matrix) %*% diag(1 / sigma_diag, nrow = n_obs) %*% F_matrix
+    vcov <- tryCatch(solve(H_foce), error = function(e) {
+      warning("FOCE variance-covariance computation failed.")
+      NULL
+    })
+  }
+
+  list(cwres = cwres, vcov = vcov)
+}
+
+
+#' Simulate predictions with a given eta vector
+#'
+#' Helper for computing the Jacobian in CWRES calculation via finite
+#' differences. Computes individual parameters from etas and runs a
+#' PKPDsim simulation.
+#'
+#' @keywords internal
+simulate_with_etas <- function(
+    eta,
+    parameters_population,
+    nonfixed,
+    as_eta,
+    model,
+    covariates,
+    regimen,
+    lagtime,
+    t_obs,
+    obs_type,
+    int_step_size,
+    iov_bins,
+    A_init,
+    t_init,
+    steady_state_analytic,
+    ...
+) {
+  # Compute individual parameters from etas
+  par <- parameters_population
+  for (i in seq_along(nonfixed)) {
+    key <- nonfixed[i]
+    if (key %in% as_eta) {
+      par[[key]] <- eta[i]
+    } else {
+      par[[key]] <- par[[key]] * exp(eta[i])
+    }
+  }
+
+  # Recompute steady-state A_init if needed
+  a_init <- A_init
+  if (!is.null(steady_state_analytic)) {
+    a_init <- PKPDsim::calc_ss_analytic(
+      f = steady_state_analytic$f,
+      dose = regimen$dose_amts[1],
+      interval = regimen$interval[1],
+      model = model,
+      parameters = par,
+      covariates = covariates,
+      map = steady_state_analytic$map,
+      n_transit_compartments = PKPDsim::ifelse0(
+        steady_state_analytic$n_transit_compartments, FALSE
+      ),
+      auc = PKPDsim::ifelse0(steady_state_analytic$auc, FALSE)
+    )
+  }
+
+  suppressMessages({
+    sim <- PKPDsim::sim_ode(
+      ode = model,
+      parameters = par,
+      mixture_group = NULL,
+      covariates = covariates,
+      n_ind = 1,
+      int_step_size = int_step_size,
+      regimen = regimen,
+      t_obs = t_obs,
+      obs_type = obs_type,
+      only_obs = TRUE,
+      checks = FALSE,
+      A_init = a_init,
+      iov_bins = iov_bins,
+      t_init = t_init,
+      lagtime = lagtime,
+      ...
+    )
+  })
+
+  sim$y
+}

R/calc_residuals.R

@@ -10,7 +10,10 @@
 #' @param A_init_individual init vector for model state (individual)
 #' @param censoring_idx vector with indices for censoring
 #' @param data_before_init data.frame with data before initial dose
-#' 
+#' @param nonfixed character vector of non-fixed parameter names
+#' @param as_eta character vector of parameters estimated directly as eta
+#' @param steady_state_analytic steady state settings (or NULL)
+#'
 calc_residuals <- function(
   obj,
   data,
@@ -33,6 +36,9 @@ calc_residuals <- function(
   censoring_idx = NULL,
   data_before_init = NULL,
   ltbs = FALSE,
+  nonfixed = NULL,
+  as_eta = c(),
+  steady_state_analytic = NULL,
   ...
 ) {
 
@@ -100,7 +106,59 @@ calc_residuals <- function(
     data_before_init,
     weights
   )
-  
+
+  ## Compute proper CWRES and FOCE vcov using the Jacobian (Hooker et al. 2007).
+  ## The Jacobian F = df/deta is computed once via finite differences (2*n_eta
+  ## simulations). Both CWRES and the FOCE Hessian/vcov are derived from F,
+  ## instead of the more expensive numDeriv::hessian computation.
+  if (!is.null(nonfixed) && length(nonfixed) > 0) {
+    n_before <- nrow(data_before_init)
+    # Extract ipred for real observations only (excluding before-init)
+    ipred_real <- sim_ipred$y[(n_before + 1):length(sim_ipred$y)]
+
+    foce_result <- tryCatch(
+      calc_cwres(
+        eta_hat = obj$fit$coef,
+        ipred_raw = ipred_real,
+        y = data$y,
+        omega_full = omega_full,
+        error = error,
+        obs_type = data$obs_type,
+        transf = transf,
+        model = model,
+        parameters_population = parameters_population,
+        nonfixed = nonfixed,
+        as_eta = as_eta,
+        covariates = covariates,
+        regimen = regimen,
+        lagtime = lagtime,
+        t_obs = data$t,
+        obs_type_sim = data$obs_type,
+        int_step_size = int_step_size,
+        iov_bins = iov_bins,
+        A_init = A_init_individual,
+        t_init = t_init,
+        steady_state_analytic = steady_state_analytic,
+        ...
+      ),
+      error = function(e) {
+        warning("CWRES computation failed: ", e$message)
+        list(cwres = rep(NA_real_, length(data$y)), vcov = NULL)
+      }
+    )
+
+    # Pad with zeros for before-init observations and apply weights
+    obj$cwres <- c(rep(0, n_before), foce_result$cwres * weights)
+
+    # Set censored observations to NA
+    if (!is.null(censoring) && any(censoring_idx)) {
+      obj$cwres[censoring_idx] <- NA_real_
+    }
+
+    # Store FOCE vcov for use in get_map_estimates
+    obj$foce_vcov <- foce_result$vcov
+  }
+
   ## Add covariates and parameters to obj
   if(output_include$covariates && !is.null(covariates)) {
     obj$covariates_time <- sim_ipred[!duplicated(sim_ipred$t), names(covariates)]

R/get_map_estimates.R

@@ -302,14 +302,18 @@ get_map_estimates <- function(
       likelihood = like
     )
   } else {
+    # When residuals=TRUE, the FOCE Jacobian computed in calc_residuals
+    # provides both CWRES and the Hessian/vcov, so we can skip the more
+    # expensive numDeriv::hessian here.
+    skip_hessian_mle <- skip_hessian || residuals
     output <- tryCatch({
       fit <- mle_wrapper(
         ll_func,
         start = omega$eta,
         method = method,
         optimizer = optimizer,
         control = control,
-        skip_hessian = skip_hessian,
+        skip_hessian = skip_hessian_mle,
         data = list(
           data = data,
           sim_object = sim_object,
@@ -427,13 +431,22 @@ get_map_estimates <- function(
       censoring_idx = censoring_idx,
       data_before_init = data_before_init,
       ltbs = ltbs,
+      nonfixed = omega$nonfixed,
+      as_eta = as_eta,
+      steady_state_analytic = steady_state_analytic,
       ...
     )
   }
 
-  ## Add variance-covariance matrix to output object
+  ## Add variance-covariance matrix to output object.
+  ## When residuals were computed, the FOCE vcov from the Jacobian is
+  ## preferred over numDeriv::hessian (faster, same FOCE approximation).
+  vcov_source <- obj$fit$vcov
+  if (!is.null(obj$foce_vcov)) {
+    vcov_source <- obj$foce_vcov
+  }
   obj$vcov_full <- get_varcov_matrix(
-    obj$fit$vcov, 
+    vcov_source,
     fallback = omega$full
   )
   if(inherits(obj$vcov_full, "matrix")) {

R/parse_residuals_from_predictions.R

@@ -39,8 +39,9 @@ parse_residuals_from_predictions <- function(
   obj$res <- (transf(y) - transf(pred))
   obj$weights <- c(rep(0, length(data_before_init$t)), weights)
   obj$wres <- (obj$res / w_pred) * obj$weights
-  obj$cwres <- obj$res / sqrt(abs(cov(transf(pred), transf(y)))) * c(rep(0, nrow(data_before_init)), obj$weights)
-  # Note: in NONMEM CWRES is on the population level, so can't really compare. NONMEM calls this CIWRES, it seems.
+  # CWRES is computed separately in calc_residuals() using the FOCE
+  # linearization method (Hooker et al. 2007). Initialize to NA here.
+  obj$cwres <- rep(NA_real_, length(y))
   obj$ires <- (transf(y) - transf(ipred))
   obj$iwres <- (obj$ires / w_ipred)
   obj$w_ipred  <- w_ipred