Add etas to results table

CLAUDE.md

@@ -46,6 +46,8 @@ The main entry point is `run_eval()` (`R/run_eval.R`), which orchestrates:
 
 4. **Statistics** (`calculate_stats.R`, `calculate_shrinkage.R`, `calculate_bayesian_impact.R`): Computes RMSE, NRMSE, MAPE, MPE, accuracy, shrinkage, and Bayesian impact metrics.
 
+Failed fits (which surface as `NA` predictions) are detected and reported by `check_failed_fits.R`. `run_eval()` calls it once after the per-subject loop to emit a single warning; `calculate_stats()` also calls it (controllable via its `warn` argument, which `run_eval()` sets to `FALSE` to avoid a duplicate warning).
+
 ### Output Structure
 
 `run_eval()` returns a list with class `"mipdeval_results"`:

NAMESPACE

@@ -22,10 +22,14 @@ export(install_default_literature_model)
 export(is_accurate)
 export(is_accurate_abs)
 export(is_accurate_rel)
+export(mape)
+export(mpe)
 export(new_ode_model)
+export(nrmse)
 export(parse_psn_proseval_results)
 export(reldiff_psn_execute_results)
 export(reldiff_psn_proseval_results)
+export(rmse)
 export(run_eval)
 export(stats_summ_options)
 export(vpc_options)

R/calculate_stats.R

@@ -7,6 +7,8 @@
 #'   providing an absolute or relative error margin. The cutoff is exclusive of
 #'   the error margin. When `NULL` (the default), accuracy will not be
 #'   calculated and will return `NA` instead.
+#' @param warn should a warning be emitted when failed fits (NA predictions) are
+#'   detected?
 #'
 #' @returns tibble
 #'
