From 5e9063e5fa0f985b46eae450a7551e2cd276024e Mon Sep 17 00:00:00 2001
From: Ron Keizer <ron@insight-rx.com>
Date: Tue, 23 Sep 2025 04:38:49 +0000
Subject: [PATCH 01/22] add more stats

---
 NAMESPACE                              |  1 +
 R/calculate_bayesian_impact.R          | 20 +++++++
 R/calculate_shrinkage.R                | 75 ++++++++++++++++++++++++++
 R/run_eval.R                           | 19 ++++++-
 R/run_eval_core.R                      |  6 +--
 man/calc_eta.Rd                        | 17 ++++++
 man/calculate_bayesian_impact.Rd       | 15 ++++++
 man/calculate_shrinkage.Rd             | 31 +++++++++++
 man/get_omega_for_parameters.Rd        | 14 +++++
 man/nm_busulfan.Rd                     |  2 +-
 man/run_eval.Rd                        |  5 +-
 tests/testthat/test-run_eval.R         |  4 +-
 tests/testthat/test-run_eval_grouped.R |  4 +-
 13 files changed, 202 insertions(+), 11 deletions(-)
 create mode 100644 R/calculate_bayesian_impact.R
 create mode 100644 R/calculate_shrinkage.R
 create mode 100644 man/calc_eta.Rd
 create mode 100644 man/calculate_bayesian_impact.Rd
 create mode 100644 man/calculate_shrinkage.Rd
 create mode 100644 man/get_omega_for_parameters.Rd

diff --git a/NAMESPACE b/NAMESPACE
index e39af50..44104e9 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -4,6 +4,7 @@ S3method(parse_model,PKPDsim)
 S3method(parse_model,character)
 S3method(parse_model,default)
 export(add_grouping_column)
+export(calculate_shrinkage)
 export(compare_psn_execute_results)
 export(compare_psn_proseval_results)
 export(group_by_dose)
