RXR-3164: Pass observation-level weights to covariance matrix calcs, and re-activate FOCE-type vcov

R/calc_cwres.R

@@ -41,13 +41,18 @@
 #' @param weight_prior_var prior weight variance scaling factor (default 1).
 #'   Used to scale omega in the FOCE Hessian computation to match the
 #'   estimation objective. CWRES uses the unscaled omega (model diagnostic).
+#' @param weights numeric vector of per-observation weights (length n_obs,
+#'   default all 1). Must match the weights used in MAP estimation so that
+#'   the FOCE Hessian matches the curvature of the actual OFV. CWRES (a model
+#'   diagnostic) always uses unit weights per Hooker et al. 2007.
 #' @param delta perturbation size for finite differences
 #'
 #' @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}
+#'     same Jacobian F used for CWRES: vcov = (Omega^-1 + F' Sigma^-1 F)^-1,
+#'     with per-observation weights applied to the likelihood term.}
 #'
 calc_cwres <- function(
     eta_hat,
@@ -72,6 +77,7 @@ calc_cwres <- function(
     t_init = 0,
     steady_state_analytic = NULL,
     weight_prior_var = 1,
+    weights = NULL,
     delta = 1e-4
 ) {
   n_obs <- length(y)
@@ -81,6 +87,8 @@ calc_cwres <- function(
     return(list(cwres = numeric(0), vcov = NULL))
   }
 
+  if (is.null(weights)) weights <- rep(1, n_obs)
+
   # Transformed individual predictions at eta_hat
   ipred_transf <- transf(ipred_raw)
 
@@ -149,11 +157,13 @@ calc_cwres <- function(
     )
   }
 
-  # FOCE Hessian: H = (Omega/w)^{-1} + F' * Sigma^{-1} * F
-  # where w = weight_prior_var, matching the scaled omega used during MAP
-  # estimation. 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).
+  # FOCE Hessian: H = (Omega/w)^{-1} + F' * diag(wt/sigma^2) * F
+  # where w = weight_prior_var (scales the prior to match the MAP OFV) and
+  # wt = per-observation weights (must match those used during MAP estimation).
+  # vcov of etas = H^{-1}.
+  # CWRES above uses unit weights (model diagnostic per Hooker et al. 2007);
+  # vcov uses the actual estimation weights so that the Hessian matches the
+  # curvature of the OFV that was minimised.
   # With weight_prior_var <= 0 the prior is effectively flat (LS fit), so
   # the FOCE vcov is not meaningful; skip and let the caller fall back.
   vcov <- NULL
@@ -162,7 +172,7 @@ calc_cwres <- function(
     omega_scaled_inv <- tryCatch(solve(omega_scaled), error = function(e) NULL)
     if (!is.null(omega_scaled_inv)) {
       H_foce <- omega_scaled_inv +
-        t(F_matrix) %*% diag(1 / sigma_diag, nrow = n_obs) %*% F_matrix
+        t(F_matrix) %*% diag(weights / sigma_diag, nrow = n_obs) %*% F_matrix
       vcov <- tryCatch(solve(H_foce), error = function(e) {
         warning("FOCE variance-covariance computation failed.")
         NULL

R/calc_residuals.R

@@ -15,6 +15,8 @@
 #' @param steady_state_analytic steady state settings (or NULL)
 #' @param weight_prior_var prior weight variance scaling factor (default 1).
 #'   Forwarded to `calc_cwres()` to scale omega in the FOCE Hessian.
+#' @param weights numeric vector of per-observation weights. Forwarded to
+#'   `calc_cwres()` so the FOCE Hessian matches the MAP OFV curvature.
 #'
 calc_residuals <- function(
   obj,
@@ -139,7 +141,8 @@ calc_residuals <- function(
         A_init = A_init_individual,
         t_init = t_init,
         steady_state_analytic = steady_state_analytic,
-        weight_prior_var = weight_prior_var
+        weight_prior_var = weight_prior_var,
+        weights = weights
       ),
       error = function(e) {
         warning("CWRES computation failed: ", e$message)

R/get_map_estimates.R

@@ -299,14 +299,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,
@@ -432,11 +436,17 @@ get_map_estimates <- function(
   }
 
   ## 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
+  }
   ## foce_vcov (computed from the FOCE Jacobian) is stored for diagnostics
   ## but not used as the primary vcov source — it has known issues with
   ## observation weights and proportional error that need to be fixed first.
   obj$vcov_full <- get_varcov_matrix(
-    obj$fit$vcov,
+    vcov_source,
     fallback = omega$full
   )
   if(inherits(obj$vcov_full, "matrix")) {

R/mle_wrapper.R

@@ -55,7 +55,10 @@ mle_wrapper <- function(minuslogl,
     try(
       {
         fit$hessian <- numDeriv::hessian(objective_function, fit$par)
-        fit$vcov <- solve(fit$hessian)
+        ## objective_function is on the -2*logLik (NONMEM OFV) scale, so its
+        ## Hessian is twice the observed Fisher information. The asymptotic
+        ## vcov is the inverse of the Fisher information, hence the factor 2.
+        fit$vcov <- 2 * solve(fit$hessian)
       }, 
       silent = FALSE
     )

tests/testthat/test-calc_cwres.R

@@ -0,0 +1,62 @@
+test_that("calc_cwres FOCE vcov uses per-observation weights", {
+  # Minimal 1-compartment IV model with two observations.
+  model <- PKPDsim::new_ode_model(
+    code = "dAdt[0] = -(CL/V) * A[0]",
+    obs = list(cmt = 1, scale = "V"),
+    dose = list(cmt = 1, bioav = 1),
+    parameters = list(CL = 5, V = 50)
+  )
+  regimen <- PKPDsim::new_regimen(
+    amt = 1000, times = 0, t_inf = 1, type = "infusion"
+  )
+
+  parameters_population <- list(CL = 5, V = 50)
+  nonfixed <- c("CL", "V")
+  omega_full <- matrix(c(0.1, 0, 0, 0.1), nrow = 2)
+  eta_hat <- c(0, 0)
+  t_obs <- c(2, 8)
+
+  # Individual predictions at eta_hat (population params, since eta = 0).
+  ipred_raw <- {
+    suppressMessages(
+      PKPDsim::sim_ode(
+        ode = model, parameters = parameters_population, n_ind = 1,
+        regimen = regimen, t_obs = t_obs, only_obs = TRUE, checks = FALSE
+      )$y
+    )
+  }
+  y <- ipred_raw  # observations equal to predictions (residuals irrelevant for vcov)
+
+  args <- list(
+    eta_hat = eta_hat,
+    ipred_raw = ipred_raw,
+    y = y,
+    omega_full = omega_full,
+    error = list(prop = 0.1, add = 0.1),
+    obs_type = c(1, 1),
+    transf = function(x) x,
+    model = model,
+    parameters_population = parameters_population,
+    nonfixed = nonfixed,
+    as_eta = c(),
+    covariates = NULL,
+    regimen = regimen,
+    lagtime = NULL,
+    t_obs = t_obs,
+    obs_type_sim = c(1, 1)
+  )
+
+  res_unit <- do.call(calc_cwres, c(args, list(weights = rep(1, 2))))
+  res_weighted <- do.call(calc_cwres, c(args, list(weights = rep(2, 2))))
+
+  # Higher observation weights increase the likelihood term in the FOCE
+  # Hessian H = Omega^-1 + F' diag(w/sigma^2) F, so vcov = H^-1 shrinks.
+  # Without weights threaded into the Hessian the two vcovs would be identical.
+  expect_false(is.null(res_unit$vcov))
+  expect_false(is.null(res_weighted$vcov))
+  expect_false(isTRUE(all.equal(res_unit$vcov, res_weighted$vcov)))
+  expect_true(det(res_weighted$vcov) < det(res_unit$vcov))
+
+  # CWRES (a model diagnostic) always uses unit weights, so it is unchanged.
+  expect_equal(res_unit$cwres, res_weighted$cwres)
+})

tests/testthat/test-get_map_estimates.R

@@ -199,14 +199,29 @@ test_that("Default MAP fits work and are equal to NONMEM", {
   fit2 <- get_map_estimates(
     model = model, 
     data = obs,
-    parameters = list(CL = 11, V = 90),
+    parameters = list(CL = 5, V = 50),
     regimen = reg, 
-    omega = c(0.1, 0.05, 0.1),
+    omega = c(0.2, 0.05, 0.2),
     error = list(prop = 0.1, add = 10), weights = c(0.25, 1), 
     residuals = T
   )
-  expect_equal(round(fit2$parameters$CL,1), 7.7)
+  expect_equal(round(fit2$parameters$CL,1), 5.9)
   expect_true(diff(abs(fit2$iwres)) < 0.1) # should be small difference. Scaling of residuals is not 100% correct, but seems close enough 
+
+  ## variance-covariance matrix calculated using FOCE vs numDeriv
+  ## are approximately the same
+  fit3 <- get_map_estimates(
+    model = model, 
+    data = obs,
+    parameters = list(CL = 5, V = 50),
+    regimen = reg, 
+    omega = c(0.2, 0.1, 0.2),
+    error = list(prop = 0.1, add = 10), weights = c(0.25, 1), 
+    residuals = FALSE, 
+    skip_hessian = FALSE # this (and residuals=FALSE) forces computation of vcov using numDeriv
+  )
+  ## deviances in vcov calculation methods are < 20%
+  expect_true(all(abs((fit3$vcov_full / fit2$vcov_full)-1) < 0.2))
 })
 
 test_that("allow_obs_before_first_dose works", {
@@ -596,8 +611,11 @@ test_that("FOCE vcov is positive definite even when numDeriv Hessian would fail"
   expect_false(is.null(fit$foce_vcov))
   expect_true(PKPDsim::is_positive_definite(fit$foce_vcov))
   expect_true(PKPDsim::is_positive_definite(fit$vcov_full))
-  skip("Skipping until weights are added to FOCE")
+  ## Since the numDeriv-based vcov fails, should use the 
+  ## FOCE-based vcov. The `vcov_fd` (from numDeriv) should 
+  ## be equal to the input omega as fallback.
   expect_equal(fit$vcov_full, fit$foce_vcov)
+  expect_equal(fit$vcov_fd, PKPDsim::triangle_to_full(omega))
 })
 
 test_that("When vcov from numDeriv is not positive definite, fallback to omega (residuals=FALSE)", {
@@ -923,7 +941,17 @@ test_that("FOCE vcov matches numDeriv Hessian vcov", {
     residuals = FALSE
   )
 
-  expect_equal(fit_foce$vcov, fit_numderiv$vcov, tolerance = 0.01)
+  # Compare on a *relative* scale (ratio ~ 1). The vcov elements (~1e-3) are
+  # smaller than any reasonable absolute tolerance, so expect_equal(a, b,
+  # tolerance = ...) silently degrades to an absolute check and passes even
+  # for a factor-2 scale error. Dividing first forces a true relative check.
+  # The ~7% residual difference is the Gauss-Newton (FOCE) vs full-Hessian
+  # (numDeriv) approximation, not a scaling error.
+  expect_equal(
+    fit_numderiv$vcov / fit_foce$vcov,
+    rep(1, length(fit_foce$vcov)),
+    tolerance = 0.1
+  )
 })
 
 test_that("Floating point precision issues don't raise warning", {