@@ -15,21 +17,14 @@ calculate_stats <- function(
     .res,
     rounding = 3,
     acc_error_abs = NULL,
-    acc_error_rel = NULL
+    acc_error_rel = NULL,
+    warn = TRUE
 ) {
   if(inherits(.res, "mipdeval_results")) {
     .res <- .res$results
   }
-  ## Check for errors during fits / predictions
-  errors <- dplyr::filter(
-    .res,
-    is.na(.data$pred) |
-    (is.na(.data$map_ipred) & !.data$apriori) |
-    is.na(.data$iter_ipred)
-  )
-  if(nrow(errors) > 0) {
-    cli::cli_warn("Errors were encountered in {nrow(errors)} out of {nrow(.res)} evaluated predictions. The problems occurred in patient(s) {unique(errors$id)}.")
-  }
+  ## Warn about any failed fits / predictions (NA), unless the caller opts out.
+  if (isTRUE(warn)) check_failed_fits(.res)
   out <- .res |>
     tidyr::pivot_longer(
       cols = c("pred", "map_ipred", "iter_ipred"), names_to = "type"

R/check_failed_fits.R

@@ -0,0 +1,25 @@
+#' Check for failed fits / predictions and (optionally) warn
+#'
+#' Detects predictions that came back as `NA`, which indicates the underlying
+#' MAP Bayesian fit failed, and---when `warn = TRUE`---emits a warning
+#' summarising how many predictions failed and in which subjects.
+#'
+#' @param .res output object (`mipdeval_results`) from [run_eval()], or a
+#'   `data.frame` with raw results.
+#'
+#' @returns invisibly, a `data.frame` of the rows with failed predictions.
+check_failed_fits <- function(.res) {
+  if(inherits(.res, "mipdeval_results")) {
+    .res <- .res$results
+  }
+  errors <- dplyr::filter(
+    .res,
+    is.na(.data$pred) |
+    (is.na(.data$map_ipred) & !.data$apriori) |
+    is.na(.data$iter_ipred)
+  )
+  if(nrow(errors) > 0) {
+    cli::cli_warn("Errors were encountered in {nrow(errors)} out of {nrow(.res)} evaluated predictions. The problems occurred in patient(s) {unique(errors$id)}.")
+  }
+  invisible(errors)
+}

R/misc.R

@@ -51,6 +51,7 @@ vec_assert_or_null <- function(
 #' @param pred predictions vector
 #'
 #' @returns A numeric vector
+#' @export
 rmse <- function (obs, pred) {
   res_sq <- (pred - obs)^2
   sqrt(mean(res_sq, na.rm = TRUE))
@@ -62,7 +63,7 @@ rmse <- function (obs, pred) {
 #' @param pred predictions vector
 #'
 #' @returns A numeric vector
-#'
+#' @export
 nrmse <- function (obs, pred) {
   res_sq <- (pred - obs)^2
   rmse <- sqrt(mean(res_sq, na.rm = TRUE))
@@ -74,6 +75,7 @@ nrmse <- function (obs, pred) {
 #' @inheritParams rmse
 #'
 #' @returns A numeric vector
+#' @export
 mape <- function (obs, pred) {
   mean(abs((obs - pred))/obs, na.rm = TRUE)
 }
@@ -83,6 +85,7 @@ mape <- function (obs, pred) {
 #' @inheritParams rmse
 #'
 #' @returns A numeric vector
+#' @export
 mpe <- function (obs, pred) {
   mean((obs - pred)/obs, na.rm = TRUE)
 }

R/run_eval.R

@@ -45,7 +45,65 @@
 #'   debugging the package it is useful to have it off, since progress bar
 #'   handlers obscure R output.
 #'
-#' @returns A named list of data frames.
+#' @details
+#' `run_eval()` evaluates predictive performance the way a model would be used
+#' for model-informed precision dosing (MIPD): for each subject it walks through
+#' the observations in time order and, at each step, refits the individual using
+#' only the data available *up to that point* to predict the next observation(s).
+#' This produces an a priori (population) prediction followed by progressively
+#' more informed a posteriori (forecasting) predictions, mirroring the iterative
+#' flow of PsN's `proseval` tool.
+#'
+#' Because of this design the `results` tibble reports quantities at three levels
+#' of individualization, which are easy to confuse. The predictions:
+#'
+#' - `pred`: the **population** prediction (typical individual, using no
+#'   individual-level data)---the a priori prediction.
+#' - `iter_ipred`: the **iterative ("forecasting")** individual prediction, using
+#'   only the data available up to that point---the a posteriori prediction. On
+#'   the a priori rows it falls back to `pred`.
+#' - `map_ipred`: the **full-data MAP** individual prediction, from a single
+#'   retrospective fit on all of a subject's data at once.
+#'
+#' The random-effect (eta) estimates follow the same structure:
+#'
+#' - The `eta<nn>` columns are the iterative ("forecasting") empirical Bayes
+#'   estimates that pair with `iter_ipred`; they evolve down the rows and are 0
+#'   on the a priori (population) rows.
+#' - The `map_eta<nn>` columns are the full-data MAP empirical Bayes estimates
+#'   that pair with `map_ipred`. They are constant per subject and appear on
+#'   every row, including the a priori row. These are the equivalent of the etas
+#'   reported in a NONMEM output table, and are what you want for an empirical
+#'   eta-distribution plot (one value per subject).
+#'
+#' @returns An `mipdeval_results` object, which is a named list with the following
+#'   elements:
+#'
+#'   - `results`: A tibble with one row per (iterative) prediction, holding the
+#'     identifiers (`id`, `_iteration`, `_grouper`, `t`), the observation (`dv`),
+#'     population predictions and residuals (`pred`, `res`, `wres`, `cwres`),
+#'     individual predictions and residuals (`iter_ipred`, `map_ipred`, `ires`,
+#'     `iwres`), the objective function value (`ofv`), the weighted sum-of-squares
+#'     (`ss_w`), an `apriori` flag, one column per model parameter, and two
+#'     families of eta (random-effect) columns: the iterative `eta<nn>` and the
+#'     full-data `map_eta<nn>` estimates (see Details).
+#'   - `mod_obj`: The parsed model object (see [parse_model()]): a named list of
+#'     model information, including `model`, `parameters`, `omega`, `ruv`,
+#'     `fixed`, `bins`, and `kappa`.
+#'   - `data`: The input data after reading and validation (see
+#'     [read_input_data()] and [check_input_data()]), as a data frame of the
+#'     NONMEM-style records used in the analysis.
+#'   - `sim`: Simulated data used for the visual predictive check (VPC) and NPDE,
+#'     or `NULL` when simulations are skipped (`vpc_options(skip = TRUE)`).
+#'   - `stats_summ`: A tibble summarising forecasting performance statistics (see
+#'     [calculate_stats()]).
+#'   - `shrinkage`: A tibble of eta-shrinkage per iteration (see
+#'     [calculate_shrinkage()]).
+#'   - `bayesian_impact`: A tibble of Bayesian-impact values based on RMSE and
+#'     MAPE (see [calculate_bayesian_impact()]).
+#'
+#'   `stats_summ`, `shrinkage`, and `bayesian_impact` are `NULL` when no
+#'   predictions are produced (e.g. `vpc_options(vpc_only = TRUE)`).
 #'
 #' @export
 run_eval <- function(
@@ -178,12 +236,14 @@ run_eval <- function(
 
   # res is NULL when vpc_options(..., vpc_only = TRUE).
   if (!is.null(res)) {
+    check_failed_fits(res)
     if(verbose) cli::cli_progress_step("Calculating forecasting statistics")
     out$stats_summ <- calculate_stats(
       out,
       rounding = .stats_summ_options$rounding,
       acc_error_abs = .stats_summ_options$acc_error_abs,
-      acc_error_rel = .stats_summ_options$acc_error_rel
+      acc_error_rel = .stats_summ_options$acc_error_rel,
+      warn = FALSE # Avoid a duplicate warning from check_failed_fits()
     )
     out$shrinkage <- calculate_shrinkage(out)
     out$bayesian_impact <- calculate_bayesian_impact(out)

R/run_eval_core.R

@@ -23,6 +23,7 @@ run_eval_core <- function(
   obs_data <- data$observations
   comb <- data.frame()
   fit_pars <- data.frame()
+  eta_names <- character(0) # captured from first successful fit below
   iterations <- unique(obs_data[["_grouper"]])
 
   for(i in seq_along(iterations)) {
@@ -95,6 +96,9 @@ run_eval_core <- function(
       )
       par_dummy <- as.data.frame(mod_upd$parameters)
       par_dummy[, 1:ncol(par_dummy)] <- NA
+      ## Also emit NA eta columns if a previous fit told us their names, so a
+      ## failed fit has the same shape as a successful one.
+      par_dummy[, eta_names] <- NA_real_
       fit_pars <- dplyr::mutate(as.data.frame(par_dummy), id = obs_data$id[1])
     } else {
       ## Data frame with predictive data
@@ -114,8 +118,9 @@ run_eval_core <- function(
         `_iteration` = iterations[i],
         `_grouper` = obs_data$`_grouper`
       )
-      ## Add parameter estimates
-      fit_pars <- dplyr::mutate(as.data.frame(fit$parameters), id = obs_data$id[1])
+      ## Add parameter estimates and etas:
+      eta_names <- names(fit$fit$par)
+      fit_pars <- dplyr::mutate(as.data.frame(c(fit$parameters, fit$fit$par)), id = obs_data$id[1])
     }
 
     comb <- dplyr::bind_rows(
@@ -163,6 +168,19 @@ run_eval_core <- function(
     map_pred_data <- pred_data
   }
 
+  ## Full-data MAP etas (empirical Bayes estimates), constant per subject and
+  ## named to parallel `map_ipred`. These differ from the iterative `eta_names`
+  ## columns, which only use the data available up to each forecast. NA if the
+  ## MAP fit failed.
+  map_etas <- if(!is.null(fit_map$fit$par)) {
+    as.list(fit_map$fit$par[eta_names])
+  } else {
+    as.list(rep(NA_real_, length(eta_names)))
+  }
+  ## guard against the no-successful-fit case where `eta_names` is empty.
+  map_eta_names <- if(length(eta_names) > 0) paste0("map_", eta_names) else character(0)
+  names(map_etas) <- map_eta_names
+
   ## pre-pend population predictions for the first observation
   # TODO: Refactor this logic into a function or functions, e.g., the first
   # argument to bind_rows() could be refactored into `get_apriori_data()`.
@@ -173,7 +191,8 @@ run_eval_core <- function(
         `_iteration` = 0,
         ipred = .data$pred,
         ofv = NA,
-        ss_w = NA
+        ss_w = NA,
+        dplyr::across(dplyr::all_of(eta_names), ~ 0) # population etas are 0
       ) |> # set to population parameters, not individual estimates
         dplyr::select(-!!names(mod_obj$parameters)) |>
         dplyr::left_join(
@@ -186,12 +205,13 @@ run_eval_core <- function(
     dplyr::mutate(
       iter_ipred = .data$ipred,
       map_ipred = map_pred_data$ipred, # ipred from full retrospective MAP
-      apriori = (.data$`_iteration` == 0)
+      apriori = (.data$`_iteration` == 0),
+      !!!map_etas # full-data MAP etas, constant per subject
     ) |>
     dplyr::select(
       "id", "_iteration", "_grouper", "t", "dv", "pred", "res", "wres", "cwres",
       "map_ipred", "ofv", "ss_w", "iter_ipred", "ires", "iwres", "apriori",
-      !!names(mod_obj$parameters)
+      !!names(mod_obj$parameters), !!eta_names, !!map_eta_names
     )
 
   out