diff --git a/R/calculate_bayesian_impact.R b/R/calculate_bayesian_impact.R
new file mode 100644
index 0000000..94c6b22
--- /dev/null
+++ b/R/calculate_bayesian_impact.R
@@ -0,0 +1,20 @@
+#' This function calculates the impact of using Bayesian updating, compared
+#' to using population estimates
+#'
+#' Bayesian impact (BI), as defined here, is the gain in accuracy that is
+#' obtained when using Bayesian updating, compared to using a priori
+#' (population) parameter estimates. The BI is computed based on accuracy
+#' measures (RMSE and MAPE), but not on bias estimates (MPE).
+#'
+calculate_bayesian_impact <- function(
+  stats_summ
+) {
+  stats_summ |>
+    dplyr::filter(!apriori) |>
+    dplyr::select(type, apriori, rmse, mape) |>
+    tidyr::pivot_wider(names_from = "type", values_from = c("rmse", "mape")) |>
+    dplyr::summarise(
+      bi_rmse = 100 * (1 - (rmse_iter_ipred / rmse_pred)),
+      bi_mape = 100 * (1 - (mape_iter_ipred / mape_pred))
+    )
+}
diff --git a/R/calculate_shrinkage.R b/R/calculate_shrinkage.R
new file mode 100644
index 0000000..dd27db8
--- /dev/null
+++ b/R/calculate_shrinkage.R
@@ -0,0 +1,75 @@
+#' Function to calculate eta-shrinkage
+#'
+#' Shrinkage for population PK model was first defined in this paper by
+#' Savic and Karlsson: https://pmc.ncbi.nlm.nih.gov/articles/PMC2758126/.
+#' It is a measure of how much information is available to update individual
+#' estimates away from the population value. In principle, if there is no
+#' information at all (i.e. in the case of population estimates only),
+#' shrinkage is 100%. In the case of fully informed Bayesian estimates
+#' (unlikely to be achieved in practice), shrinkage is 0%. In most practical
+#' scenarios with limited sampling, shrinkage will be between 10-30%. When
+#' shrinkage is higher than 50% or so, one could argue there is limited benefit
+#' of sampling at all. So the shrinkage reported in this package can be used
+#' to evaluate whether the sampling was efficient, and how shrinkage could be
+#' reduced with additional sampling.
+#'
+#' @inheritParams run_eval_core
+#' @param res_df data.frame or tibble with raw results (including parameters)
+#' from a run with `mipdeval::run_eval()`.
+#'
+#' @returns tibble
+#' @export
+calculate_shrinkage <- function(
+  res_df,
+  mod_obj
+) {
+  om <- get_omega_for_parameters(mod_obj)
+  shr_df <- res_df |>
+    dplyr::mutate(across(
+      .cols = names(om),
+      .fns  = ~ calc_eta(.x, dplyr::cur_column(), mod_obj$parameters),
+      .names = "eta_{.col}"
+    )) |>
+    dplyr::group_by(`_iteration`) |>
+    dplyr::summarise(across(
+      .cols = paste0("eta_", names(om)),
+      .fns = ~ calc_shrinkage(.x, dplyr::cur_column(), om),
+      .names = "shr_{.col}"
+    ))
+  names(shr_df) <- gsub("^shr_eta_", "", names(shr_df))
+  shr_df
+}
+
+#' Calculate the "eta"-value for a parameters, assuming an exponential
+#' shape for IIV, i.e. `PAR = TV_PAR * EXP(ETA(n))`, and
+#' additive shape for IOV (kappa's), i.e. `PAR = TV_PAR + ETA(n)` or when
+#' final parameter estimate is 0.
+#'
+calc_eta <- function(ind, par_name, parameters) {
+  pop <- parameters[[par_name]]
+  if(stringr::str_detect("^kappa_", par_name) || pop == 0) { ## PAR = TV_PAR + ETA()
+    return(ind - pop)
+  } else {
+    return(log(ind/pop)) ## PAR = TV_PAR + ETA()
+  }
+}
+
+calc_shrinkage <- function(eta, par_name, om) {
+  100 * (1 - (sd(eta) / sqrt(om[[gsub("^eta_", "", par_name)]])))
+}
+
+#' Return all diagonal om^2 elements for each non-fixed parameter, as a list
+#'
+#' @inheritParams run_core
+#' @returns list
+get_omega_for_parameters <- function(mod_obj) {
+  pars <- names(mod_obj$parameters)
+  pars <- pars[! (pars %in% mod_obj$fixed)]
+  om <- as.list(
+    diag(
+      PKPDsim::triangle_to_full(mod_obj$omega)
+    )[1:length(pars)]
+  )
+  names(om) <- pars
+  om
+}
diff --git a/R/run_eval.R b/R/run_eval.R
index 91cc821..5414c5c 100644
--- a/R/run_eval.R
+++ b/R/run_eval.R
@@ -58,7 +58,8 @@ run_eval <- function(
   censor_covariates = TRUE,
   incremental = FALSE,
   threads = 1,
-  progress = TRUE
+  progress = TRUE,
+  verbose = TRUE
 ) {
   # TODO: Refactor to S3 method like in parse_model() to make required and
   # optional arguments clearer for the different types of `model` inputs.
@@ -74,6 +75,7 @@ run_eval <- function(
   )
 
   ## 1. read NONMEM data from file or data.frame. Do some simple checks
+  if(verbose) cli::cli_alert_info("Reading and parsing input data")
   input_data <- read_input_data(data) |>
     check_input_data(
       dictionary = dictionary
@@ -103,6 +105,7 @@ run_eval <- function(
   }
 
   ## 3. run the core function on each individual-level dataset in the list
+  if(verbose) cli::cli_alert_info("Running forecasts for subjects in dataset")
   if(threads > 1) {
     # TODO: consider using purrr::in_parallel() in the future when it's stable.
     future::plan(future::multisession, workers = threads)
@@ -126,6 +129,7 @@ run_eval <- function(
   }
 
   ## 4. Combine results and basic stats into return object
+  if(verbose) cli::cli_alert_info("Calculating forecasting statistics")
   res_df <- dplyr::bind_rows(res)
   stats_summ <- res_df |>
     tidyr::pivot_longer(
@@ -137,13 +141,24 @@ run_eval <- function(
       mpe = mpe(.data$dv, .data$value),
       mape = mape(.data$dv, .data$value)
     )
+  shrinkage_summ <- calculate_shrinkage(
+    res_df,
+    mod_obj = mod_obj
+  )
+  bayesian_impact <- calculate_bayesian_impact(
+    stats_summ
+  )
+
   out <- list(
     results = tibble::as_tibble(res_df),
-    stats = tibble::as_tibble(stats_summ)
+    stats = tibble::as_tibble(stats_summ),
+    shrinkage = shrinkage_summ,
+    bayesian_impact = bayesian_impact
   )
 
   # TODO: Turn out into an S3 class, so we can give it methods like print(), etc.
 
   ## 5. Return results
+  if(verbose) cli::cli_alert_info("Done")
   out
 }
diff --git a/R/run_eval_core.R b/R/run_eval_core.R
index 7e82fb0..0e2bd43 100644
--- a/R/run_eval_core.R
+++ b/R/run_eval_core.R
@@ -89,9 +89,9 @@ run_eval_core <- function(
         `_iteration` = 0,
         ipred = .data$pred
       ) |> # set to population parameters, not individual estimates
-        dplyr::select(-!!names(parameters)) |>
+        dplyr::select(-!!names(mod_obj$parameters)) |>
         dplyr::left_join(
-          dplyr::mutate(as.data.frame(parameters), id = obs_data$id[1]),
+          dplyr::mutate(as.data.frame(mod_obj$parameters), id = obs_data$id[1]),
           by = "id"
         ),
     comb
@@ -105,7 +105,7 @@ run_eval_core <- function(
     dplyr::select(
       "id", "_iteration", "_grouper", "t", "dv", "pred", "map_ipred",
       "iter_ipred", "apriori",
-      !!names(parameters)
+      !!names(mod_obj$parameters)
     )
 
   out
diff --git a/man/calc_eta.Rd b/man/calc_eta.Rd
new file mode 100644
index 0000000..54caa07
--- /dev/null
+++ b/man/calc_eta.Rd
@@ -0,0 +1,17 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/calculate_shrinkage.R
+\name{calc_eta}
+\alias{calc_eta}
+\title{Calculate the "eta"-value for a parameters, assuming an exponential
+shape for IIV, i.e. \code{PAR = TV_PAR * EXP(ETA(n))}, and
+additive shape for IOV (kappa's), i.e. \code{PAR = TV_PAR + ETA(n)} or when
+final parameter estimate is 0.}
+\usage{
+calc_eta(ind, par_name, parameters)
+}
+\description{
+Calculate the "eta"-value for a parameters, assuming an exponential
+shape for IIV, i.e. \code{PAR = TV_PAR * EXP(ETA(n))}, and
+additive shape for IOV (kappa's), i.e. \code{PAR = TV_PAR + ETA(n)} or when
+final parameter estimate is 0.
+}
diff --git a/man/calculate_bayesian_impact.Rd b/man/calculate_bayesian_impact.Rd
new file mode 100644
index 0000000..cec8568
--- /dev/null
+++ b/man/calculate_bayesian_impact.Rd
@@ -0,0 +1,15 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/calculate_bayesian_impact.R
+\name{calculate_bayesian_impact}
+\alias{calculate_bayesian_impact}
+\title{This function calculates the impact of using Bayesian updating, compared
+to using population estimates}
+\usage{
+calculate_bayesian_impact(stats_summ)
+}
+\description{
+Bayesian impact (BI), as defined here, is the gain in accuracy that is
+obtained when using Bayesian updating, compared to using a priori
+(population) parameter estimates. The BI is computed based on accuracy
+measures (RMSE and MAPE), but not on bias estimates (MPE).
+}
diff --git a/man/calculate_shrinkage.Rd b/man/calculate_shrinkage.Rd
new file mode 100644
index 0000000..9605df3
--- /dev/null
+++ b/man/calculate_shrinkage.Rd
@@ -0,0 +1,31 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/calculate_shrinkage.R
+\name{calculate_shrinkage}
+\alias{calculate_shrinkage}
+\title{Function to calculate eta-shrinkage}
+\usage{
+calculate_shrinkage(res_df, mod_obj)
+}
+\arguments{
+\item{res_df}{data.frame or tibble with raw results (including parameters)
+from a run with \code{mipdeval::run_eval()}.}
+
+\item{mod_obj}{list object with model information}
+}
+\value{
+tibble
+}
+\description{
+Shrinkage for population PK model was first defined in this paper by
+Savic and Karlsson: https://pmc.ncbi.nlm.nih.gov/articles/PMC2758126/.
+It is a measure of how much information is available to update individual
+estimates away from the population value. In principle, if there is no
+information at all (i.e. in the case of population estimates only),
+shrinkage is 100\%. In the case of fully informed Bayesian estimates
+(unlikely to be achieved in practice), shrinkage is 0\%. In most practical
+scenarios with limited sampling, shrinkage will be between 10-30\%. When
+shrinkage is higher than 50\% or so, one could argue there is limited benefit
+of sampling at all. So the shrinkage reported in this package can be used
+to evaluate whether the sampling was efficient, and how shrinkage could be
+reduced with additional sampling.
+}
diff --git a/man/get_omega_for_parameters.Rd b/man/get_omega_for_parameters.Rd
new file mode 100644
index 0000000..e85bada
--- /dev/null
+++ b/man/get_omega_for_parameters.Rd
@@ -0,0 +1,14 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/calculate_shrinkage.R
+\name{get_omega_for_parameters}
+\alias{get_omega_for_parameters}
+\title{Return all diagonal om^2 elements for each non-fixed parameter, as a list}
+\usage{
+get_omega_for_parameters(mod_obj)
+}
+\value{
+list
+}
+\description{
+Return all diagonal om^2 elements for each non-fixed parameter, as a list
+}
diff --git a/man/nm_busulfan.Rd b/man/nm_busulfan.Rd
index a041f12..71ccbd9 100644
--- a/man/nm_busulfan.Rd
+++ b/man/nm_busulfan.Rd
@@ -6,7 +6,7 @@
 \title{Busulfan data}
 \format{
 A NONMEM-style data frame with 240 rows and
-13 columns:
+14 columns:
 
 \describe{
 \item{\code{ID}}{Patient ID.}
diff --git a/man/run_eval.Rd b/man/run_eval.Rd
index 51f3dd3..d462754 100644
--- a/man/run_eval.Rd
+++ b/man/run_eval.Rd
@@ -20,7 +20,8 @@ run_eval(
   censor_covariates = TRUE,
   incremental = FALSE,
   threads = 1,
-  progress = TRUE
+  progress = TRUE,
+  verbose = TRUE
 )
 }
 \arguments{
@@ -85,6 +86,8 @@ i.e. no parallel execution}
 \item{progress}{should a progress bar be shown? Default is \code{TRUE}, but when
 debugging the package it is useful to have it off, since progress bar
 handlers obscure R output.}
+
+\item{verbose}{show more output}
 }
 \value{
 A named list of data frames.
diff --git a/tests/testthat/test-run_eval.R b/tests/testthat/test-run_eval.R
index 80b4aff..3c321b1 100644
--- a/tests/testthat/test-run_eval.R
+++ b/tests/testthat/test-run_eval.R
@@ -11,11 +11,11 @@ test_that("Basic run with vanco data + model works", {
     ruv = mod_obj$ruv,
     fixed = mod_obj$fixed,
     censor_covariates = FALSE, # shouldn't matter, since no time-varying covs
-    progress = FALSE
+    progress = TRUE
   )
 
   ## Expected structure:
-  expect_equal(names(res), c("results", "stats"))
+  expect_equal(names(res), c("results", "stats", "shrinkage", "bayesian_impact"))
   expect_s3_class(res$results, c("tbl_df", "tbl", "data.frame"))
   expect_s3_class(res$stats, c("tbl_df", "tbl", "data.frame"))
 
diff --git a/tests/testthat/test-run_eval_grouped.R b/tests/testthat/test-run_eval_grouped.R
index 646d2fa..4b29090 100644
--- a/tests/testthat/test-run_eval_grouped.R
+++ b/tests/testthat/test-run_eval_grouped.R
@@ -10,10 +10,10 @@ test_that("Grouped run for busulfan/shukla works", {
     model = "pkbusulfanshukla",
     data = nm_busulfan,
     group = "group",
-    progress = FALSE,
+    progress = TRUE,
     ids = c(1:n_ids)
   )
-  expect_equal(names(res), c("results", "stats"))
+  expect_equal(names(res), c("results", "stats", "shrinkage", "bayesian_impact"))
 
   ## Reference results from PsN proseval:
   raw_proseval <- read.csv(

From 9bb7eb7e03c4f01e3bbc8f6b1e909a5b653c25c8 Mon Sep 17 00:00:00 2001
From: Ron Keizer <ron@insight-rx.com>
Date: Tue, 23 Sep 2025 12:22:17 -0700
Subject: [PATCH 02/22] add tests + update docs

---
 R/calculate_shrinkage.R                       |  11 +-
 man/calculate_shrinkage.Rd                    |   6 +-
 man/get_omega_for_parameters.Rd               |   3 +
 man/nm_busulfan.Rd                            |   2 +-
 .../testthat/test-calculate_bayesian_impact.R | 110 +++++++++
 tests/testthat/test-calculate_shrinkage.R     | 208 ++++++++++++++++++
 tests/testthat/test-run_eval.R                |   4 +-
 7 files changed, 333 insertions(+), 11 deletions(-)
 create mode 100644 tests/testthat/test-calculate_bayesian_impact.R
 create mode 100644 tests/testthat/test-calculate_shrinkage.R

diff --git a/R/calculate_shrinkage.R b/R/calculate_shrinkage.R
index dd27db8..382f231 100644
--- a/R/calculate_shrinkage.R
+++ b/R/calculate_shrinkage.R
@@ -1,13 +1,13 @@
 #' Function to calculate eta-shrinkage
 #'
-#' Shrinkage for population PK model was first defined in this paper by
+#' Shrinkage for population PK models was first defined in this paper by
 #' Savic and Karlsson: https://pmc.ncbi.nlm.nih.gov/articles/PMC2758126/.
 #' It is a measure of how much information is available to update individual
 #' estimates away from the population value. In principle, if there is no
 #' information at all (i.e. in the case of population estimates only),
-#' shrinkage is 100%. In the case of fully informed Bayesian estimates
+#' shrinkage will be 100%. In the case of fully informed Bayesian estimates
 #' (unlikely to be achieved in practice), shrinkage is 0%. In most practical
-#' scenarios with limited sampling, shrinkage will be between 10-30%. When
+#' scenarios with limited sampling, shrinkage will be between 10-40%. When
 #' shrinkage is higher than 50% or so, one could argue there is limited benefit
 #' of sampling at all. So the shrinkage reported in this package can be used
 #' to evaluate whether the sampling was efficient, and how shrinkage could be
@@ -47,7 +47,7 @@ calculate_shrinkage <- function(
 #'
 calc_eta <- function(ind, par_name, parameters) {
   pop <- parameters[[par_name]]
-  if(stringr::str_detect("^kappa_", par_name) || pop == 0) { ## PAR = TV_PAR + ETA()
+  if(stringr::str_detect(par_name, "^kappa_") || pop == 0) { ## PAR = TV_PAR + ETA()
     return(ind - pop)
   } else {
     return(log(ind/pop)) ## PAR = TV_PAR + ETA()
@@ -60,7 +60,8 @@ calc_shrinkage <- function(eta, par_name, om) {
 
 #' Return all diagonal om^2 elements for each non-fixed parameter, as a list
 #'
-#' @inheritParams run_core
+#' @inheritParams run_eval_core
+#' 
 #' @returns list
 get_omega_for_parameters <- function(mod_obj) {
   pars <- names(mod_obj$parameters)
diff --git a/man/calculate_shrinkage.Rd b/man/calculate_shrinkage.Rd
index 9605df3..03d289c 100644
--- a/man/calculate_shrinkage.Rd
+++ b/man/calculate_shrinkage.Rd
@@ -16,14 +16,14 @@ from a run with \code{mipdeval::run_eval()}.}
 tibble
 }
 \description{
-Shrinkage for population PK model was first defined in this paper by
+Shrinkage for population PK models was first defined in this paper by
 Savic and Karlsson: https://pmc.ncbi.nlm.nih.gov/articles/PMC2758126/.
 It is a measure of how much information is available to update individual
 estimates away from the population value. In principle, if there is no
 information at all (i.e. in the case of population estimates only),
-shrinkage is 100\%. In the case of fully informed Bayesian estimates
+shrinkage will be 100\%. In the case of fully informed Bayesian estimates
 (unlikely to be achieved in practice), shrinkage is 0\%. In most practical
-scenarios with limited sampling, shrinkage will be between 10-30\%. When
+scenarios with limited sampling, shrinkage will be between 10-40\%. When
 shrinkage is higher than 50\% or so, one could argue there is limited benefit
 of sampling at all. So the shrinkage reported in this package can be used
 to evaluate whether the sampling was efficient, and how shrinkage could be
diff --git a/man/get_omega_for_parameters.Rd b/man/get_omega_for_parameters.Rd
index e85bada..6adbdf9 100644
--- a/man/get_omega_for_parameters.Rd
+++ b/man/get_omega_for_parameters.Rd
@@ -6,6 +6,9 @@
 \usage{
 get_omega_for_parameters(mod_obj)
 }
+\arguments{
+\item{mod_obj}{list object with model information}
+}
 \value{
 list
 }
diff --git a/man/nm_busulfan.Rd b/man/nm_busulfan.Rd
index 71ccbd9..a041f12 100644
--- a/man/nm_busulfan.Rd
+++ b/man/nm_busulfan.Rd
@@ -6,7 +6,7 @@
 \title{Busulfan data}
 \format{
 A NONMEM-style data frame with 240 rows and
-14 columns:
+13 columns:
 
 \describe{
 \item{\code{ID}}{Patient ID.}
diff --git a/tests/testthat/test-calculate_bayesian_impact.R b/tests/testthat/test-calculate_bayesian_impact.R
new file mode 100644
index 0000000..c1ef1d8
--- /dev/null
+++ b/tests/testthat/test-calculate_bayesian_impact.R
@@ -0,0 +1,110 @@
+test_that("calculate_bayesian_impact returns correct structure", {
+  # Create test data
+  test_data <- data.frame(
+    type = c("pred", "iter_ipred", "pred", "iter_ipred"),
+    apriori = c(FALSE, FALSE, TRUE, TRUE),
+    rmse = c(10, 5, 15, 8),
+    mape = c(20, 10, 25, 12)
+  )
+  
+  result <- calculate_bayesian_impact(test_data)
+  
+  expect_s3_class(result, "data.frame")
+  expect_named(result, c("bi_rmse", "bi_mape"))
+  expect_equal(nrow(result), 1)
+})
+
+test_that("calculate_bayesian_impact calculates correct values", {
+  # Test data with known expected results
+  test_data <- data.frame(
+    type = c("pred", "iter_ipred"),
+    apriori = c(FALSE, FALSE),
+    rmse = c(10, 5),  # iter_ipred is 50% better
+    mape = c(20, 10)  # iter_ipred is 50% better
+  )
+  
+  result <- calculate_bayesian_impact(test_data)
+  
+  # BI = 100 * (1 - (iter_ipred / pred))
+  # For RMSE: 100 * (1 - (5/10)) = 50%
+  # For MAPE: 100 * (1 - (10/20)) = 50%
+  expect_equal(result$bi_rmse, 50)
+  expect_equal(result$bi_mape, 50)
+})
+
+test_that("calculate_bayesian_impact handles negative impact correctly", {
+  # Test case where Bayesian updating performs worse
+  test_data <- data.frame(
+    type = c("pred", "iter_ipred"),
+    apriori = c(FALSE, FALSE),
+    rmse = c(5, 10),   # iter_ipred is worse
+    mape = c(10, 20)   # iter_ipred is worse
+  )
+  
+  result <- calculate_bayesian_impact(test_data)
+  
+  # BI = 100 * (1 - (10/5)) = -100%
+  expect_equal(result$bi_rmse, -100)
+  expect_equal(result$bi_mape, -100)
+})
+
+test_that("calculate_bayesian_impact filters apriori correctly", {
+  # Include apriori = TRUE rows that should be filtered out
+  test_data <- data.frame(
+    type = c("pred", "iter_ipred", "pred", "iter_ipred"),
+    apriori = c(TRUE, TRUE, FALSE, FALSE),
+    rmse = c(100, 50, 10, 5),
+    mape = c(200, 100, 20, 10)
+  )
+  
+  result <- calculate_bayesian_impact(test_data)
+  
+  # Should only use the apriori = FALSE rows (10, 5 and 20, 10)
+  expect_equal(result$bi_rmse, 50)
+  expect_equal(result$bi_mape, 50)
+})
+
+test_that("calculate_bayesian_impact handles edge cases", {
+  # Test with identical values (no improvement)
+  test_data <- data.frame(
+    type = c("pred", "iter_ipred"),
+    apriori = c(FALSE, FALSE),
+    rmse = c(10, 10),
+    mape = c(20, 20)
+  )
+  
+  result <- calculate_bayesian_impact(test_data)
+  
+  expect_equal(result$bi_rmse, 0)
+  expect_equal(result$bi_mape, 0)
+})
+
+test_that("calculate_bayesian_impact handles division by zero", {
+  # Test with pred values of zero
+  test_data <- data.frame(
+    type = c("pred", "iter_ipred"),
+    apriori = c(FALSE, FALSE),
+    rmse = c(0, 5),
+    mape = c(0, 10)
+  )
+  
+  # This should handle division by zero gracefully
+  expect_warning({
+    result <- calculate_bayesian_impact(test_data)
+  }, NA) # No warning expected, but Inf values are possible
+  
+  result <- calculate_bayesian_impact(test_data)
+  expect_true(is.infinite(result$bi_rmse) || is.nan(result$bi_rmse))
+  expect_true(is.infinite(result$bi_mape) || is.nan(result$bi_mape))
+})
+
+test_that("calculate_bayesian_impact requires correct column names", {
+  # Missing required columns
+  incomplete_data <- data.frame(
+    type = c("pred", "iter_ipred"),
+    rmse = c(10, 5)
+    # Missing apriori and mape columns
+  )
+  
+  expect_error(calculate_bayesian_impact(incomplete_data))
+})
diff --git a/tests/testthat/test-calculate_shrinkage.R b/tests/testthat/test-calculate_shrinkage.R
new file mode 100644
index 0000000..3788036
--- /dev/null
+++ b/tests/testthat/test-calculate_shrinkage.R
@@ -0,0 +1,208 @@
+test_that("calculate_shrinkage returns correct structure", {
+  # Mock mod_obj with required components
+  mod_obj <- list(
+    parameters = list(cl = 1.5, v = 10.2, ka = 0.8),
+    fixed = character(0),  # No fixed parameters
+    omega = c(0.1, 0.2, 0.15)  # Triangle format: var(cl), var(v), var(ka)
+  )
+  
+  # Mock res_df with parameter estimates
+  res_df <- data.frame(
+    `_iteration` = c(1, 1, 2, 2),
+    cl = c(1.2, 1.8, 1.0, 2.0),
+    v = c(9.5, 11.0, 8.8, 12.1),
+    ka = c(0.7, 0.9, 0.6, 1.1),
+    check.names = FALSE
+  )
+  
+  result <- calculate_shrinkage(res_df, mod_obj)
+  
+  expect_s3_class(result, "data.frame")
+  expect_true("_iteration" %in% names(result))
+  expect_true(all(c("cl", "v", "ka") %in% names(result)))
+  expect_equal(nrow(result), 2)  # Two iterations
+})
+
+test_that("calculate_shrinkage handles exponential parameters correctly", {
+  mod_obj <- list(
+    parameters = list(cl = 2.0, v = 15.0),
+    fixed = character(0),
+    omega = c(0.09, 0.01, 0.16)  # var(cl) = 0.09, var(v) = 0.16
+  )
+  
+  # Create data where individual estimates vary from population
+  res_df <- data.frame(
+    `_iteration` = c(1, 1, 1),
+    cl = c(2.0, 2.2, 1.8),  # Around pop value of 2.0
+    v = c(15.0, 16.5, 13.5),  # Around pop value of 15.0
+    check.names = FALSE
+  )
+  
+  result <- calculate_shrinkage(res_df, mod_obj)
+  
+  expect_equal(nrow(result), 1)
+  expect_true(is.numeric(result$cl))
+  expect_true(is.numeric(result$v))
+  expect_true(result$cl >= 0 && result$cl <= 100)
+  expect_true(result$v >= 0 && result$v <= 100)
+})
+
+test_that("calculate_shrinkage handles kappa parameters (additive IIV)", {
+  mod_obj <- list(
+    parameters = list(kappa_cl = 0.1, kappa_v = 0.05),
+    fixed = character(0),
+    omega = c(0.01, 0.0025)  # var(kappa_cl), var(kappa_v)
+  )
+  
+  res_df <- data.frame(
+    `_iteration` = c(1, 1, 1),
+    kappa_cl = c(0.12, 0.08, 0.11),
+    kappa_v = c(0.06, 0.04, 0.055),
+    check.names = FALSE
+  )
+  
+  result <- calculate_shrinkage(res_df, mod_obj)
+  
+  expect_equal(nrow(result), 1)
+  expect_true(is.numeric(result$kappa_cl))
+  expect_true(is.numeric(result$kappa_v))
+})
+
+test_that("calculate_shrinkage handles zero population parameter", {
+  mod_obj <- list(
+    parameters = list(cl = 1.5, baseline = 0),  # Zero baseline
+    fixed = character(0),
+    omega = c(0.1, 0.01, 0.05)
+  )
+  
+  res_df <- data.frame(
+    `_iteration` = c(1, 1, 1),
+    cl = c(1.2, 1.8, 1.0),
+    baseline = c(-0.1, 0.1, 0.05),  # Should use additive for zero pop
+    check.names = FALSE
+  )
+  
+  result <- calculate_shrinkage(res_df, mod_obj)
+  
+  expect_equal(nrow(result), 1)
+  expect_true(is.numeric(result$cl))
+  expect_true(is.numeric(result$baseline))
+})
+
+test_that("calculate_shrinkage handles fixed parameters correctly", {
+  mod_obj <- list(
+    parameters = list(cl = 1.5, ka = 0.8, v = 10.2),
+    fixed = c("ka"),  # ka is fixed
+    omega = c(0.1, 0.01, 0.2)  # Only cl and v have random effects
+  )
+  
+  res_df <- data.frame(
+    `_iteration` = c(1, 1, 2, 2),
+    cl = c(1.2, 1.8, 1.0, 2.0),
+    v = c(9.5, 11.0, 8.8, 12.1),
+    ka = c(0.8, 0.8, 0.8, 0.8),  # Fixed value
+    check.names = FALSE
+  )
+  
+  result <- calculate_shrinkage(res_df, mod_obj)
+  
+  expect_true("cl" %in% names(result))
+  expect_true("v" %in% names(result))
+  expect_false("ka" %in% names(result))  # Should not appear (fixed)
+})
+
+test_that("calculate_shrinkage handles multiple iterations", {
+  mod_obj <- list(
+    parameters = list(cl = 1.5, v = 10.2),
+    fixed = character(0),
+    omega = c(0.1, 0.01, 0.2)
+  )
+  
+  res_df <- data.frame(
+    `_iteration` = c(1, 1, 1, 2, 2, 2, 3, 3, 3),
+    cl = c(1.2, 1.8, 1.0, 1.6, 1.4, 1.7, 1.1, 1.9, 1.3),
+    v = c(9.5, 11.0, 8.8, 10.8, 9.2, 11.5, 8.9, 12.0, 9.8),
+    check.names = FALSE
+  )
+  
+  result <- calculate_shrinkage(res_df, mod_obj)
+  
+  expect_equal(nrow(result), 3)  # Three iterations
+  expect_equal(result$`_iteration`, c(1, 2, 3))
+})
+
+test_that("calc_eta helper function works correctly", {
+  parameters <- list(cl = 2.0, v = 15.0, kappa_cl = 0.1, baseline = 0)
+  
+  # Test exponential parameter (cl)
+  eta_cl <- calc_eta(2.2, "cl", parameters)
+  expect_equal(eta_cl, log(2.2/2.0), tolerance = 1e-10)
+  
+  # Test kappa parameter (additive)
+  eta_kappa <- calc_eta(0.12, "kappa_cl", parameters)
+  expect_equal(eta_kappa, 0.12 - 0.1)
+  
+  # Test zero population parameter (additive)
+  eta_baseline <- calc_eta(0.05, "baseline", parameters)
+  expect_equal(eta_baseline, 0.05 - 0)
+})
+
+test_that("calc_shrinkage helper function works correctly", {
+  # Test shrinkage calculation
+  omega_list <- list(cl = 0.09)  # omega^2 = 0.09, so omega = 0.3
+  eta_values <- c(0.1, -0.1, 0.05, -0.05, 0.02)
+  
+  shrinkage <- calc_shrinkage(eta_values, "eta_cl", omega_list)
+  expected_shrinkage <- 100 * (1 - (sd(eta_values) / sqrt(0.09)))
+  
+  expect_equal(shrinkage, expected_shrinkage, tolerance = 1e-10)
+  expect_true(shrinkage >= -100 && shrinkage <= 100)  # Reasonable range
+})
+
+test_that("get_omega_for_parameters works correctly", {
+  mod_obj <- list(
+    parameters = list(cl = 1.5, v = 10.2, ka = 0.8, f1 = 1.0),
+    fixed = c("f1"),  # f1 is fixed
+    omega = c(
+      0.1,
+      0.05, 0.02, 
+      0.15, 0.01, 0.25)  # Triangle: 3x3 matrix for cl, v, ka
+  )
+  
+  omega_list <- get_omega_for_parameters(mod_obj)
+  
+  expect_equal(omega_list, list(cl = 0.1, v = 0.02, ka = 0.25))
+})
+
+test_that("calculate_shrinkage handles edge case with no variability", {
+  mod_obj <- list(
+    parameters = list(cl = 1.5),
+    fixed = character(0),
+    omega = c(0.1)
+  )
+  
+  # All individual estimates are the same (no variability)
+  res_df <- data.frame(
+    `_iteration` = c(1, 1, 1),
+    cl = c(1.5, 1.5, 1.5),  # No variation
+    check.names = FALSE
+  )
+  
+  result <- calculate_shrinkage(res_df, mod_obj)
+  
+  expect_equal(nrow(result), 1)
+  expect_equal(result$cl, 100, tolerance = 1e-10)  # Perfect shrinkage
+})
+
+test_that("calculate_shrinkage validates input structure", {
+  mod_obj <- list(
+    parameters = list(cl = 1.5),
+    fixed = character(0),
+    omega = c(0.1)
+  )
+  
+  # Missing _iteration column
+  bad_res_df <- data.frame(cl = c(1.2, 1.8, 1.0))
+  
+  expect_error(calculate_shrinkage(bad_res_df, mod_obj))
+})
diff --git a/tests/testthat/test-run_eval.R b/tests/testthat/test-run_eval.R
index 3c321b1..78d4f1f 100644
--- a/tests/testthat/test-run_eval.R
+++ b/tests/testthat/test-run_eval.R
@@ -60,7 +60,7 @@ test_that("Run also works when `model` argument just references the package", {
     progress = F,
     ids = c(1:3)
   )
-  expect_equal(names(res), c("results", "stats"))
+  expect_equal(names(res), c("results", "stats", "shrinkage", "bayesian_impact"))
   # TODO: test outputs
 })
 
@@ -79,7 +79,7 @@ test_that("Flattening of prior results in different predictions", {
     progress = F,
     ids = c(1:3)
   )
-  expect_equal(names(res_flat), c("results", "stats"))
+  expect_equal(names(res_flat), c("results", "stats", "shrinkage", "bayesian_impact"))
   expect_true(
     all(
       res$results |>

From c0d09f15120b9833d85eb37a26f2d5c7cf0414f9 Mon Sep 17 00:00:00 2001
From: Ron Keizer <ron@insight-rx.com>
Date: Tue, 23 Sep 2025 22:59:29 +0000
Subject: [PATCH 03/22] fix tests

---
 tests/testthat/test-run_eval.R         | 4 ++--
 tests/testthat/test-run_eval_grouped.R | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/tests/testthat/test-run_eval.R b/tests/testthat/test-run_eval.R
index 78d4f1f..6aae168 100644
--- a/tests/testthat/test-run_eval.R
+++ b/tests/testthat/test-run_eval.R
@@ -11,7 +11,7 @@ test_that("Basic run with vanco data + model works", {
     ruv = mod_obj$ruv,
     fixed = mod_obj$fixed,
     censor_covariates = FALSE, # shouldn't matter, since no time-varying covs
-    progress = TRUE
+    progress = FALSE
   )
 
   ## Expected structure:
@@ -104,6 +104,6 @@ test_that("Run also works when `dictionary` is used", {
     progress = F,
     ids = 1
   )
-  expect_equal(names(res), c("results", "stats"))
+  expect_equal(names(res), c("results", "stats", "shrinkage", "bayesian_impact"))
 })
 
diff --git a/tests/testthat/test-run_eval_grouped.R b/tests/testthat/test-run_eval_grouped.R
index 4b29090..3982cb7 100644
--- a/tests/testthat/test-run_eval_grouped.R
+++ b/tests/testthat/test-run_eval_grouped.R
@@ -10,7 +10,7 @@ test_that("Grouped run for busulfan/shukla works", {
     model = "pkbusulfanshukla",
     data = nm_busulfan,
     group = "group",
-    progress = TRUE,
+    progress = FALSE,
     ids = c(1:n_ids)
   )
   expect_equal(names(res), c("results", "stats", "shrinkage", "bayesian_impact"))

From 2df4445b1a9442fdd5ce1565325f912eecb280cb Mon Sep 17 00:00:00 2001
From: Ron Keizer <ron@insight-rx.com>
Date: Wed, 24 Sep 2025 18:01:08 +0000
Subject: [PATCH 04/22] minor build fixes

---
 NAMESPACE                                     |  1 +
 R/calculate_bayesian_impact.R                 | 13 +++++++---
 R/calculate_shrinkage.R                       | 16 +++++++-----
 man/calc_eta.Rd                               |  8 ++++++
 man/calculate_bayesian_impact.Rd              |  6 +++++
 .../testthat/test-calculate_bayesian_impact.R | 26 +++++++++----------
 6 files changed, 47 insertions(+), 23 deletions(-)

diff --git a/NAMESPACE b/NAMESPACE
index 44104e9..55800c2 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -4,6 +4,7 @@ S3method(parse_model,PKPDsim)
 S3method(parse_model,character)
 S3method(parse_model,default)
 export(add_grouping_column)
+export(calculate_bayesian_impact)
 export(calculate_shrinkage)
 export(compare_psn_execute_results)
 export(compare_psn_proseval_results)
diff --git a/R/calculate_bayesian_impact.R b/R/calculate_bayesian_impact.R
index 94c6b22..fdaf456 100644
--- a/R/calculate_bayesian_impact.R
+++ b/R/calculate_bayesian_impact.R
@@ -6,15 +6,20 @@
 #' (population) parameter estimates. The BI is computed based on accuracy
 #' measures (RMSE and MAPE), but not on bias estimates (MPE).
 #'
+#' @param stats_summ stats summary, a data.frame created in `run_eval()`
+#'
+#' @returns a tibble with bayesian impact values, based on RMSE and MAPE
+#'
+#' @export
 calculate_bayesian_impact <- function(
   stats_summ
 ) {
   stats_summ |>
-    dplyr::filter(!apriori) |>
-    dplyr::select(type, apriori, rmse, mape) |>
+    dplyr::filter(! .data$apriori) |>
+    dplyr::select("type", "apriori", "rmse", "mape") |>
     tidyr::pivot_wider(names_from = "type", values_from = c("rmse", "mape")) |>
     dplyr::summarise(
-      bi_rmse = 100 * (1 - (rmse_iter_ipred / rmse_pred)),
-      bi_mape = 100 * (1 - (mape_iter_ipred / mape_pred))
+      bi_rmse = 100 * (1 - (.data$rmse_iter_ipred / .data$rmse_pred)),
+      bi_mape = 100 * (1 - (.data$mape_iter_ipred / .data$mape_pred))
     )
 }
diff --git a/R/calculate_shrinkage.R b/R/calculate_shrinkage.R
index 382f231..8017283 100644
--- a/R/calculate_shrinkage.R
+++ b/R/calculate_shrinkage.R
@@ -25,13 +25,13 @@ calculate_shrinkage <- function(
 ) {
   om <- get_omega_for_parameters(mod_obj)
   shr_df <- res_df |>
-    dplyr::mutate(across(
+    dplyr::mutate(dplyr::across(
       .cols = names(om),
       .fns  = ~ calc_eta(.x, dplyr::cur_column(), mod_obj$parameters),
       .names = "eta_{.col}"
     )) |>
-    dplyr::group_by(`_iteration`) |>
-    dplyr::summarise(across(
+    dplyr::group_by(.data$`_iteration`) |>
+    dplyr::summarise(dplyr::across(
       .cols = paste0("eta_", names(om)),
       .fns = ~ calc_shrinkage(.x, dplyr::cur_column(), om),
       .names = "shr_{.col}"
@@ -45,9 +45,13 @@ calculate_shrinkage <- function(
 #' additive shape for IOV (kappa's), i.e. `PAR = TV_PAR + ETA(n)` or when
 #' final parameter estimate is 0.
 #'
+#' @param ind individual parameter estimate
+#' @param par_name name of parameter in model
+#' @param parameters list of parameter values, should include entry for
+#'   `par_name`
 calc_eta <- function(ind, par_name, parameters) {
   pop <- parameters[[par_name]]
-  if(stringr::str_detect(par_name, "^kappa_") || pop == 0) { ## PAR = TV_PAR + ETA()
+  if(grepl("^kappa_", par_name) || pop == 0) { ## PAR = TV_PAR + ETA()
     return(ind - pop)
   } else {
     return(log(ind/pop)) ## PAR = TV_PAR + ETA()
@@ -55,13 +59,13 @@ calc_eta <- function(ind, par_name, parameters) {
 }
 
 calc_shrinkage <- function(eta, par_name, om) {
-  100 * (1 - (sd(eta) / sqrt(om[[gsub("^eta_", "", par_name)]])))
+  100 * (1 - (stats::sd(eta) / sqrt(om[[gsub("^eta_", "", par_name)]])))
 }
 
 #' Return all diagonal om^2 elements for each non-fixed parameter, as a list
 #'
 #' @inheritParams run_eval_core
-#' 
+#'
 #' @returns list
 get_omega_for_parameters <- function(mod_obj) {
   pars <- names(mod_obj$parameters)
diff --git a/man/calc_eta.Rd b/man/calc_eta.Rd
index 54caa07..bcd1210 100644
--- a/man/calc_eta.Rd
+++ b/man/calc_eta.Rd
@@ -9,6 +9,14 @@ final parameter estimate is 0.}
 \usage{
 calc_eta(ind, par_name, parameters)
 }
+\arguments{
+\item{ind}{individual parameter estimate}
+
+\item{par_name}{name of parameter in model}
+
+\item{parameters}{list of parameter values, should include entry for
+\code{par_name}}
+}
 \description{
 Calculate the "eta"-value for a parameters, assuming an exponential
 shape for IIV, i.e. \code{PAR = TV_PAR * EXP(ETA(n))}, and
diff --git a/man/calculate_bayesian_impact.Rd b/man/calculate_bayesian_impact.Rd
index cec8568..a552fe8 100644
--- a/man/calculate_bayesian_impact.Rd
+++ b/man/calculate_bayesian_impact.Rd
@@ -7,6 +7,12 @@ to using population estimates}
 \usage{
 calculate_bayesian_impact(stats_summ)
 }
+\arguments{
+\item{stats_summ}{stats summary, a data.frame created in \code{run_eval()}}
+}
+\value{
+a tibble with bayesian impact values, based on RMSE and MAPE
+}
 \description{
 Bayesian impact (BI), as defined here, is the gain in accuracy that is
 obtained when using Bayesian updating, compared to using a priori
diff --git a/tests/testthat/test-calculate_bayesian_impact.R b/tests/testthat/test-calculate_bayesian_impact.R
index c1ef1d8..cca45b0 100644
--- a/tests/testthat/test-calculate_bayesian_impact.R
+++ b/tests/testthat/test-calculate_bayesian_impact.R
@@ -6,9 +6,9 @@ test_that("calculate_bayesian_impact returns correct structure", {
     rmse = c(10, 5, 15, 8),
     mape = c(20, 10, 25, 12)
   )
-  
+
   result <- calculate_bayesian_impact(test_data)
-  
+
   expect_s3_class(result, "data.frame")
   expect_named(result, c("bi_rmse", "bi_mape"))
   expect_equal(nrow(result), 1)
@@ -22,9 +22,9 @@ test_that("calculate_bayesian_impact calculates correct values", {
     rmse = c(10, 5),  # iter_ipred is 50% better
     mape = c(20, 10)  # iter_ipred is 50% better
   )
-  
+
   result <- calculate_bayesian_impact(test_data)
-  
+
   # BI = 100 * (1 - (iter_ipred / pred))
   # For RMSE: 100 * (1 - (5/10)) = 50%
   # For MAPE: 100 * (1 - (10/20)) = 50%
@@ -40,9 +40,9 @@ test_that("calculate_bayesian_impact handles negative impact correctly", {
     rmse = c(5, 10),   # iter_ipred is worse
     mape = c(10, 20)   # iter_ipred is worse
   )
-  
+
   result <- calculate_bayesian_impact(test_data)
-  
+
   # BI = 100 * (1 - (10/5)) = -100%
   expect_equal(result$bi_rmse, -100)
   expect_equal(result$bi_mape, -100)
@@ -56,9 +56,9 @@ test_that("calculate_bayesian_impact filters apriori correctly", {
     rmse = c(100, 50, 10, 5),
     mape = c(200, 100, 20, 10)
   )
-  
+
   result <- calculate_bayesian_impact(test_data)
-  
+
   # Should only use the apriori = FALSE rows (10, 5 and 20, 10)
   expect_equal(result$bi_rmse, 50)
   expect_equal(result$bi_mape, 50)
@@ -72,9 +72,9 @@ test_that("calculate_bayesian_impact handles edge cases", {
     rmse = c(10, 10),
     mape = c(20, 20)
   )
-  
+
   result <- calculate_bayesian_impact(test_data)
-  
+
   expect_equal(result$bi_rmse, 0)
   expect_equal(result$bi_mape, 0)
 })
@@ -87,12 +87,12 @@ test_that("calculate_bayesian_impact handles division by zero", {
     rmse = c(0, 5),
     mape = c(0, 10)
   )
-  
+
   # This should handle division by zero gracefully
   expect_warning({
     result <- calculate_bayesian_impact(test_data)
   }, NA) # No warning expected, but Inf values are possible
-  
+
   result <- calculate_bayesian_impact(test_data)
   expect_true(is.infinite(result$bi_rmse) || is.nan(result$bi_rmse))
   expect_true(is.infinite(result$bi_mape) || is.nan(result$bi_mape))
@@ -105,6 +105,6 @@ test_that("calculate_bayesian_impact requires correct column names", {
     rmse = c(10, 5)
     # Missing apriori and mape columns
   )
-  
+
   expect_error(calculate_bayesian_impact(incomplete_data))
 })

From 683b08ea7571359745f77110df7323a664f1ebc7 Mon Sep 17 00:00:00 2001
From: roninsightrx <ron@insight-rx.com>
Date: Thu, 25 Sep 2025 12:23:54 -0700
Subject: [PATCH 05/22] Apply suggestion from @mccarthy-m-g

Co-authored-by: Michael McCarthy <51542091+mccarthy-m-g@users.noreply.github.com>
---
 R/calculate_bayesian_impact.R | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/R/calculate_bayesian_impact.R b/R/calculate_bayesian_impact.R
index fdaf456..b663bb3 100644
--- a/R/calculate_bayesian_impact.R
+++ b/R/calculate_bayesian_impact.R
@@ -1,5 +1,4 @@
-#' This function calculates the impact of using Bayesian updating, compared
-#' to using population estimates
+#' Calculate the impact of using Bayesian updating compared to population estimates
 #'
 #' Bayesian impact (BI), as defined here, is the gain in accuracy that is
 #' obtained when using Bayesian updating, compared to using a priori

From c805aa94d534e0192f50f39d4949e50e4913e3ca Mon Sep 17 00:00:00 2001
From: roninsightrx <ron@insight-rx.com>
Date: Thu, 25 Sep 2025 12:24:03 -0700
Subject: [PATCH 06/22] Apply suggestion from @mccarthy-m-g

Co-authored-by: Michael McCarthy <51542091+mccarthy-m-g@users.noreply.github.com>
---
 R/calculate_bayesian_impact.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/calculate_bayesian_impact.R b/R/calculate_bayesian_impact.R
index b663bb3..5eff180 100644
--- a/R/calculate_bayesian_impact.R
+++ b/R/calculate_bayesian_impact.R
@@ -14,7 +14,7 @@ calculate_bayesian_impact <- function(
   stats_summ
 ) {
   stats_summ |>
-    dplyr::filter(! .data$apriori) |>
+    dplyr::filter(!.data$apriori) |>
     dplyr::select("type", "apriori", "rmse", "mape") |>
     tidyr::pivot_wider(names_from = "type", values_from = c("rmse", "mape")) |>
     dplyr::summarise(

From 984349735e19dc14a2f4d68ecf5456f73a7fa5f8 Mon Sep 17 00:00:00 2001
From: roninsightrx <ron@insight-rx.com>
Date: Thu, 25 Sep 2025 12:25:27 -0700
Subject: [PATCH 07/22] Apply suggestion from @mccarthy-m-g

Co-authored-by: Michael McCarthy <51542091+mccarthy-m-g@users.noreply.github.com>
---
 R/calculate_shrinkage.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/calculate_shrinkage.R b/R/calculate_shrinkage.R
index 8017283..8587dcc 100644
--- a/R/calculate_shrinkage.R
+++ b/R/calculate_shrinkage.R
@@ -1,4 +1,4 @@
-#' Function to calculate eta-shrinkage
+#' Calculate eta-shrinkage
 #'
 #' Shrinkage for population PK models was first defined in this paper by
 #' Savic and Karlsson: https://pmc.ncbi.nlm.nih.gov/articles/PMC2758126/.

From 73edcbe122d841bec6ffcddf9883505c3f4dcd32 Mon Sep 17 00:00:00 2001
From: roninsightrx <ron@insight-rx.com>
Date: Thu, 25 Sep 2025 12:25:57 -0700
Subject: [PATCH 08/22] Apply suggestion from @mccarthy-m-g

Co-authored-by: Michael McCarthy <51542091+mccarthy-m-g@users.noreply.github.com>
---
 R/calculate_shrinkage.R | 14 ++------------
 1 file changed, 2 insertions(+), 12 deletions(-)

diff --git a/R/calculate_shrinkage.R b/R/calculate_shrinkage.R
index 8587dcc..37948ae 100644
--- a/R/calculate_shrinkage.R
+++ b/R/calculate_shrinkage.R
@@ -1,17 +1,7 @@
 #' Calculate eta-shrinkage
 #'
-#' Shrinkage for population PK models was first defined in this paper by
-#' Savic and Karlsson: https://pmc.ncbi.nlm.nih.gov/articles/PMC2758126/.
-#' It is a measure of how much information is available to update individual
-#' estimates away from the population value. In principle, if there is no
-#' information at all (i.e. in the case of population estimates only),
-#' shrinkage will be 100%. In the case of fully informed Bayesian estimates
-#' (unlikely to be achieved in practice), shrinkage is 0%. In most practical
-#' scenarios with limited sampling, shrinkage will be between 10-40%. When
-#' shrinkage is higher than 50% or so, one could argue there is limited benefit
-#' of sampling at all. So the shrinkage reported in this package can be used
-#' to evaluate whether the sampling was efficient, and how shrinkage could be
-#' reduced with additional sampling.
+#' Calculate eta-shrinkage, measure of how much information is available to
+#' update individual estimates away from the population value.
 #'
 #' @inheritParams run_eval_core
 #' @param res_df data.frame or tibble with raw results (including parameters)

From 1590ea9efdaa06d8a513c5800e50b51d430f021a Mon Sep 17 00:00:00 2001
From: roninsightrx <ron@insight-rx.com>
Date: Thu, 25 Sep 2025 12:26:13 -0700
Subject: [PATCH 09/22] Apply suggestion from @mccarthy-m-g

Co-authored-by: Michael McCarthy <51542091+mccarthy-m-g@users.noreply.github.com>
---
 R/calculate_shrinkage.R | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/R/calculate_shrinkage.R b/R/calculate_shrinkage.R
index 37948ae..ec58895 100644
--- a/R/calculate_shrinkage.R
+++ b/R/calculate_shrinkage.R
@@ -7,6 +7,20 @@
 #' @param res_df data.frame or tibble with raw results (including parameters)
 #' from a run with `mipdeval::run_eval()`.
 #'
+#' @details
+#' Shrinkage for population PK models was first defined in this paper by
+#' Savic and Karlsson: https://pmc.ncbi.nlm.nih.gov/articles/PMC2758126/.
+#' It is a measure of how much information is available to update individual
+#' estimates away from the population value. In principle, if there is no
+#' information at all (i.e. in the case of population estimates only),
+#' shrinkage will be 100%. In the case of fully informed Bayesian estimates
+#' (unlikely to be achieved in practice), shrinkage is 0%. In most practical
+#' scenarios with limited sampling, shrinkage will be between 10-40%. When
+#' shrinkage is higher than 50% or so, one could argue there is limited benefit
+#' of sampling at all. So the shrinkage reported in this package can be used
+#' to evaluate whether the sampling was efficient, and how shrinkage could be
+#' reduced with additional sampling.
+#'
 #' @returns tibble
 #' @export
 calculate_shrinkage <- function(

From ab21eb10d76d2fd8285eb6cb9cbf7aa892a5ca5e Mon Sep 17 00:00:00 2001
From: roninsightrx <ron@insight-rx.com>
Date: Thu, 25 Sep 2025 12:26:26 -0700
Subject: [PATCH 10/22] Apply suggestion from @mccarthy-m-g

Co-authored-by: Michael McCarthy <51542091+mccarthy-m-g@users.noreply.github.com>
---
 R/calculate_shrinkage.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/calculate_shrinkage.R b/R/calculate_shrinkage.R
index ec58895..789577a 100644
--- a/R/calculate_shrinkage.R
+++ b/R/calculate_shrinkage.R
@@ -73,7 +73,7 @@ calc_shrinkage <- function(eta, par_name, om) {
 #' @returns list
 get_omega_for_parameters <- function(mod_obj) {
   pars <- names(mod_obj$parameters)
-  pars <- pars[! (pars %in% mod_obj$fixed)]
+  pars <- pars[!(pars %in% mod_obj$fixed)]
   om <- as.list(
     diag(
       PKPDsim::triangle_to_full(mod_obj$omega)

From 47de80603d76706c8347dbfcd531d87aa429801f Mon Sep 17 00:00:00 2001
From: roninsightrx <ron@insight-rx.com>
Date: Thu, 25 Sep 2025 12:26:33 -0700
Subject: [PATCH 11/22] Apply suggestion from @mccarthy-m-g

Co-authored-by: Michael McCarthy <51542091+mccarthy-m-g@users.noreply.github.com>
---
 R/calculate_shrinkage.R | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/R/calculate_shrinkage.R b/R/calculate_shrinkage.R
index 789577a..b306ae2 100644
--- a/R/calculate_shrinkage.R
+++ b/R/calculate_shrinkage.R
@@ -74,11 +74,7 @@ calc_shrinkage <- function(eta, par_name, om) {
 get_omega_for_parameters <- function(mod_obj) {
   pars <- names(mod_obj$parameters)
   pars <- pars[!(pars %in% mod_obj$fixed)]
-  om <- as.list(
-    diag(
-      PKPDsim::triangle_to_full(mod_obj$omega)
-    )[1:length(pars)]
-  )
+    om <- as.list(diag(PKPDsim::triangle_to_full(mod_obj$omega))[1:length(pars)])
   names(om) <- pars
   om
 }

From 9d692be7f4c3c55a1db1bf93f4a8344765fbbfcf Mon Sep 17 00:00:00 2001
From: roninsightrx <ron@insight-rx.com>
Date: Thu, 25 Sep 2025 12:26:42 -0700
Subject: [PATCH 12/22] Apply suggestion from @mccarthy-m-g

Co-authored-by: Michael McCarthy <51542091+mccarthy-m-g@users.noreply.github.com>
---
 R/run_eval.R | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/R/run_eval.R b/R/run_eval.R
index 5414c5c..a6b93ec 100644
--- a/R/run_eval.R
+++ b/R/run_eval.R
@@ -141,10 +141,7 @@ run_eval <- function(
       mpe = mpe(.data$dv, .data$value),
       mape = mape(.data$dv, .data$value)
     )
-  shrinkage_summ <- calculate_shrinkage(
-    res_df,
-    mod_obj = mod_obj
-  )
+  shrinkage_summ <- calculate_shrinkage(res_df, mod_obj = mod_obj)
   bayesian_impact <- calculate_bayesian_impact(
     stats_summ
   )

From 59dad4102df42f7669a73444d4660d7809a77305 Mon Sep 17 00:00:00 2001
From: roninsightrx <ron@insight-rx.com>
Date: Thu, 25 Sep 2025 12:27:18 -0700
Subject: [PATCH 13/22] Apply suggestion from @mccarthy-m-g

Co-authored-by: Michael McCarthy <51542091+mccarthy-m-g@users.noreply.github.com>
---
 tests/testthat/test-calculate_bayesian_impact.R | 9 +++------
 1 file changed, 3 insertions(+), 6 deletions(-)

diff --git a/tests/testthat/test-calculate_bayesian_impact.R b/tests/testthat/test-calculate_bayesian_impact.R
index cca45b0..3d2b2f1 100644
--- a/tests/testthat/test-calculate_bayesian_impact.R
+++ b/tests/testthat/test-calculate_bayesian_impact.R
@@ -88,12 +88,9 @@ test_that("calculate_bayesian_impact handles division by zero", {
     mape = c(0, 10)
   )
 
-  # This should handle division by zero gracefully
-  expect_warning({
-    result <- calculate_bayesian_impact(test_data)
-  }, NA) # No warning expected, but Inf values are possible
-
-  result <- calculate_bayesian_impact(test_data)
+  # This should handle division by zero gracefully. No warning expected, but Inf
+  # values are possible.
+  result <- expect_no_warning(calculate_bayesian_impact(test_data))
   expect_true(is.infinite(result$bi_rmse) || is.nan(result$bi_rmse))
   expect_true(is.infinite(result$bi_mape) || is.nan(result$bi_mape))
 })

From 4b22613be0494db430ed3fe54d75070e2d0119e2 Mon Sep 17 00:00:00 2001
From: roninsightrx <ron@insight-rx.com>
Date: Thu, 25 Sep 2025 12:27:41 -0700
Subject: [PATCH 14/22] Apply suggestion from @mccarthy-m-g

Co-authored-by: Michael McCarthy <51542091+mccarthy-m-g@users.noreply.github.com>
---
 R/run_eval.R | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/R/run_eval.R b/R/run_eval.R
index a6b93ec..ef86c14 100644
--- a/R/run_eval.R
+++ b/R/run_eval.R
@@ -142,9 +142,7 @@ run_eval <- function(
       mape = mape(.data$dv, .data$value)
     )
   shrinkage_summ <- calculate_shrinkage(res_df, mod_obj = mod_obj)
-  bayesian_impact <- calculate_bayesian_impact(
-    stats_summ
-  )
+  bayesian_impact <- calculate_bayesian_impact(stats_summ)
 
   out <- list(
     results = tibble::as_tibble(res_df),

From 4c7087564ee54d858885a2abb36f0cb31645e9f6 Mon Sep 17 00:00:00 2001
From: Ron Keizer <ron@insight-rx.com>
Date: Thu, 25 Sep 2025 21:03:31 +0000
Subject: [PATCH 15/22] move to S3 + print function + few more stats

---
 NAMESPACE                                     |  5 ++
 R/calculate_bayesian_impact.R                 | 12 ++--
 R/calculate_shrinkage.R                       | 27 ++++-----
 R/calculate_stats.R                           | 32 +++++++++++
 R/misc.R                                      | 17 +++++-
 R/print_functions.R                           | 56 +++++++++++++++++++
 R/run_eval.R                                  | 29 +++-------
 R/run_eval_core.R                             |  3 -
 man/calculate_bayesian_impact.Rd              |  8 +--
 man/calculate_shrinkage.Rd                    | 14 +++--
 man/calculate_stats.Rd                        | 20 +++++++
 man/mape.Rd                                   |  4 +-
 man/mpe.Rd                                    |  4 +-
 man/nm_busulfan.Rd                            |  2 +-
 man/nrmse.Rd                                  | 19 +++++++
 man/print.mipdeval_results.Rd                 | 16 ++++++
 man/print.mipdeval_results_bayesian_impact.Rd | 16 ++++++
 man/print.mipdeval_results_shrinkage.Rd       | 16 ++++++
 man/print.mipdeval_results_stats_summ.Rd      | 16 ++++++
 man/run_eval_core.Rd                          | 10 ----
 20 files changed, 258 insertions(+), 68 deletions(-)
 create mode 100644 R/calculate_stats.R
 create mode 100644 R/print_functions.R
 create mode 100644 man/calculate_stats.Rd
 create mode 100644 man/nrmse.Rd
 create mode 100644 man/print.mipdeval_results.Rd
 create mode 100644 man/print.mipdeval_results_bayesian_impact.Rd
 create mode 100644 man/print.mipdeval_results_shrinkage.Rd
 create mode 100644 man/print.mipdeval_results_stats_summ.Rd

diff --git a/NAMESPACE b/NAMESPACE
index 55800c2..79ebd02 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -3,9 +3,14 @@
 S3method(parse_model,PKPDsim)
 S3method(parse_model,character)
 S3method(parse_model,default)
+S3method(print,mipdeval_results)
+S3method(print,mipdeval_results_bayesian_impact)
+S3method(print,mipdeval_results_shrinkage)
+S3method(print,mipdeval_results_stats_summ)
 export(add_grouping_column)
 export(calculate_bayesian_impact)
 export(calculate_shrinkage)
+export(calculate_stats)
 export(compare_psn_execute_results)
 export(compare_psn_proseval_results)
 export(group_by_dose)
diff --git a/R/calculate_bayesian_impact.R b/R/calculate_bayesian_impact.R
index 5eff180..29332b0 100644
--- a/R/calculate_bayesian_impact.R
+++ b/R/calculate_bayesian_impact.R
@@ -5,20 +5,22 @@
 #' (population) parameter estimates. The BI is computed based on accuracy
 #' measures (RMSE and MAPE), but not on bias estimates (MPE).
 #'
-#' @param stats_summ stats summary, a data.frame created in `run_eval()`
+#' @inheritParams calculate_stats
 #'
 #' @returns a tibble with bayesian impact values, based on RMSE and MAPE
 #'
 #' @export
 calculate_bayesian_impact <- function(
-  stats_summ
+  res
 ) {
-  stats_summ |>
+  out <- res$stats_summ |>
     dplyr::filter(!.data$apriori) |>
     dplyr::select("type", "apriori", "rmse", "mape") |>
     tidyr::pivot_wider(names_from = "type", values_from = c("rmse", "mape")) |>
     dplyr::summarise(
-      bi_rmse = 100 * (1 - (.data$rmse_iter_ipred / .data$rmse_pred)),
-      bi_mape = 100 * (1 - (.data$mape_iter_ipred / .data$mape_pred))
+      bi_rmse = round(100 * (1 - (.data$rmse_iter_ipred / .data$rmse_pred)), 1),
+      bi_mape = round(100 * (1 - (.data$mape_iter_ipred / .data$mape_pred)), 1)
     )
+  class(out) <- c("mipdeval_results_bayesian_impact", class(out))
+  out
 }
diff --git a/R/calculate_shrinkage.R b/R/calculate_shrinkage.R
index b306ae2..6604634 100644
--- a/R/calculate_shrinkage.R
+++ b/R/calculate_shrinkage.R
@@ -3,9 +3,7 @@
 #' Calculate eta-shrinkage, measure of how much information is available to
 #' update individual estimates away from the population value.
 #'
-#' @inheritParams run_eval_core
-#' @param res_df data.frame or tibble with raw results (including parameters)
-#' from a run with `mipdeval::run_eval()`.
+#' @inheritParams calculate_stats
 #'
 #' @details
 #' Shrinkage for population PK models was first defined in this paper by
@@ -22,16 +20,14 @@
 #' reduced with additional sampling.
 #'
 #' @returns tibble
+#'
 #' @export
-calculate_shrinkage <- function(
-  res_df,
-  mod_obj
-) {
-  om <- get_omega_for_parameters(mod_obj)
-  shr_df <- res_df |>
+calculate_shrinkage <- function(res) {
+  om <- get_omega_for_parameters(res$mod_obj)
+  out <- res$results |>
     dplyr::mutate(dplyr::across(
       .cols = names(om),
-      .fns  = ~ calc_eta(.x, dplyr::cur_column(), mod_obj$parameters),
+      .fns  = ~ calc_eta(.x, dplyr::cur_column(), res$mod_obj$parameters),
       .names = "eta_{.col}"
     )) |>
     dplyr::group_by(.data$`_iteration`) |>
@@ -40,8 +36,10 @@ calculate_shrinkage <- function(
       .fns = ~ calc_shrinkage(.x, dplyr::cur_column(), om),
       .names = "shr_{.col}"
     ))
-  names(shr_df) <- gsub("^shr_eta_", "", names(shr_df))
-  shr_df
+  out <- out |>
+    rlang::set_names(gsub("^shr_eta_", "", names(out)))
+  class(out) <- c("mipdeval_results_shrinkage", class(out))
+  out
 }
 
 #' Calculate the "eta"-value for a parameters, assuming an exponential
@@ -63,7 +61,10 @@ calc_eta <- function(ind, par_name, parameters) {
 }
 
 calc_shrinkage <- function(eta, par_name, om) {
-  100 * (1 - (stats::sd(eta) / sqrt(om[[gsub("^eta_", "", par_name)]])))
+  round(
+    100 * (1 - (stats::sd(eta) / sqrt(om[[gsub("^eta_", "", par_name)]]))),
+    1
+  )
 }
 
 #' Return all diagonal om^2 elements for each non-fixed parameter, as a list
diff --git a/R/calculate_stats.R b/R/calculate_stats.R
new file mode 100644
index 0000000..7de3a71
--- /dev/null
+++ b/R/calculate_stats.R
@@ -0,0 +1,32 @@
+#' Calculate basic statistics, like RMSE, MPE, MAPE for forecasted data
+#'
+#' @param res output object (`mipdeval_results`) from `run_eval()`, or
+#' `data.frame` with raw results.
+#' @param rounding number of decimals to round to.
+#'
+#' @returns tibble
+#'
+#' @export
+calculate_stats <- function(
+    res,
+    rounding = 3
+) {
+  if(inherits(res, "mipdeval_results")) {
+    res <- res$results
+  }
+  out <- res |>
+    tidyr::pivot_longer(
+      cols = c("pred", "map_ipred", "iter_ipred"), names_to = "type"
+    ) |>
+    dplyr::group_by(.data$type, .data$apriori) |>
+    dplyr::summarise(
+      rmse = rmse(.data$dv, .data$value),
+      nrmse = nrmse(.data$dv, .data$value),
+      mpe = mpe(.data$dv, .data$value),
+      mape = mape(.data$dv, .data$value)
+    ) |>
+    dplyr::mutate(dplyr::across(dplyr::everything(), round, rounding)) |>
+    dplyr::as_tibble()
+  class(out) <- c("mipdeval_results_stats_summ", class(out))
+  out
+}
diff --git a/R/misc.R b/R/misc.R
index 27aff52..ac76d0e 100644
--- a/R/misc.R
+++ b/R/misc.R
@@ -37,7 +37,20 @@ rmse <- function (obs, pred) {
   sqrt(mean(res_sq, na.rm = TRUE))
 }
 
-#' Mean absolute prediction error
+#' Normalized root-mean-squared error
+#'
+#' @param obs observations vector
+#' @param pred predictions vector
+#'
+#' @returns A numeric vector
+#'
+nrmse <- function (obs, pred) {
+  res_sq <- (pred - obs)^2
+  rmse <- sqrt(mean(res_sq, na.rm = T))
+  rmse/mean(obs, na.rm = T)
+}
+
+#' Mean absolute percentage error
 #'
 #' @inheritParams rmse
 #'
@@ -46,7 +59,7 @@ mape <- function (obs, pred) {
   sum(abs((obs - pred))/obs)/length(obs)
 }
 
-#' Mean prediction error
+#' Mean percentage error
 #'
 #' @inheritParams rmse
 #'
diff --git a/R/print_functions.R b/R/print_functions.R
new file mode 100644
index 0000000..e15dd2b
--- /dev/null
+++ b/R/print_functions.R
@@ -0,0 +1,56 @@
+#' Print results from run_eval()
+#'
+#' @param x output object from `run_eval()`.
+#' @param ... further arguments passed to or from other methods.
+#'
+#' @export
+print.mipdeval_results <- function(x, ...) {
+  print(x$stats_summ, ...)
+  print(x$shrinkage, ...)
+  print(x$bayesian_impact, ...)
+}
+
+#' Print predictive performance statistics from run_eval()
+#'
+#' @inheritParams print.mipdeval_results
+#' @param x stats summary object in output from `run_eval()`.
+#'
+#' @export
+print.mipdeval_results_stats_summ <- function(x, ...) {
+  x |>
+    dplyr::mutate(
+      type = dplyr::case_when(
+        type == "iter_ipred" ~ "Iterative forecasting",
+        type == "map_ipred" ~ "Retrospective MAP Bayesian",
+        type == "pred" ~ "Population prediction"
+      )
+    ) |>
+    knitr::kable(format = "simple", caption = "Predictive performance") |>
+    print(...)
+}
+
+#' Print shrinkage results from run_eval()
+#'
+#' @inheritParams print.mipdeval_results
+#' @param x shrinkage object in output from `run_eval()`.
+#'
+#' @export
+print.mipdeval_results_shrinkage <- function(x, ...) {
+  x |>
+    dplyr::rename("Iteration" = .data$`_iteration`) |>
+    knitr::kable(format = "simple", caption = "Shrinkage (%)") |>
+    print(...)
+}
+
+#' Print Bayesian impact results from run_eval()
+#'
+#' @inheritParams print.mipdeval_results
+#' @param x bayesian impact object in output from `run_eval()`.
+#'
+#' @export
+print.mipdeval_results_bayesian_impact <- function(x, ...) {
+  x |>
+    rlang::set_names(gsub("bi_", "", names(x))) |>
+    knitr::kable(format = "simple", caption = "Bayesian impact (%), based on RMSE / MAPE") |>
+    print(...)
+}
diff --git a/R/run_eval.R b/R/run_eval.R
index ef86c14..be6cc20 100644
--- a/R/run_eval.R
+++ b/R/run_eval.R
@@ -130,28 +130,17 @@ run_eval <- function(
 
   ## 4. Combine results and basic stats into return object
   if(verbose) cli::cli_alert_info("Calculating forecasting statistics")
-  res_df <- dplyr::bind_rows(res)
-  stats_summ <- res_df |>
-    tidyr::pivot_longer(
-      cols = c("pred", "map_ipred", "iter_ipred"), names_to = "type"
-    ) |>
-    dplyr::group_by(.data$type, .data$apriori) |>
-    dplyr::summarise(
-      rmse = rmse(.data$dv, .data$value),
-      mpe = mpe(.data$dv, .data$value),
-      mape = mape(.data$dv, .data$value)
-    )
-  shrinkage_summ <- calculate_shrinkage(res_df, mod_obj = mod_obj)
-  bayesian_impact <- calculate_bayesian_impact(stats_summ)
-
+  res_tbl <- dplyr::bind_rows(res) |>
+    tibble::as_tibble()
   out <- list(
-    results = tibble::as_tibble(res_df),
-    stats = tibble::as_tibble(stats_summ),
-    shrinkage = shrinkage_summ,
-    bayesian_impact = bayesian_impact
+    results = res_tbl,
+    mod_obj = mod_obj,
+    data = data
   )
-
-  # TODO: Turn out into an S3 class, so we can give it methods like print(), etc.
+  class(out) <- c("mipdeval_results", "list")
+  out$stats_summ <- calculate_stats(out)
+  out$shrinkage <- calculate_shrinkage(out)
+  out$bayesian_impact <- calculate_bayesian_impact(out)
 
   ## 5. Return results
   if(verbose) cli::cli_alert_info("Done")
diff --git a/R/run_eval_core.R b/R/run_eval_core.R
index 0e2bd43..b1ef6ac 100644
--- a/R/run_eval_core.R
+++ b/R/run_eval_core.R
@@ -10,9 +10,6 @@
 run_eval_core <- function(
   mod_obj,
   data,
-  parameters = NULL,
-  omega = NULL,
-  ruv = NULL,
   weights = NULL,
   weight_prior = 1,
   censor_covariates = TRUE,
diff --git a/man/calculate_bayesian_impact.Rd b/man/calculate_bayesian_impact.Rd
index a552fe8..b4c0edf 100644
--- a/man/calculate_bayesian_impact.Rd
+++ b/man/calculate_bayesian_impact.Rd
@@ -2,13 +2,13 @@
 % Please edit documentation in R/calculate_bayesian_impact.R
 \name{calculate_bayesian_impact}
 \alias{calculate_bayesian_impact}
-\title{This function calculates the impact of using Bayesian updating, compared
-to using population estimates}
+\title{Calculate the impact of using Bayesian updating compared to population estimates}
 \usage{
-calculate_bayesian_impact(stats_summ)
+calculate_bayesian_impact(res)
 }
 \arguments{
-\item{stats_summ}{stats summary, a data.frame created in \code{run_eval()}}
+\item{res}{output object (\code{mipdeval_results}) from \code{run_eval()}, or
+\code{data.frame} with raw results.}
 }
 \value{
 a tibble with bayesian impact values, based on RMSE and MAPE
diff --git a/man/calculate_shrinkage.Rd b/man/calculate_shrinkage.Rd
index 03d289c..98c48f6 100644
--- a/man/calculate_shrinkage.Rd
+++ b/man/calculate_shrinkage.Rd
@@ -2,20 +2,22 @@
 % Please edit documentation in R/calculate_shrinkage.R
 \name{calculate_shrinkage}
 \alias{calculate_shrinkage}
-\title{Function to calculate eta-shrinkage}
+\title{Calculate eta-shrinkage}
 \usage{
-calculate_shrinkage(res_df, mod_obj)
+calculate_shrinkage(res)
 }
 \arguments{
-\item{res_df}{data.frame or tibble with raw results (including parameters)
-from a run with \code{mipdeval::run_eval()}.}
-
-\item{mod_obj}{list object with model information}
+\item{res}{output object (\code{mipdeval_results}) from \code{run_eval()}, or
+\code{data.frame} with raw results.}
 }
 \value{
 tibble
 }
 \description{
+Calculate eta-shrinkage, measure of how much information is available to
+update individual estimates away from the population value.
+}
+\details{
 Shrinkage for population PK models was first defined in this paper by
 Savic and Karlsson: https://pmc.ncbi.nlm.nih.gov/articles/PMC2758126/.
 It is a measure of how much information is available to update individual
diff --git a/man/calculate_stats.Rd b/man/calculate_stats.Rd
new file mode 100644
index 0000000..29f44b7
--- /dev/null
+++ b/man/calculate_stats.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/calculate_stats.R
+\name{calculate_stats}
+\alias{calculate_stats}
+\title{Calculate basic statistics, like RMSE, MPE, MAPE for forecasted data}
+\usage{
+calculate_stats(res, rounding = 3)
+}
+\arguments{
+\item{res}{output object (\code{mipdeval_results}) from \code{run_eval()}, or
+\code{data.frame} with raw results.}
+
+\item{rounding}{number of decimals to round to.}
+}
+\value{
+tibble
+}
+\description{
+Calculate basic statistics, like RMSE, MPE, MAPE for forecasted data
+}
diff --git a/man/mape.Rd b/man/mape.Rd
index 152df39..c002469 100644
--- a/man/mape.Rd
+++ b/man/mape.Rd
@@ -2,7 +2,7 @@
 % Please edit documentation in R/misc.R
 \name{mape}
 \alias{mape}
-\title{Mean absolute prediction error}
+\title{Mean absolute percentage error}
 \usage{
 mape(obs, pred)
 }
@@ -15,5 +15,5 @@ mape(obs, pred)
 A numeric vector
 }
 \description{
-Mean absolute prediction error
+Mean absolute percentage error
 }
diff --git a/man/mpe.Rd b/man/mpe.Rd
index b55efac..8cf8872 100644
--- a/man/mpe.Rd
+++ b/man/mpe.Rd
@@ -2,7 +2,7 @@
 % Please edit documentation in R/misc.R
 \name{mpe}
 \alias{mpe}
-\title{Mean prediction error}
+\title{Mean percentage error}
 \usage{
 mpe(obs, pred)
 }
@@ -15,5 +15,5 @@ mpe(obs, pred)
 A numeric vector
 }
 \description{
-Mean prediction error
+Mean percentage error
 }
diff --git a/man/nm_busulfan.Rd b/man/nm_busulfan.Rd
index a041f12..71ccbd9 100644
--- a/man/nm_busulfan.Rd
+++ b/man/nm_busulfan.Rd
@@ -6,7 +6,7 @@
 \title{Busulfan data}
 \format{
 A NONMEM-style data frame with 240 rows and
-13 columns:
+14 columns:
 
 \describe{
 \item{\code{ID}}{Patient ID.}
diff --git a/man/nrmse.Rd b/man/nrmse.Rd
new file mode 100644
index 0000000..e922a0f
--- /dev/null
+++ b/man/nrmse.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/misc.R
+\name{nrmse}
+\alias{nrmse}
+\title{Normalized root-mean-squared error}
+\usage{
+nrmse(obs, pred)
+}
+\arguments{
+\item{obs}{observations vector}
+
+\item{pred}{predictions vector}
+}
+\value{
+A numeric vector
+}
+\description{
+Normalized root-mean-squared error
+}
diff --git a/man/print.mipdeval_results.Rd b/man/print.mipdeval_results.Rd
new file mode 100644
index 0000000..45de204
--- /dev/null
+++ b/man/print.mipdeval_results.Rd
@@ -0,0 +1,16 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/print_functions.R
+\name{print.mipdeval_results}
+\alias{print.mipdeval_results}
+\title{Print results from run_eval()}
+\usage{
+\method{print}{mipdeval_results}(x, ...)
+}
+\arguments{
+\item{x}{output object from \code{run_eval()}.}
+
+\item{...}{further arguments passed to or from other methods.}
+}
+\description{
+Print results from run_eval()
+}
diff --git a/man/print.mipdeval_results_bayesian_impact.Rd b/man/print.mipdeval_results_bayesian_impact.Rd
new file mode 100644
index 0000000..926bc89
--- /dev/null
+++ b/man/print.mipdeval_results_bayesian_impact.Rd
@@ -0,0 +1,16 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/print_functions.R
+\name{print.mipdeval_results_bayesian_impact}
+\alias{print.mipdeval_results_bayesian_impact}
+\title{Print Bayesian impact results from run_eval()}
+\usage{
+\method{print}{mipdeval_results_bayesian_impact}(x, ...)
+}
+\arguments{
+\item{x}{bayesian impact object in output from \code{run_eval()}.}
+
+\item{...}{further arguments passed to or from other methods.}
+}
+\description{
+Print Bayesian impact results from run_eval()
+}
diff --git a/man/print.mipdeval_results_shrinkage.Rd b/man/print.mipdeval_results_shrinkage.Rd
new file mode 100644
index 0000000..22bd404
--- /dev/null
+++ b/man/print.mipdeval_results_shrinkage.Rd
@@ -0,0 +1,16 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/print_functions.R
+\name{print.mipdeval_results_shrinkage}
+\alias{print.mipdeval_results_shrinkage}
+\title{Print shrinkage results from run_eval()}
+\usage{
+\method{print}{mipdeval_results_shrinkage}(x, ...)
+}
+\arguments{
+\item{x}{shrinkage object in output from \code{run_eval()}.}
+
+\item{...}{further arguments passed to or from other methods.}
+}
+\description{
+Print shrinkage results from run_eval()
+}
diff --git a/man/print.mipdeval_results_stats_summ.Rd b/man/print.mipdeval_results_stats_summ.Rd
new file mode 100644
index 0000000..1807c19
--- /dev/null
+++ b/man/print.mipdeval_results_stats_summ.Rd
@@ -0,0 +1,16 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/print_functions.R
+\name{print.mipdeval_results_stats_summ}
+\alias{print.mipdeval_results_stats_summ}
+\title{Print predictive performance statistics from run_eval()}
+\usage{
+\method{print}{mipdeval_results_stats_summ}(x, ...)
+}
+\arguments{
+\item{x}{stats summary object in output from \code{run_eval()}.}
+
+\item{...}{further arguments passed to or from other methods.}
+}
+\description{
+Print predictive performance statistics from run_eval()
+}
diff --git a/man/run_eval_core.Rd b/man/run_eval_core.Rd
index dff7ba0..cc0195a 100644
--- a/man/run_eval_core.Rd
+++ b/man/run_eval_core.Rd
@@ -8,9 +8,6 @@ an individual's dataset}
 run_eval_core(
   mod_obj,
   data,
-  parameters = NULL,
-  omega = NULL,
-  ruv = NULL,
   weights = NULL,
   weight_prior = 1,
   censor_covariates = TRUE,
@@ -24,13 +21,6 @@ run_eval_core(
 
 \item{data}{NONMEM-style data.frame, or path to CSV file with NONMEM data}
 
-\item{parameters}{list of parameters}
-
-\item{omega}{between subject variability, supplied as vector specifiying the
-lower triangle of the covariance matrix of random effects}
-
-\item{ruv}{residual error variability magnitude, specified as list.}
-
 \item{weights}{vector of weights for error. Length of vector should be same 
 as length of observation vector. If NULL (default), all weights are equal. 
 Used in both MAP and NP methods. Note that `weights` argument will also 

From 08b9d39d8bff7fb0b3abb8b28983abb83a837a95 Mon Sep 17 00:00:00 2001
From: Ron Keizer <ron@insight-rx.com>
Date: Thu, 25 Sep 2025 22:53:18 +0000
Subject: [PATCH 16/22] fix tests

---
 R/run_eval.R                                  |  9 +-
 .../testthat/test-calculate_bayesian_impact.R | 18 ++--
 tests/testthat/test-calculate_shrinkage.R     | 82 ++++++++++---------
 tests/testthat/test-run_eval.R                | 10 +--
 tests/testthat/test-run_eval_grouped.R        |  2 +-
 5 files changed, 62 insertions(+), 59 deletions(-)

diff --git a/R/run_eval.R b/R/run_eval.R
index be6cc20..5ffb758 100644
--- a/R/run_eval.R
+++ b/R/run_eval.R
@@ -61,8 +61,6 @@ run_eval <- function(
   progress = TRUE,
   verbose = TRUE
 ) {
-  # TODO: Refactor to S3 method like in parse_model() to make required and
-  # optional arguments clearer for the different types of `model` inputs.
 
   ## 0. Gather model information in an object
   mod_obj <- parse_model(
@@ -75,7 +73,7 @@ run_eval <- function(
   )
 
   ## 1. read NONMEM data from file or data.frame. Do some simple checks
-  if(verbose) cli::cli_alert_info("Reading and parsing input data")
+  if(verbose) cli::cli_progress_step("Reading and parsing input data")
   input_data <- read_input_data(data) |>
     check_input_data(
       dictionary = dictionary
@@ -105,7 +103,7 @@ run_eval <- function(
   }
 
   ## 3. run the core function on each individual-level dataset in the list
-  if(verbose) cli::cli_alert_info("Running forecasts for subjects in dataset")
+  if(verbose) cli::cli_progress_step("Running forecasts for subjects in dataset")
   if(threads > 1) {
     # TODO: consider using purrr::in_parallel() in the future when it's stable.
     future::plan(future::multisession, workers = threads)
@@ -129,7 +127,7 @@ run_eval <- function(
   }
 
   ## 4. Combine results and basic stats into return object
-  if(verbose) cli::cli_alert_info("Calculating forecasting statistics")
+  if(verbose) cli::cli_progress_step("Calculating forecasting statistics")
   res_tbl <- dplyr::bind_rows(res) |>
     tibble::as_tibble()
   out <- list(
@@ -143,6 +141,5 @@ run_eval <- function(
   out$bayesian_impact <- calculate_bayesian_impact(out)
 
   ## 5. Return results
-  if(verbose) cli::cli_alert_info("Done")
   out
 }
diff --git a/tests/testthat/test-calculate_bayesian_impact.R b/tests/testthat/test-calculate_bayesian_impact.R
index 3d2b2f1..50b3fe7 100644
--- a/tests/testthat/test-calculate_bayesian_impact.R
+++ b/tests/testthat/test-calculate_bayesian_impact.R
@@ -7,7 +7,7 @@ test_that("calculate_bayesian_impact returns correct structure", {
     mape = c(20, 10, 25, 12)
   )
 
-  result <- calculate_bayesian_impact(test_data)
+  result <- calculate_bayesian_impact(list(stats_summ = test_data))
 
   expect_s3_class(result, "data.frame")
   expect_named(result, c("bi_rmse", "bi_mape"))
@@ -23,7 +23,7 @@ test_that("calculate_bayesian_impact calculates correct values", {
     mape = c(20, 10)  # iter_ipred is 50% better
   )
 
-  result <- calculate_bayesian_impact(test_data)
+  result <- calculate_bayesian_impact(list(stats_summ = test_data))
 
   # BI = 100 * (1 - (iter_ipred / pred))
   # For RMSE: 100 * (1 - (5/10)) = 50%
@@ -41,7 +41,7 @@ test_that("calculate_bayesian_impact handles negative impact correctly", {
     mape = c(10, 20)   # iter_ipred is worse
   )
 
-  result <- calculate_bayesian_impact(test_data)
+  result <- calculate_bayesian_impact(list(stats_summ = test_data))
 
   # BI = 100 * (1 - (10/5)) = -100%
   expect_equal(result$bi_rmse, -100)
@@ -57,7 +57,7 @@ test_that("calculate_bayesian_impact filters apriori correctly", {
     mape = c(200, 100, 20, 10)
   )
 
-  result <- calculate_bayesian_impact(test_data)
+  result <- calculate_bayesian_impact(list(stats_summ = test_data))
 
   # Should only use the apriori = FALSE rows (10, 5 and 20, 10)
   expect_equal(result$bi_rmse, 50)
@@ -73,7 +73,7 @@ test_that("calculate_bayesian_impact handles edge cases", {
     mape = c(20, 20)
   )
 
-  result <- calculate_bayesian_impact(test_data)
+  result <- calculate_bayesian_impact(list(stats_summ = test_data))
 
   expect_equal(result$bi_rmse, 0)
   expect_equal(result$bi_mape, 0)
@@ -90,7 +90,9 @@ test_that("calculate_bayesian_impact handles division by zero", {
 
   # This should handle division by zero gracefully. No warning expected, but Inf
   # values are possible.
-  result <- expect_no_warning(calculate_bayesian_impact(test_data))
+  result <- expect_no_warning(
+    calculate_bayesian_impact(list(stats_summ = test_data))
+  )
   expect_true(is.infinite(result$bi_rmse) || is.nan(result$bi_rmse))
   expect_true(is.infinite(result$bi_mape) || is.nan(result$bi_mape))
 })
@@ -103,5 +105,7 @@ test_that("calculate_bayesian_impact requires correct column names", {
     # Missing apriori and mape columns
   )
 
-  expect_error(calculate_bayesian_impact(incomplete_data))
+  expect_error(
+    calculate_bayesian_impact(list(stats_summ = test_data))
+  )
 })
diff --git a/tests/testthat/test-calculate_shrinkage.R b/tests/testthat/test-calculate_shrinkage.R
index 3788036..8f1a97a 100644
--- a/tests/testthat/test-calculate_shrinkage.R
+++ b/tests/testthat/test-calculate_shrinkage.R
@@ -5,7 +5,7 @@ test_that("calculate_shrinkage returns correct structure", {
     fixed = character(0),  # No fixed parameters
     omega = c(0.1, 0.2, 0.15)  # Triangle format: var(cl), var(v), var(ka)
   )
-  
+
   # Mock res_df with parameter estimates
   res_df <- data.frame(
     `_iteration` = c(1, 1, 2, 2),
@@ -14,9 +14,9 @@ test_that("calculate_shrinkage returns correct structure", {
     ka = c(0.7, 0.9, 0.6, 1.1),
     check.names = FALSE
   )
-  
-  result <- calculate_shrinkage(res_df, mod_obj)
-  
+
+  result <- calculate_shrinkage(list(results = res_df, mod_obj = mod_obj))
+
   expect_s3_class(result, "data.frame")
   expect_true("_iteration" %in% names(result))
   expect_true(all(c("cl", "v", "ka") %in% names(result)))
@@ -29,7 +29,7 @@ test_that("calculate_shrinkage handles exponential parameters correctly", {
     fixed = character(0),
     omega = c(0.09, 0.01, 0.16)  # var(cl) = 0.09, var(v) = 0.16
   )
-  
+
   # Create data where individual estimates vary from population
   res_df <- data.frame(
     `_iteration` = c(1, 1, 1),
@@ -37,9 +37,9 @@ test_that("calculate_shrinkage handles exponential parameters correctly", {
     v = c(15.0, 16.5, 13.5),  # Around pop value of 15.0
     check.names = FALSE
   )
-  
-  result <- calculate_shrinkage(res_df, mod_obj)
-  
+
+  result <- calculate_shrinkage(list(results = res_df, mod_obj = mod_obj))
+
   expect_equal(nrow(result), 1)
   expect_true(is.numeric(result$cl))
   expect_true(is.numeric(result$v))
@@ -53,16 +53,16 @@ test_that("calculate_shrinkage handles kappa parameters (additive IIV)", {
     fixed = character(0),
     omega = c(0.01, 0.0025)  # var(kappa_cl), var(kappa_v)
   )
-  
+
   res_df <- data.frame(
     `_iteration` = c(1, 1, 1),
     kappa_cl = c(0.12, 0.08, 0.11),
     kappa_v = c(0.06, 0.04, 0.055),
     check.names = FALSE
   )
-  
-  result <- calculate_shrinkage(res_df, mod_obj)
-  
+
+  result <- calculate_shrinkage(list(results = res_df, mod_obj = mod_obj))
+
   expect_equal(nrow(result), 1)
   expect_true(is.numeric(result$kappa_cl))
   expect_true(is.numeric(result$kappa_v))
@@ -74,16 +74,16 @@ test_that("calculate_shrinkage handles zero population parameter", {
     fixed = character(0),
     omega = c(0.1, 0.01, 0.05)
   )
-  
+
   res_df <- data.frame(
     `_iteration` = c(1, 1, 1),
     cl = c(1.2, 1.8, 1.0),
     baseline = c(-0.1, 0.1, 0.05),  # Should use additive for zero pop
     check.names = FALSE
   )
-  
-  result <- calculate_shrinkage(res_df, mod_obj)
-  
+
+  result <- calculate_shrinkage(list(results = res_df, mod_obj = mod_obj))
+
   expect_equal(nrow(result), 1)
   expect_true(is.numeric(result$cl))
   expect_true(is.numeric(result$baseline))
@@ -95,7 +95,7 @@ test_that("calculate_shrinkage handles fixed parameters correctly", {
     fixed = c("ka"),  # ka is fixed
     omega = c(0.1, 0.01, 0.2)  # Only cl and v have random effects
   )
-  
+
   res_df <- data.frame(
     `_iteration` = c(1, 1, 2, 2),
     cl = c(1.2, 1.8, 1.0, 2.0),
@@ -103,9 +103,9 @@ test_that("calculate_shrinkage handles fixed parameters correctly", {
     ka = c(0.8, 0.8, 0.8, 0.8),  # Fixed value
     check.names = FALSE
   )
-  
-  result <- calculate_shrinkage(res_df, mod_obj)
-  
+
+  result <- calculate_shrinkage(list(results = res_df, mod_obj = mod_obj))
+
   expect_true("cl" %in% names(result))
   expect_true("v" %in% names(result))
   expect_false("ka" %in% names(result))  # Should not appear (fixed)
@@ -117,31 +117,31 @@ test_that("calculate_shrinkage handles multiple iterations", {
     fixed = character(0),
     omega = c(0.1, 0.01, 0.2)
   )
-  
+
   res_df <- data.frame(
     `_iteration` = c(1, 1, 1, 2, 2, 2, 3, 3, 3),
     cl = c(1.2, 1.8, 1.0, 1.6, 1.4, 1.7, 1.1, 1.9, 1.3),
     v = c(9.5, 11.0, 8.8, 10.8, 9.2, 11.5, 8.9, 12.0, 9.8),
     check.names = FALSE
   )
-  
-  result <- calculate_shrinkage(res_df, mod_obj)
-  
+
+  result <- calculate_shrinkage(list(results = res_df, mod_obj = mod_obj))
+
   expect_equal(nrow(result), 3)  # Three iterations
   expect_equal(result$`_iteration`, c(1, 2, 3))
 })
 
 test_that("calc_eta helper function works correctly", {
   parameters <- list(cl = 2.0, v = 15.0, kappa_cl = 0.1, baseline = 0)
-  
+
   # Test exponential parameter (cl)
   eta_cl <- calc_eta(2.2, "cl", parameters)
   expect_equal(eta_cl, log(2.2/2.0), tolerance = 1e-10)
-  
+
   # Test kappa parameter (additive)
   eta_kappa <- calc_eta(0.12, "kappa_cl", parameters)
   expect_equal(eta_kappa, 0.12 - 0.1)
-  
+
   # Test zero population parameter (additive)
   eta_baseline <- calc_eta(0.05, "baseline", parameters)
   expect_equal(eta_baseline, 0.05 - 0)
@@ -151,10 +151,10 @@ test_that("calc_shrinkage helper function works correctly", {
   # Test shrinkage calculation
   omega_list <- list(cl = 0.09)  # omega^2 = 0.09, so omega = 0.3
   eta_values <- c(0.1, -0.1, 0.05, -0.05, 0.02)
-  
+
   shrinkage <- calc_shrinkage(eta_values, "eta_cl", omega_list)
-  expected_shrinkage <- 100 * (1 - (sd(eta_values) / sqrt(0.09)))
-  
+  expected_shrinkage <- round(100 * (1 - (sd(eta_values) / sqrt(0.09))), 1)
+
   expect_equal(shrinkage, expected_shrinkage, tolerance = 1e-10)
   expect_true(shrinkage >= -100 && shrinkage <= 100)  # Reasonable range
 })
@@ -165,12 +165,12 @@ test_that("get_omega_for_parameters works correctly", {
     fixed = c("f1"),  # f1 is fixed
     omega = c(
       0.1,
-      0.05, 0.02, 
+      0.05, 0.02,
       0.15, 0.01, 0.25)  # Triangle: 3x3 matrix for cl, v, ka
   )
-  
+
   omega_list <- get_omega_for_parameters(mod_obj)
-  
+
   expect_equal(omega_list, list(cl = 0.1, v = 0.02, ka = 0.25))
 })
 
@@ -180,16 +180,16 @@ test_that("calculate_shrinkage handles edge case with no variability", {
     fixed = character(0),
     omega = c(0.1)
   )
-  
+
   # All individual estimates are the same (no variability)
   res_df <- data.frame(
     `_iteration` = c(1, 1, 1),
     cl = c(1.5, 1.5, 1.5),  # No variation
     check.names = FALSE
   )
-  
-  result <- calculate_shrinkage(res_df, mod_obj)
-  
+
+  result <- calculate_shrinkage(list(results = res_df, mod_obj = mod_obj))
+
   expect_equal(nrow(result), 1)
   expect_equal(result$cl, 100, tolerance = 1e-10)  # Perfect shrinkage
 })
@@ -200,9 +200,11 @@ test_that("calculate_shrinkage validates input structure", {
     fixed = character(0),
     omega = c(0.1)
   )
-  
+
   # Missing _iteration column
   bad_res_df <- data.frame(cl = c(1.2, 1.8, 1.0))
-  
-  expect_error(calculate_shrinkage(bad_res_df, mod_obj))
+
+  expect_error(
+    calculate_shrinkage(list(results = bad_res_df, mod_obj = mod_obj))
+  )
 })
diff --git a/tests/testthat/test-run_eval.R b/tests/testthat/test-run_eval.R
index 6aae168..cddf91e 100644
--- a/tests/testthat/test-run_eval.R
+++ b/tests/testthat/test-run_eval.R
@@ -15,9 +15,9 @@ test_that("Basic run with vanco data + model works", {
   )
 
   ## Expected structure:
-  expect_equal(names(res), c("results", "stats", "shrinkage", "bayesian_impact"))
+  expect_equal(names(res), c("results", "mod_obj", "data", "stats_summ", "shrinkage", "bayesian_impact"))
   expect_s3_class(res$results, c("tbl_df", "tbl", "data.frame"))
-  expect_s3_class(res$stats, c("tbl_df", "tbl", "data.frame"))
+  expect_s3_class(res$stats_summ, c("tbl_df", "tbl", "data.frame"))
 
   # Using PKPDsim model library:
   res_2 <- run_eval(
@@ -60,7 +60,7 @@ test_that("Run also works when `model` argument just references the package", {
     progress = F,
     ids = c(1:3)
   )
-  expect_equal(names(res), c("results", "stats", "shrinkage", "bayesian_impact"))
+  expect_equal(names(res), c("results", "mod_obj", "data", "stats_summ", "shrinkage", "bayesian_impact"))
   # TODO: test outputs
 })
 
@@ -79,7 +79,7 @@ test_that("Flattening of prior results in different predictions", {
     progress = F,
     ids = c(1:3)
   )
-  expect_equal(names(res_flat), c("results", "stats", "shrinkage", "bayesian_impact"))
+  expect_equal(names(res), c("results", "mod_obj", "data", "stats_summ", "shrinkage", "bayesian_impact"))
   expect_true(
     all(
       res$results |>
@@ -104,6 +104,6 @@ test_that("Run also works when `dictionary` is used", {
     progress = F,
     ids = 1
   )
-  expect_equal(names(res), c("results", "stats", "shrinkage", "bayesian_impact"))
+  expect_equal(names(res), c("results", "mod_obj", "data", "stats_summ", "shrinkage", "bayesian_impact"))
 })
 
diff --git a/tests/testthat/test-run_eval_grouped.R b/tests/testthat/test-run_eval_grouped.R
index 3982cb7..10a5ebf 100644
--- a/tests/testthat/test-run_eval_grouped.R
+++ b/tests/testthat/test-run_eval_grouped.R
@@ -13,7 +13,7 @@ test_that("Grouped run for busulfan/shukla works", {
     progress = FALSE,
     ids = c(1:n_ids)
   )
-  expect_equal(names(res), c("results", "stats", "shrinkage", "bayesian_impact"))
+  expect_equal(names(res), c("results", "mod_obj", "data", "stats_summ", "shrinkage", "bayesian_impact"))
 
   ## Reference results from PsN proseval:
   raw_proseval <- read.csv(

From 3e0e9e4765d6f2f793279bbf0c8d3dc3d130a862 Mon Sep 17 00:00:00 2001
From: Ron Keizer <ron@insight-rx.com>
Date: Fri, 26 Sep 2025 21:36:48 +0000
Subject: [PATCH 17/22] add incremental option

---
 R/run_eval.R                   |  2 ++
 R/run_eval_core.R              | 46 +++++++++++++++++++++++++++++++---
 man/nm_busulfan.Rd             |  2 +-
 tests/testthat/test-run_eval.R | 12 +++++++++
 4 files changed, 58 insertions(+), 4 deletions(-)

diff --git a/R/run_eval.R b/R/run_eval.R
index 5ffb758..d128199 100644
--- a/R/run_eval.R
+++ b/R/run_eval.R
@@ -113,6 +113,7 @@ run_eval <- function(
       mod_obj = mod_obj,
       censor_covariates = censor_covariates,
       weight_prior = weight_prior,
+      incremental = incremental,
       progress_function = p
     )
   } else {
@@ -122,6 +123,7 @@ run_eval <- function(
       mod_obj = mod_obj,
       censor_covariates = censor_covariates,
       weight_prior = weight_prior,
+      incremental = incremental,
       progress_function = p
     )
   }
diff --git a/R/run_eval_core.R b/R/run_eval_core.R
index b1ef6ac..01259be 100644
--- a/R/run_eval_core.R
+++ b/R/run_eval_core.R
@@ -42,10 +42,17 @@ run_eval_core <- function(
     ## Do a fit with PKPDmap, using only the points
     ## with weight=1. The rest of the points will get a predicted value
     ## in the output, but they were not weighted in the fit.
+    mod_upd <- mod_obj
+    if(i > 1) {
+      if(incremental) {
+        mod_upd$parameters <- fit$parameters # take params from previous fit
+        mod_upd$omega <- fit$vcov
+      }
+    }
     fit <- PKPDmap::get_map_estimates(
       model = mod_obj$model,
-      parameters = mod_obj$parameters,
-      omega = mod_obj$omega,
+      parameters = mod_upd$parameters,
+      omega = mod_upd$omega,
       error = mod_obj$ruv,
       fixed = mod_obj$fixed,
       as_eta = mod_obj$kappa,
@@ -76,6 +83,39 @@ run_eval_core <- function(
     )
   }
 
+  ## Get data for MAP fit
+  if(incremental) {
+    ## need to do an extra MAP Bayesian fit, because we can't use the
+    ## incrementally updated parameters + omega
+    fit_map <- PKPDmap::get_map_estimates(
+      model = mod_obj$model,
+      parameters = mod_obj$parameters, # use original model params!
+      omega = mod_obj$omega, # use original model params!
+      error = mod_obj$ruv,
+      fixed = mod_obj$fixed,
+      as_eta = mod_obj$kappa,
+      data = data$observations,
+      covariates = cov_data,
+      regimen = data$regimen,
+      weight_prior = weight_prior,
+      weights = weights,
+      iov_bins = mod_obj$bins,
+      verbose = FALSE
+    )
+    map_pred_data <- tibble::tibble(
+      id = obs_data$id,
+      t = obs_data$t,
+      dv = fit_map$dv,
+      ipred = fit_map$ipred,
+      pred = fit_map$pred,
+      `_iteration` = iterations[i],
+      `_grouper` = obs_data$`_grouper`
+    )
+  } else {  # just take last fit object and pred_data
+    fit_map <- fit
+    map_pred_data <- pred_data
+  }
+
   ## 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()`.
@@ -96,7 +136,7 @@ run_eval_core <- function(
     dplyr::filter(.data$`_iteration` == (.data$`_grouper` - 1)) |>
     dplyr::mutate(
       iter_ipred = .data$ipred,
-      map_ipred = pred_data$ipred, # ipred from last fit (full MAP)
+      map_ipred = map_pred_data$ipred, # ipred from full retrospective MAP
       apriori = (.data$`_iteration` == 0)
     ) |>
     dplyr::select(
diff --git a/man/nm_busulfan.Rd b/man/nm_busulfan.Rd
index 71ccbd9..a041f12 100644
--- a/man/nm_busulfan.Rd
+++ b/man/nm_busulfan.Rd
@@ -6,7 +6,7 @@
 \title{Busulfan data}
 \format{
 A NONMEM-style data frame with 240 rows and
-14 columns:
+13 columns:
 
 \describe{
 \item{\code{ID}}{Patient ID.}
diff --git a/tests/testthat/test-run_eval.R b/tests/testthat/test-run_eval.R
index cddf91e..9a3e838 100644
--- a/tests/testthat/test-run_eval.R
+++ b/tests/testthat/test-run_eval.R
@@ -107,3 +107,15 @@ test_that("Run also works when `dictionary` is used", {
   expect_equal(names(res), c("results", "mod_obj", "data", "stats_summ", "shrinkage", "bayesian_impact"))
 })
 
+test_that("Incremental Bayes method works", {
+  local_mipdeval_options()
+  res <- run_eval(
+    model = "pkvancothomson",
+    data = nm_vanco,
+    incremental = TRUE,
+    progress = F,
+    ids = c(1:3)
+  )
+  expect_equal(names(res), c("results", "mod_obj", "data", "stats_summ", "shrinkage", "bayesian_impact"))
+  # TODO: test outputs
+})

From 8339566e8ad9ed99bb0cf57d39433901a1289eed Mon Sep 17 00:00:00 2001
From: Ron Keizer <ron@insight-rx.com>
Date: Fri, 26 Sep 2025 23:11:28 +0000
Subject: [PATCH 18/22] split off weighting function

---
 R/run_eval_core.R              | 50 +++++++++++++++++++++++++++-------
 man/handle_sample_weighting.Rd | 29 ++++++++++++++++++++
 tests/testthat/test-run_eval.R |  2 +-
 3 files changed, 70 insertions(+), 11 deletions(-)
 create mode 100644 man/handle_sample_weighting.Rd

diff --git a/R/run_eval_core.R b/R/run_eval_core.R
index 01259be..cd6fcdc 100644
--- a/R/run_eval_core.R
+++ b/R/run_eval_core.R
@@ -26,10 +26,15 @@ run_eval_core <- function(
 
   for(i in seq_along(iterations)) {
 
-    ## TODO: this will change once we support sample-weighting strategies
-    ##       For now this is just simple full-weighting for all previous points
-    weights <- rep(0, nrow(obs_data))
-    weights[obs_data[["_grouper"]] %in% iterations[1:i]] <- 1
+    ## Select which samples should be used in fit, for regular iterative
+    ## forecasting and incremental.
+    ## TODO: handle weighting down of earlier samples
+    weights <- handle_sample_weighting(
+      obs_data,
+      iterations,
+      incremental,
+      i
+    )
 
     ## Should covariate data be leaked? PsN::proseval does this,
     ## but for use in MIPD they should be censored.
@@ -43,11 +48,9 @@ run_eval_core <- function(
     ## with weight=1. The rest of the points will get a predicted value
     ## in the output, but they were not weighted in the fit.
     mod_upd <- mod_obj
-    if(i > 1) {
-      if(incremental) {
-        mod_upd$parameters <- fit$parameters # take params from previous fit
-        mod_upd$omega <- fit$vcov
-      }
+    if(incremental & i > 1) {
+      mod_upd$parameters <- fit$parameters # take params from previous fit
+      mod_upd$omega <- fit$vcov
     }
     fit <- PKPDmap::get_map_estimates(
       model = mod_obj$model,
@@ -98,7 +101,7 @@ run_eval_core <- function(
       covariates = cov_data,
       regimen = data$regimen,
       weight_prior = weight_prior,
-      weights = weights,
+      weights = NULL, # no sample weighting, full MAP Bayesian on all samples
       iov_bins = mod_obj$bins,
       verbose = FALSE
     )
@@ -176,3 +179,30 @@ handle_covariate_censoring <- function(
   }
   cov_data
 }
+
+#' Handle weighting of samples
+#'
+#' This function is used to select the samples used in the fit (1 or 0),
+#' but also to select their weight, if a sample weighting strategy is
+#' selected.
+#'
+#' @inheritParams run_eval_core
+#' @param iterations numeric vector of groups
+#' @param i index
+#'
+#' @returns vector of weights (numeric)
+#'
+handle_sample_weighting <- function(
+  obs_data,
+  iterations,
+  incremental,
+  i
+) {
+  weights <- rep(0, nrow(obs_data))
+  if(incremental) { # just fit current sample or group
+    weights[obs_data[["_grouper"]] %in% iterations[i]] <- 1
+  } else { # fit all samples up until current sample
+    weights[obs_data[["_grouper"]] %in% iterations[1:i]] <- 1
+  }
+  weights
+}
diff --git a/man/handle_sample_weighting.Rd b/man/handle_sample_weighting.Rd
new file mode 100644
index 0000000..5c22293
--- /dev/null
+++ b/man/handle_sample_weighting.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/run_eval_core.R
+\name{handle_sample_weighting}
+\alias{handle_sample_weighting}
+\title{Handle weighting of samples}
+\usage{
+handle_sample_weighting(obs_data, iterations, incremental, i)
+}
+\arguments{
+\item{iterations}{numeric vector of groups}
+
+\item{incremental}{should MAP Bayesian do incremental fits in the iterative
+loop. I.e. in this case it would use the first iterations MAP Bayesian
+estimates as input for the second iteration, and so forth. The uncertainty
+around the MAP estimates would be used as the new \code{omega} matrix. This
+approach has been called "model predictive control (MPC)"
+(www.page-meeting.org/?abstract=9076) and may be more predictive than
+"regular" MAP in some scenarios. Default is \code{FALSE}.}
+
+\item{i}{index}
+}
+\value{
+vector of weights (numeric)
+}
+\description{
+This function is used to select the samples used in the fit (1 or 0),
+but also to select their weight, if a sample weighting strategy is
+selected.
+}
diff --git a/tests/testthat/test-run_eval.R b/tests/testthat/test-run_eval.R
index 9a3e838..0642b3e 100644
--- a/tests/testthat/test-run_eval.R
+++ b/tests/testthat/test-run_eval.R
@@ -114,7 +114,7 @@ test_that("Incremental Bayes method works", {
     data = nm_vanco,
     incremental = TRUE,
     progress = F,
-    ids = c(1:3)
+    ids = c(1:10)
   )
   expect_equal(names(res), c("results", "mod_obj", "data", "stats_summ", "shrinkage", "bayesian_impact"))
   # TODO: test outputs

From dff0780f6247dca73df0689bd3a803e788ab807d Mon Sep 17 00:00:00 2001
From: Ron Keizer <ron@insight-rx.com>
Date: Tue, 30 Sep 2025 19:48:39 +0000
Subject: [PATCH 19/22] fix docs + print function

---
 R/print_functions.R            | 4 ++++
 R/run_eval_core.R              | 1 +
 man/handle_sample_weighting.Rd | 2 ++
 tests/testthat/test-run_eval.R | 2 +-
 4 files changed, 8 insertions(+), 1 deletion(-)

diff --git a/R/print_functions.R b/R/print_functions.R
index e15dd2b..6a4d989 100644
--- a/R/print_functions.R
+++ b/R/print_functions.R
@@ -5,6 +5,10 @@
 #'
 #' @export
 print.mipdeval_results <- function(x, ...) {
+  if(is.null(x$results)) {
+    cli::cli_alert_info("No forecasting info available, did you run with `vpc_only=TRUE`?")
+    return(invisible())
+  }
   print(x$stats_summ, ...)
   print(x$shrinkage, ...)
   print(x$bayesian_impact, ...)
diff --git a/R/run_eval_core.R b/R/run_eval_core.R
index cd6fcdc..e31fd64 100644
--- a/R/run_eval_core.R
+++ b/R/run_eval_core.R
@@ -187,6 +187,7 @@ handle_covariate_censoring <- function(
 #' selected.
 #'
 #' @inheritParams run_eval_core
+#' @param obs_data tibble or data.frame with observed data for individual
 #' @param iterations numeric vector of groups
 #' @param i index
 #'
diff --git a/man/handle_sample_weighting.Rd b/man/handle_sample_weighting.Rd
index 3476549..40e009d 100644
--- a/man/handle_sample_weighting.Rd
+++ b/man/handle_sample_weighting.Rd
@@ -7,6 +7,8 @@
 handle_sample_weighting(obs_data, iterations, incremental, i)
 }
 \arguments{
+\item{obs_data}{tibble or data.frame with observed data for individual}
+
 \item{iterations}{numeric vector of groups}
 
 \item{incremental}{should MAP Bayesian do incremental fits in the iterative
diff --git a/tests/testthat/test-run_eval.R b/tests/testthat/test-run_eval.R
index e6e7a3f..fa40762 100644
--- a/tests/testthat/test-run_eval.R
+++ b/tests/testthat/test-run_eval.R
@@ -122,6 +122,6 @@ test_that("Incremental Bayes method works", {
     progress = F,
     ids = c(1:10)
   )
-  expect_equal(names(res), c("results", "mod_obj", "data", "stats_summ", "shrinkage", "bayesian_impact"))
+  expect_equal(names(res), c("results", "mod_obj", "data", "sim", "stats_summ", "shrinkage", "bayesian_impact"))
   # TODO: test outputs
 })

From 689d7411a0bb8edd383b40881d0237bc4ce704d3 Mon Sep 17 00:00:00 2001
From: Ron Keizer <ron@insight-rx.com>
Date: Thu, 9 Oct 2025 21:13:15 +0000
Subject: [PATCH 20/22] add more testing

---
 tests/testthat/test-run_eval.R | 28 ++++++++++++++++++++++++++--
 1 file changed, 26 insertions(+), 2 deletions(-)

diff --git a/tests/testthat/test-run_eval.R b/tests/testthat/test-run_eval.R
index fa40762..7721b17 100644
--- a/tests/testthat/test-run_eval.R
+++ b/tests/testthat/test-run_eval.R
@@ -19,6 +19,10 @@ test_that("Basic run with vanco data + model works", {
   expect_equal(names(res), c("results", "mod_obj", "data", "sim", "stats_summ", "shrinkage", "bayesian_impact"))
   expect_s3_class(res$results, c("tbl_df", "tbl", "data.frame"))
   expect_s3_class(res$stats_summ, c("tbl_df", "tbl", "data.frame"))
+  expect_equal(
+    round(res$results$CL[1:5], 3),
+    c(2.99, 2.685, 2.462, 2.430, 2.439)
+  )
 
   # Using PKPDsim model library:
   res_2 <- run_eval(
@@ -120,8 +124,28 @@ test_that("Incremental Bayes method works", {
     data = nm_vanco,
     incremental = TRUE,
     progress = F,
-    ids = c(1:10)
+    ids = c(1:3)
   )
   expect_equal(names(res), c("results", "mod_obj", "data", "sim", "stats_summ", "shrinkage", "bayesian_impact"))
-  # TODO: test outputs
+  expect_equal(
+    res$stats_summ,
+    structure(list(
+      type = c("iter_ipred", "iter_ipred", "map_ipred",
+"map_ipred", "pred", "pred"),
+      apriori = c(0, 1, 0, 1, 0, 1),
+      rmse = c(4.051, 8.299, 1.96, 4.319, 4.789, 8.299),
+      nrmse = c(0.411,
+    0.232, 0.199, 0.121, 0.486, 0.232),
+      mpe = c(-0.339, -0.004,
+    -0.045, -0.002, -0.415, -0.004),
+      mape = c(0.445, 0.246, 0.166,
+    0.118, 0.558, 0.246)
+    ),
+    class = c("mipdeval_results_stats_summ","tbl_df", "tbl", "data.frame"),
+    row.names = c(NA, -6L))
+  )
+  expect_equal(
+    round(res$results$CL[1:5], 3),
+    c(2.99, 2.685, 2.529, 2.668, 2.666) # should be different CL progression (expect first 2 values) than from test above (non-incremental)
+  )
 })

From db610f39b391562f45cd9b6246d68e7b44cfd93b Mon Sep 17 00:00:00 2001
From: Ron Keizer <ron@insight-rx.com>
Date: Thu, 9 Oct 2025 21:20:04 +0000
Subject: [PATCH 21/22] add tests for handle_sample_weighting()

---
 tests/testthat/test-handle_sample_weighting.R | 111 ++++++++++++++++++
 1 file changed, 111 insertions(+)
 create mode 100644 tests/testthat/test-handle_sample_weighting.R

diff --git a/tests/testthat/test-handle_sample_weighting.R b/tests/testthat/test-handle_sample_weighting.R
new file mode 100644
index 0000000..27921f5
--- /dev/null
+++ b/tests/testthat/test-handle_sample_weighting.R
@@ -0,0 +1,111 @@
+test_that("handle_sample_weighting returns correct structure", {
+
+  # Create test data
+  obs_data <- data.frame(
+    id = c(1, 1, 1, 1),
+    t = c(0, 12, 24, 48),
+    dv = c(10, 8, 6, 4),
+    `_grouper` = c(1, 2, 3, 4),
+    check.names = FALSE
+  )
+
+  iterations <- c(1, 2, 3, 4)
+
+  result <- handle_sample_weighting(obs_data, iterations, incremental = FALSE, i = 2)
+
+  expect_type(result, "double")
+  expect_length(result, nrow(obs_data))
+  expect_true(all(result %in% c(0, 1)))  # Weights should be 0 or 1
+})
+
+test_that("handle_sample_weighting works correctly for non-incremental mode", {
+  obs_data <- data.frame(
+    `_grouper` = c(1, 2, 3, 4, 5),
+    check.names = FALSE
+  )
+
+  iterations <- c(1, 2, 3, 4, 5)
+
+  # Test i = 1 (first iteration)
+  weights_i1 <- handle_sample_weighting(obs_data, iterations, incremental = FALSE, i = 1)
+  expect_equal(weights_i1, c(1, 0, 0, 0, 0))
+
+  # Test i = 3 (third iteration)
+  weights_i3 <- handle_sample_weighting(obs_data, iterations, incremental = FALSE, i = 3)
+  expect_equal(weights_i3, c(1, 1, 1, 0, 0))
+
+  # Test i = 5 (last iteration)
+  weights_i5 <- handle_sample_weighting(obs_data, iterations, incremental = FALSE, i = 5)
+  expect_equal(weights_i5, c(1, 1, 1, 1, 1))
+})
+
+test_that("handle_sample_weighting works correctly for incremental mode", {
+  obs_data <- data.frame(
+    `_grouper` = c(1, 2, 3, 4, 5),
+    check.names = FALSE
+  )
+
+  iterations <- c(1, 2, 3, 4, 5)
+
+  # Test i = 1 (first iteration) - should only use first group
+  weights_i1 <- handle_sample_weighting(obs_data, iterations, incremental = TRUE, i = 1)
+  expect_equal(weights_i1, c(1, 0, 0, 0, 0))
+
+  # Test i = 3 (third iteration) - should only use third group
+  weights_i3 <- handle_sample_weighting(obs_data, iterations, incremental = TRUE, i = 3)
+  expect_equal(weights_i3, c(0, 0, 1, 0, 0))
+
+  # Test i = 5 (last iteration) - should only use fifth group
+  weights_i5 <- handle_sample_weighting(obs_data, iterations, incremental = TRUE, i = 5)
+  expect_equal(weights_i5, c(0, 0, 0, 0, 1))
+})
+
+test_that("handle_sample_weighting handles multiple observations per group", {
+  # Multiple observations can have the same _grouper value
+  obs_data <- data.frame(
+    `_grouper` = c(1, 1, 2, 2, 3, 3),
+    check.names = FALSE
+  )
+
+  iterations <- c(1, 2, 3)
+
+  # Non-incremental: i = 2 should include groups 1 and 2
+  weights_non_inc <- handle_sample_weighting(obs_data, iterations, incremental = FALSE, i = 2)
+  expect_equal(weights_non_inc, c(1, 1, 1, 1, 0, 0))
+
+  # Incremental: i = 2 should only include group 2
+  weights_inc <- handle_sample_weighting(obs_data, iterations, incremental = TRUE, i = 2)
+  expect_equal(weights_inc, c(0, 0, 1, 1, 0, 0))
+})
+
+
+test_that("handle_sample_weighting handles edge case with single observation", {
+  obs_data <- data.frame(
+    `_grouper` = c(1),
+    check.names = FALSE
+  )
+
+  iterations <- c(1)
+
+  # Both modes should give same result for single observation
+  weights_non_inc <- handle_sample_weighting(obs_data, iterations, incremental = FALSE, i = 1)
+  weights_inc <- handle_sample_weighting(obs_data, iterations, incremental = TRUE, i = 1)
+
+  expect_equal(weights_non_inc, c(1))
+  expect_equal(weights_inc, c(1))
+})
+
+test_that("handle_sample_weighting maintains correct vector length", {
+  # Test with various data sizes
+  for(n_obs in c(1, 5, 10, 100)) {
+    obs_data <- data.frame(
+      `_grouper` = rep(1:min(n_obs, 10), length.out = n_obs),
+      check.names = FALSE
+    )
+
+    iterations <- unique(obs_data[["_grouper"]])
+
+    weights <- handle_sample_weighting(obs_data, iterations, incremental = FALSE, i = 1)
+    expect_length(weights, n_obs)
+  }
+})

From 4f53158976529985ecb8187795e99fe513aebd57 Mon Sep 17 00:00:00 2001
From: roninsightrx <ron@insight-rx.com>
Date: Fri, 10 Oct 2025 12:13:51 -0700
Subject: [PATCH 22/22] Add more statistics to the output table (#28)

* add ofv and ss

* add test

* Update R/run_eval_core.R

Co-authored-by: Michael McCarthy <51542091+mccarthy-m-g@users.noreply.github.com>

---------

Co-authored-by: Michael McCarthy <51542091+mccarthy-m-g@users.noreply.github.com>
---
 R/misc.R                       | 16 ++++++++++++++++
 R/run_eval_core.R              | 10 ++++++++--
 man/ss.Rd                      | 18 ++++++++++++++++++
 tests/testthat/test-run_eval.R | 14 ++++++++++++++
 4 files changed, 56 insertions(+), 2 deletions(-)
 create mode 100644 man/ss.Rd

diff --git a/R/misc.R b/R/misc.R
index ac76d0e..d0a1ab0 100644
--- a/R/misc.R
+++ b/R/misc.R
@@ -67,3 +67,19 @@ mape <- function (obs, pred) {
 mpe <- function (obs, pred) {
   sum((obs - pred)/obs)/length(obs)
 }
+
+#' Weighted sum-of-squares of residuals
+#'
+#' @inheritParams rmse
+#' @param w weights
+#'
+ss <- function(obs, pred, w = NULL) {
+  if(is.null(w)) {
+    w <- rep(1, length(obs))
+  }
+  if (length(obs) != length(pred) || length(obs) != length(w)) {
+    cli::cli_abort("`obs`, `pred`, and `w` must have the same length")
+  }
+  if(sum(w) == 0) return(NA)
+  sum(w * (obs - pred)^2)
+}
diff --git a/R/run_eval_core.R b/R/run_eval_core.R
index e31fd64..08bbdf9 100644
--- a/R/run_eval_core.R
+++ b/R/run_eval_core.R
@@ -75,6 +75,8 @@ run_eval_core <- function(
       dv = fit$dv,
       ipred = fit$ipred,
       pred = fit$pred,
+      ofv = fit$fit$value,
+      ss_w = ss(fit$dv, fit$ipred, weights),
       `_iteration` = iterations[i],
       `_grouper` = obs_data$`_grouper`
     )
@@ -111,6 +113,8 @@ run_eval_core <- function(
       dv = fit_map$dv,
       ipred = fit_map$ipred,
       pred = fit_map$pred,
+      ofv = fit_map$fit$value,
+      ss_w = ss(fit_map$dv, fit_map$ipred, w = NULL),
       `_iteration` = iterations[i],
       `_grouper` = obs_data$`_grouper`
     )
@@ -127,7 +131,9 @@ run_eval_core <- function(
       dplyr::filter(.data$`_iteration` == 1) |>
       dplyr::mutate(
         `_iteration` = 0,
-        ipred = .data$pred
+        ipred = .data$pred,
+        ofv = NA,
+        ss_w = NA
       ) |> # set to population parameters, not individual estimates
         dplyr::select(-!!names(mod_obj$parameters)) |>
         dplyr::left_join(
@@ -144,12 +150,12 @@ run_eval_core <- function(
     ) |>
     dplyr::select(
       "id", "_iteration", "_grouper", "t", "dv", "pred", "map_ipred",
+      "ofv", "ss_w",
       "iter_ipred", "apriori",
       !!names(mod_obj$parameters)
     )
 
   out
-
 }
 
 #' Handle covariate censoring
diff --git a/man/ss.Rd b/man/ss.Rd
new file mode 100644
index 0000000..aef365c
--- /dev/null
+++ b/man/ss.Rd
@@ -0,0 +1,18 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/misc.R
+\name{ss}
+\alias{ss}
+\title{Weighted sum-of-squares of residuals}
+\usage{
+ss(obs, pred, w = NULL)
+}
+\arguments{
+\item{obs}{observations vector}
+
+\item{pred}{predictions vector}
+
+\item{w}{weights}
+}
+\description{
+Weighted sum-of-squares of residuals
+}
diff --git a/tests/testthat/test-run_eval.R b/tests/testthat/test-run_eval.R
index 7721b17..4f50953 100644
--- a/tests/testthat/test-run_eval.R
+++ b/tests/testthat/test-run_eval.R
@@ -19,10 +19,24 @@ test_that("Basic run with vanco data + model works", {
   expect_equal(names(res), c("results", "mod_obj", "data", "sim", "stats_summ", "shrinkage", "bayesian_impact"))
   expect_s3_class(res$results, c("tbl_df", "tbl", "data.frame"))
   expect_s3_class(res$stats_summ, c("tbl_df", "tbl", "data.frame"))
+  expect_equal(
+    names(res$results),
+    c("id", "_iteration", "_grouper", "t", "dv", "pred", "map_ipred",
+      "ofv", "ss_w", "iter_ipred", "apriori", "CL", "V", "TH_CRCL",
+      "Q", "V2")
+  )
   expect_equal(
     round(res$results$CL[1:5], 3),
     c(2.99, 2.685, 2.462, 2.430, 2.439)
   )
+  expect_equal(
+    round(res$results$ofv[1:5], 3),
+    c(NA, 6.196, 11.846, 17.046, 21.563)
+  )
+  expect_equal(
+    round(res$results$ss_w[1:5], 3),
+    c(NA, 16.839, 11.241, 27.358, 30.944)
+  )
 
   # Using PKPDsim model library:
   res_2 <- run_eval(