|
| 1 | +#' Calculate CWRES (Conditional Weighted Residuals) |
| 2 | +#' |
| 3 | +#' Computes CWRES following the FOCE approximation method described in |
| 4 | +#' Hooker et al. (2007). The Jacobian of model predictions with respect to |
| 5 | +#' random effects (etas) is computed via central finite differences. |
| 6 | +#' |
| 7 | +#' CWRES = L^{-1} * (y - f(eta_hat) + F * eta_hat) |
| 8 | +#' |
| 9 | +#' where: |
| 10 | +#' - f(eta_hat) = individual predictions (IPRED) in the transformed space |
| 11 | +#' - F = Jacobian df/deta evaluated at eta_hat |
| 12 | +#' - L * L^T = V = F * Omega * F^T + Sigma (FOCE variance) |
| 13 | +#' - Sigma = diagonal residual error variance matrix |
| 14 | +#' |
| 15 | +#' @references Hooker AC, Staatz CE, Karlsson MO. Conditional weighted |
| 16 | +#' residuals (CWRES): a model diagnostic for the FOCE method. |
| 17 | +#' Pharm Res. 2007;24(12):2187-2197. |
| 18 | +#' |
| 19 | +#' @param eta_hat numeric vector of estimated etas (MAP estimates) |
| 20 | +#' @param ipred_raw numeric vector of raw (untransformed) individual |
| 21 | +#' predictions at eta_hat |
| 22 | +#' @param y numeric vector of observed values |
| 23 | +#' @param omega_full full omega covariance matrix |
| 24 | +#' @param error list with `prop` and `add` residual error components |
| 25 | +#' @param obs_type integer vector of observation types |
| 26 | +#' @param transf transformation function (identity or log for LTBS) |
| 27 | +#' @param model PKPDsim model object |
| 28 | +#' @param parameters_population list of population parameters |
| 29 | +#' @param nonfixed character vector of non-fixed parameter names |
| 30 | +#' @param as_eta character vector of parameters estimated directly as eta |
| 31 | +#' @param covariates list of PKPDsim covariates |
| 32 | +#' @param regimen PKPDsim regimen object |
| 33 | +#' @param lagtime lagtime specification |
| 34 | +#' @param t_obs numeric vector of observation times |
| 35 | +#' @param obs_type_sim observation type vector for simulation |
| 36 | +#' @param int_step_size integrator step size |
| 37 | +#' @param iov_bins IOV bin specification |
| 38 | +#' @param A_init initial state vector |
| 39 | +#' @param t_init initialization time |
| 40 | +#' @param steady_state_analytic steady state settings (or NULL) |
| 41 | +#' @param weight_prior_var prior weight variance scaling factor (default 1). |
| 42 | +#' Used to scale omega in the FOCE Hessian computation to match the |
| 43 | +#' estimation objective. CWRES uses the unscaled omega (model diagnostic). |
| 44 | +#' @param delta perturbation size for finite differences |
| 45 | +#' |
| 46 | +#' @return list with components: |
| 47 | +#' \item{cwres}{numeric vector of CWRES values} |
| 48 | +#' \item{vcov}{FOCE variance-covariance matrix of eta estimates |
| 49 | +#' (n_eta x n_eta), or NULL if computation failed. Derived from the |
| 50 | +#' same Jacobian F used for CWRES: vcov = (Omega^-1 + F' Sigma^-1 F)^-1} |
| 51 | +#' |
| 52 | +calc_cwres <- function( |
| 53 | + eta_hat, |
| 54 | + ipred_raw, |
| 55 | + y, |
| 56 | + omega_full, |
| 57 | + error, |
| 58 | + obs_type, |
| 59 | + transf, |
| 60 | + model, |
| 61 | + parameters_population, |
| 62 | + nonfixed, |
| 63 | + as_eta, |
| 64 | + covariates, |
| 65 | + regimen, |
| 66 | + lagtime, |
| 67 | + t_obs, |
| 68 | + obs_type_sim, |
| 69 | + int_step_size = 0.1, |
| 70 | + iov_bins = NULL, |
| 71 | + A_init = NULL, |
| 72 | + t_init = 0, |
| 73 | + steady_state_analytic = NULL, |
| 74 | + weight_prior_var = 1, |
| 75 | + delta = 1e-4 |
| 76 | +) { |
| 77 | + n_obs <- length(y) |
| 78 | + n_eta <- length(eta_hat) |
| 79 | + |
| 80 | + if (n_obs == 0 || n_eta == 0) { |
| 81 | + return(list(cwres = numeric(0), vcov = NULL)) |
| 82 | + } |
| 83 | + |
| 84 | + # Transformed individual predictions at eta_hat |
| 85 | + ipred_transf <- transf(ipred_raw) |
| 86 | + |
| 87 | + # Compute Jacobian F = d(transf(f))/d(eta) via central finite differences |
| 88 | + F_matrix <- matrix(0, nrow = n_obs, ncol = n_eta) |
| 89 | + for (j in seq_len(n_eta)) { |
| 90 | + eta_plus <- eta_hat |
| 91 | + eta_minus <- eta_hat |
| 92 | + eta_plus[j] <- eta_hat[j] + delta |
| 93 | + eta_minus[j] <- eta_hat[j] - delta |
| 94 | + |
| 95 | + pred_plus <- simulate_with_etas( |
| 96 | + eta_plus, parameters_population, nonfixed, as_eta, |
| 97 | + model, covariates, regimen, lagtime, t_obs, obs_type_sim, |
| 98 | + int_step_size, iov_bins, A_init, t_init, |
| 99 | + steady_state_analytic |
| 100 | + ) |
| 101 | + pred_minus <- simulate_with_etas( |
| 102 | + eta_minus, parameters_population, nonfixed, as_eta, |
| 103 | + model, covariates, regimen, lagtime, t_obs, obs_type_sim, |
| 104 | + int_step_size, iov_bins, A_init, t_init, |
| 105 | + steady_state_analytic |
| 106 | + ) |
| 107 | + |
| 108 | + F_matrix[, j] <- (transf(pred_plus) - transf(pred_minus)) / (2 * delta) |
| 109 | + } |
| 110 | + |
| 111 | + # Residual error variance diagonal (in transformed space) |
| 112 | + sigma_diag <- error$prop[obs_type]^2 * ipred_transf^2 + |
| 113 | + error$add[obs_type]^2 |
| 114 | + |
| 115 | + if (any(!is.finite(sigma_diag) | sigma_diag <= 0)) { |
| 116 | + bad <- which(!is.finite(sigma_diag) | sigma_diag <= 0) |
| 117 | + warning( |
| 118 | + "CWRES/vcov computation failed: residual error variance is zero, ", |
| 119 | + "negative, or non-finite for observation(s) ", |
| 120 | + paste(bad, collapse = ", "), |
| 121 | + " (obs_type: ", paste(unique(obs_type[bad]), collapse = ", "), ")." |
| 122 | + ) |
| 123 | + return(list(cwres = rep(NA_real_, n_obs), vcov = NULL)) |
| 124 | + } |
| 125 | + |
| 126 | + # Omega submatrix for estimated parameters |
| 127 | + omega_est <- omega_full[seq_len(n_eta), seq_len(n_eta), drop = FALSE] |
| 128 | + |
| 129 | + # FOCE approximate marginal variance: V = F * Omega * F^T + Sigma |
| 130 | + V <- F_matrix %*% omega_est %*% t(F_matrix) + |
| 131 | + diag(sigma_diag, nrow = n_obs) |
| 132 | + |
| 133 | + # Cholesky decomposition (lower triangular) |
| 134 | + L <- tryCatch( |
| 135 | + t(chol(V)), |
| 136 | + error = function(e) NULL |
| 137 | + ) |
| 138 | + cwres <- if (is.null(L)) { |
| 139 | + warning("CWRES computation failed: variance matrix is not positive definite.") |
| 140 | + rep(NA_real_, n_obs) |
| 141 | + } else { |
| 142 | + # CWRES = L^{-1} * (y_transf - ipred_transf + F * eta_hat) |
| 143 | + # Derivation: |
| 144 | + # E_FOCE(y) = f(eta_hat) - F * eta_hat (since E[eta] = 0) |
| 145 | + # y - E_FOCE(y) = y - f(eta_hat) + F * eta_hat |
| 146 | + y_transf <- transf(y) |
| 147 | + as.numeric( |
| 148 | + solve(L, y_transf - ipred_transf + F_matrix %*% eta_hat) |
| 149 | + ) |
| 150 | + } |
| 151 | + |
| 152 | + # FOCE Hessian: H = (Omega/w)^{-1} + F' * Sigma^{-1} * F |
| 153 | + # where w = weight_prior_var, matching the scaled omega used during MAP |
| 154 | + # estimation. vcov of etas = H^{-1}. |
| 155 | + # This reuses the Jacobian F already computed for CWRES, so no extra |
| 156 | + # simulations are needed (replaces the numDeriv::hessian computation). |
| 157 | + # With weight_prior_var <= 0 the prior is effectively flat (LS fit), so |
| 158 | + # the FOCE vcov is not meaningful; skip and let the caller fall back. |
| 159 | + vcov <- NULL |
| 160 | + if (weight_prior_var > 0) { |
| 161 | + omega_scaled <- omega_est / weight_prior_var |
| 162 | + omega_scaled_inv <- tryCatch(solve(omega_scaled), error = function(e) NULL) |
| 163 | + if (!is.null(omega_scaled_inv)) { |
| 164 | + H_foce <- omega_scaled_inv + |
| 165 | + t(F_matrix) %*% diag(1 / sigma_diag, nrow = n_obs) %*% F_matrix |
| 166 | + vcov <- tryCatch(solve(H_foce), error = function(e) { |
| 167 | + warning("FOCE variance-covariance computation failed.") |
| 168 | + NULL |
| 169 | + }) |
| 170 | + } |
| 171 | + } |
| 172 | + |
| 173 | + list(cwres = cwres, vcov = vcov) |
| 174 | +} |
| 175 | + |
| 176 | + |
| 177 | +#' Simulate predictions with a given eta vector |
| 178 | +#' |
| 179 | +#' Helper for computing the Jacobian in CWRES calculation via finite |
| 180 | +#' differences. Computes individual parameters from etas and runs a |
| 181 | +#' PKPDsim simulation. |
| 182 | +#' |
| 183 | +#' @keywords internal |
| 184 | +simulate_with_etas <- function( |
| 185 | + eta, |
| 186 | + parameters_population, |
| 187 | + nonfixed, |
| 188 | + as_eta, |
| 189 | + model, |
| 190 | + covariates, |
| 191 | + regimen, |
| 192 | + lagtime, |
| 193 | + t_obs, |
| 194 | + obs_type, |
| 195 | + int_step_size, |
| 196 | + iov_bins, |
| 197 | + A_init, |
| 198 | + t_init, |
| 199 | + steady_state_analytic |
| 200 | +) { |
| 201 | + # Compute individual parameters from etas |
| 202 | + par <- parameters_population |
| 203 | + for (i in seq_along(nonfixed)) { |
| 204 | + key <- nonfixed[i] |
| 205 | + if (key %in% as_eta) { |
| 206 | + par[[key]] <- eta[i] |
| 207 | + } else { |
| 208 | + par[[key]] <- par[[key]] * exp(eta[i]) |
| 209 | + } |
| 210 | + } |
| 211 | + |
| 212 | + # Recompute steady-state A_init if needed |
| 213 | + a_init <- A_init |
| 214 | + if (!is.null(steady_state_analytic)) { |
| 215 | + a_init <- PKPDsim::calc_ss_analytic( |
| 216 | + f = steady_state_analytic$f, |
| 217 | + dose = regimen$dose_amts[1], |
| 218 | + interval = regimen$interval[1], |
| 219 | + model = model, |
| 220 | + parameters = par, |
| 221 | + covariates = covariates, |
| 222 | + map = steady_state_analytic$map, |
| 223 | + n_transit_compartments = PKPDsim::ifelse0( |
| 224 | + steady_state_analytic$n_transit_compartments, FALSE |
| 225 | + ), |
| 226 | + auc = PKPDsim::ifelse0(steady_state_analytic$auc, FALSE) |
| 227 | + ) |
| 228 | + } |
| 229 | + |
| 230 | + suppressMessages({ |
| 231 | + sim <- PKPDsim::sim_ode( |
| 232 | + ode = model, |
| 233 | + parameters = par, |
| 234 | + mixture_group = NULL, |
| 235 | + covariates = covariates, |
| 236 | + n_ind = 1, |
| 237 | + int_step_size = int_step_size, |
| 238 | + regimen = regimen, |
| 239 | + t_obs = t_obs, |
| 240 | + obs_type = obs_type, |
| 241 | + only_obs = TRUE, |
| 242 | + checks = FALSE, |
| 243 | + A_init = a_init, |
| 244 | + iov_bins = iov_bins, |
| 245 | + t_init = t_init, |
| 246 | + lagtime = lagtime |
| 247 | + ) |
| 248 | + }) |
| 249 | + |
| 250 | + sim$y |
| 251 | +} |
0 commit comments