From 4d6b44c32048c141ef7c76113e42bc41808fd3e8 Mon Sep 17 00:00:00 2001 From: Michael McCarthy <51542091+mccarthy-m-g@users.noreply.github.com> Date: Thu, 4 Jun 2026 09:32:47 -0700 Subject: [PATCH 1/5] add etas to results table --- R/run_eval.R | 60 +++++++++++++++++++++++++++++++- R/run_eval_core.R | 26 +++++++++++--- man/run_eval.Rd | 63 +++++++++++++++++++++++++++++++++- tests/testthat/test-run_eval.R | 3 +- 4 files changed, 144 insertions(+), 8 deletions(-) diff --git a/R/run_eval.R b/R/run_eval.R index c04c3cf..24d85eb 100644 --- a/R/run_eval.R +++ b/R/run_eval.R @@ -45,7 +45,65 @@ #' debugging the package it is useful to have it off, since progress bar #' handlers obscure R output. #' -#' @returns A named list of data frames. +#' @details +#' `run_eval()` evaluates predictive performance the way a model would be used +#' for model-informed precision dosing (MIPD): for each subject it walks through +#' the observations in time order and, at each step, refits the individual using +#' only the data available *up to that point* to predict the next observation(s). +#' This produces an a priori (population) prediction followed by progressively +#' more informed a posteriori (forecasting) predictions, mirroring the iterative +#' flow of PsN's `proseval` tool. +#' +#' Because of this design the `results` tibble reports quantities at three levels +#' of individualization, which are easy to confuse. The predictions: +#' +#' - `pred`: the **population** prediction (typical individual, using no +#' individual-level data)---the a priori prediction. +#' - `iter_ipred`: the **iterative ("forecasting")** individual prediction, using +#' only the data available up to that point---the a posteriori prediction. On +#' the a priori rows it falls back to `pred`. +#' - `map_ipred`: the **full-data MAP** individual prediction, from a single +#' retrospective fit on all of a subject's data at once. +#' +#' The random-effect (eta) estimates follow the same structure: +#' +#' - The `eta` columns are the iterative ("forecasting") empirical Bayes +#' estimates that pair with `iter_ipred`; they evolve down the rows and are 0 +#' on the a priori (population) rows. +#' - The `map_eta` columns are the full-data MAP empirical Bayes estimates +#' that pair with `map_ipred`. They are constant per subject and appear on +#' every row, including the a priori row. These are the equivalent of the etas +#' reported in a NONMEM output table, and are what you want for an empirical +#' eta-distribution plot (one value per subject). +#' +#' @returns An `mipdeval_results` object, which is a named list with the following +#' elements: +#' +#' - `results`: A tibble with one row per (iterative) prediction, holding the +#' identifiers (`id`, `_iteration`, `_grouper`, `t`), the observation (`dv`), +#' population predictions and residuals (`pred`, `res`, `wres`, `cwres`), +#' individual predictions and residuals (`iter_ipred`, `map_ipred`, `ires`, +#' `iwres`), the objective function value (`ofv`), the weighted sum-of-squares +#' (`ss_w`), an `apriori` flag, one column per model parameter, and two +#' families of eta (random-effect) columns: the iterative `eta` and the +#' full-data `map_eta` estimates (see Details). +#' - `mod_obj`: The parsed model object (see [parse_model()]): a named list of +#' model information, including `model`, `parameters`, `omega`, `ruv`, +#' `fixed`, `bins`, and `kappa`. +#' - `data`: The input data after reading and validation (see +#' [read_input_data()] and [check_input_data()]), as a data frame of the +#' NONMEM-style records used in the analysis. +#' - `sim`: Simulated data used for the visual predictive check (VPC) and NPDE, +#' or `NULL` when simulations are skipped (`vpc_options(skip = TRUE)`). +#' - `stats_summ`: A tibble summarising forecasting performance statistics (see +#' [calculate_stats()]). +#' - `shrinkage`: A tibble of eta-shrinkage per iteration (see +#' [calculate_shrinkage()]). +#' - `bayesian_impact`: A tibble of Bayesian-impact values based on RMSE and +#' MAPE (see [calculate_bayesian_impact()]). +#' +#' `stats_summ`, `shrinkage`, and `bayesian_impact` are `NULL` when no +#' predictions are produced (e.g. `vpc_options(vpc_only = TRUE)`). #' #' @export run_eval <- function( diff --git a/R/run_eval_core.R b/R/run_eval_core.R index 7f105a1..965367d 100644 --- a/R/run_eval_core.R +++ b/R/run_eval_core.R @@ -23,6 +23,7 @@ run_eval_core <- function( obs_data <- data$observations comb <- data.frame() fit_pars <- data.frame() + eta_names <- character(0) # captured from first successful fit below iterations <- unique(obs_data[["_grouper"]]) for(i in seq_along(iterations)) { @@ -126,8 +127,9 @@ run_eval_core <- function( `_iteration` = iterations[i], `_grouper` = obs_data$`_grouper` ) - ## Add parameter estimates - fit_pars <- dplyr::mutate(as.data.frame(fit$parameters), id = obs_data$id[1]) + ## Add parameter estimates and etas: + eta_names <- names(fit$fit$par) + fit_pars <- dplyr::mutate(as.data.frame(c(fit$parameters, fit$fit$par)), id = obs_data$id[1]) } comb <- dplyr::bind_rows( @@ -175,6 +177,18 @@ run_eval_core <- function( map_pred_data <- pred_data } + ## Full-data MAP etas (empirical Bayes estimates), constant per subject and + ## named to parallel `map_ipred`. These differ from the iterative `eta_names` + ## columns, which only use the data available up to each forecast. NA if the + ## MAP fit failed. + map_eta_names <- paste0("map_", eta_names) + map_etas <- if(!is.null(fit_map$fit$par)) { + as.list(fit_map$fit$par[eta_names]) + } else { + as.list(rep(NA_real_, length(eta_names))) + } + names(map_etas) <- map_eta_names + ## pre-pend population predictions for the first observation # TODO: Refactor this logic into a function or functions, e.g., the first # argument to bind_rows() could be refactored into `get_apriori_data()`. @@ -185,7 +199,8 @@ run_eval_core <- function( `_iteration` = 0, ipred = .data$pred, ofv = NA, - ss_w = NA + ss_w = NA, + dplyr::across(dplyr::all_of(eta_names), ~ 0) # population etas are 0 ) |> # set to population parameters, not individual estimates dplyr::select(-!!names(mod_obj$parameters)) |> dplyr::left_join( @@ -198,12 +213,13 @@ run_eval_core <- function( dplyr::mutate( iter_ipred = .data$ipred, map_ipred = map_pred_data$ipred, # ipred from full retrospective MAP - apriori = (.data$`_iteration` == 0) + apriori = (.data$`_iteration` == 0), + !!!map_etas # full-data MAP etas, constant per subject ) |> dplyr::select( "id", "_iteration", "_grouper", "t", "dv", "pred", "res", "wres", "cwres", "map_ipred", "ofv", "ss_w", "iter_ipred", "ires", "iwres", "apriori", - !!names(mod_obj$parameters) + !!names(mod_obj$parameters), !!eta_names, !!map_eta_names ) out diff --git a/man/run_eval.Rd b/man/run_eval.Rd index 9c0c94c..8518afc 100644 --- a/man/run_eval.Rd +++ b/man/run_eval.Rd @@ -102,8 +102,69 @@ handlers obscure R output.} \item{verbose}{show more output} } \value{ -A named list of data frames. +An \code{mipdeval_results} object, which is a named list with the following +elements: +\itemize{ +\item \code{results}: A tibble with one row per (iterative) prediction, holding the +identifiers (\code{id}, \verb{_iteration}, \verb{_grouper}, \code{t}), the observation (\code{dv}), +population predictions and residuals (\code{pred}, \code{res}, \code{wres}, \code{cwres}), +individual predictions and residuals (\code{iter_ipred}, \code{map_ipred}, \code{ires}, +\code{iwres}), the objective function value (\code{ofv}), the weighted sum-of-squares +(\code{ss_w}), an \code{apriori} flag, one column per model parameter, and two +families of eta (random-effect) columns: the iterative \verb{eta} and the +full-data \verb{map_eta} estimates (see Details). +\item \code{mod_obj}: The parsed model object (see \code{\link[=parse_model]{parse_model()}}): a named list of +model information, including \code{model}, \code{parameters}, \code{omega}, \code{ruv}, +\code{fixed}, \code{bins}, and \code{kappa}. +\item \code{data}: The input data after reading and validation (see +\code{\link[=read_input_data]{read_input_data()}} and \code{\link[=check_input_data]{check_input_data()}}), as a data frame of the +NONMEM-style records used in the analysis. +\item \code{sim}: Simulated data used for the visual predictive check (VPC) and NPDE, +or \code{NULL} when simulations are skipped (\code{vpc_options(skip = TRUE)}). +\item \code{stats_summ}: A tibble summarising forecasting performance statistics (see +\code{\link[=calculate_stats]{calculate_stats()}}). +\item \code{shrinkage}: A tibble of eta-shrinkage per iteration (see +\code{\link[=calculate_shrinkage]{calculate_shrinkage()}}). +\item \code{bayesian_impact}: A tibble of Bayesian-impact values based on RMSE and +MAPE (see \code{\link[=calculate_bayesian_impact]{calculate_bayesian_impact()}}). +} + +\code{stats_summ}, \code{shrinkage}, and \code{bayesian_impact} are \code{NULL} when no +predictions are produced (e.g. \code{vpc_options(vpc_only = TRUE)}). } \description{ Run iterative predictive analysis, looping over each individual's data } +\details{ +\code{run_eval()} evaluates predictive performance the way a model would be used +for model-informed precision dosing (MIPD): for each subject it walks through +the observations in time order and, at each step, refits the individual using +only the data available \emph{up to that point} to predict the next observation(s). +This produces an a priori (population) prediction followed by progressively +more informed a posteriori (forecasting) predictions, mirroring the iterative +flow of PsN's \code{proseval} tool. + +Because of this design the \code{results} tibble reports quantities at three levels +of individualization, which are easy to confuse. The predictions: +\itemize{ +\item \code{pred}: the \strong{population} prediction (typical individual, using no +individual-level data)---the a priori prediction. +\item \code{iter_ipred}: the \strong{iterative ("forecasting")} individual prediction, using +only the data available up to that point---the a posteriori prediction. On +the a priori rows it falls back to \code{pred}. +\item \code{map_ipred}: the \strong{full-data MAP} individual prediction, from a single +retrospective fit on all of a subject's data at once. +} + +The random-effect (eta) estimates follow the same structure: +\itemize{ +\item The \verb{eta} columns are the iterative ("forecasting") empirical Bayes +estimates that pair with \code{iter_ipred}; they evolve down the rows and are 0 +on the a priori (population) rows. +\item The \verb{map_eta} columns are the full-data MAP empirical Bayes estimates +that pair with \code{map_ipred}. They are constant per subject and appear on +every row, including the a priori row. These are the equivalent of the etas +reported in a NONMEM output table, and are what you want for an empirical +eta-distribution plot (one value per subject). +} +} diff --git a/tests/testthat/test-run_eval.R b/tests/testthat/test-run_eval.R index f45f595..02a4f57 100644 --- a/tests/testthat/test-run_eval.R +++ b/tests/testthat/test-run_eval.R @@ -24,7 +24,8 @@ test_that("Basic run with vanco data + model works", { names(res$results), c("id", "_iteration", "_grouper", "t", "dv", "pred", "res", "wres", "cwres", "map_ipred", "ofv", "ss_w", "iter_ipred", "ires", "iwres", "apriori", "CL", - "V", "TH_CRCL", "Q", "V2") + "V", "TH_CRCL", "Q", "V2", "TDM_INIT", "eta01", "eta02", "eta03", "eta04", + "map_eta01", "map_eta02", "map_eta03", "map_eta04") ) expect_equal( round(res$results$CL[1:5], 3), From 2a85686fa7dd6564632528bf0d9e110c43f2bc29 Mon Sep 17 00:00:00 2001 From: Michael McCarthy <51542091+mccarthy-m-g@users.noreply.github.com> Date: Thu, 4 Jun 2026 15:54:27 -0700 Subject: [PATCH 2/5] remove "TDM_INIT" from expected output --- tests/testthat/test-run_eval.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-run_eval.R b/tests/testthat/test-run_eval.R index 02a4f57..a0c5cad 100644 --- a/tests/testthat/test-run_eval.R +++ b/tests/testthat/test-run_eval.R @@ -24,7 +24,7 @@ test_that("Basic run with vanco data + model works", { names(res$results), c("id", "_iteration", "_grouper", "t", "dv", "pred", "res", "wres", "cwres", "map_ipred", "ofv", "ss_w", "iter_ipred", "ires", "iwres", "apriori", "CL", - "V", "TH_CRCL", "Q", "V2", "TDM_INIT", "eta01", "eta02", "eta03", "eta04", + "V", "TH_CRCL", "Q", "V2", "eta01", "eta02", "eta03", "eta04", "map_eta01", "map_eta02", "map_eta03", "map_eta04") ) expect_equal( From 41eb886615018736ebad157d3075b00fa80c3705 Mon Sep 17 00:00:00 2001 From: Michael McCarthy <51542091+mccarthy-m-g@users.noreply.github.com> Date: Thu, 4 Jun 2026 15:58:43 -0700 Subject: [PATCH 3/5] export error metrics --- NAMESPACE | 4 ++++ R/misc.R | 5 ++++- 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/NAMESPACE b/NAMESPACE index 683b32d..56d7d2e 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -22,10 +22,14 @@ export(install_default_literature_model) export(is_accurate) export(is_accurate_abs) export(is_accurate_rel) +export(mape) +export(mpe) export(new_ode_model) +export(nrmse) export(parse_psn_proseval_results) export(reldiff_psn_execute_results) export(reldiff_psn_proseval_results) +export(rmse) export(run_eval) export(stats_summ_options) export(vpc_options) diff --git a/R/misc.R b/R/misc.R index 40b13ec..429e51c 100644 --- a/R/misc.R +++ b/R/misc.R @@ -51,6 +51,7 @@ vec_assert_or_null <- function( #' @param pred predictions vector #' #' @returns A numeric vector +#' @export rmse <- function (obs, pred) { res_sq <- (pred - obs)^2 sqrt(mean(res_sq, na.rm = TRUE)) @@ -62,7 +63,7 @@ rmse <- function (obs, pred) { #' @param pred predictions vector #' #' @returns A numeric vector -#' +#' @export nrmse <- function (obs, pred) { res_sq <- (pred - obs)^2 rmse <- sqrt(mean(res_sq, na.rm = T)) @@ -74,6 +75,7 @@ nrmse <- function (obs, pred) { #' @inheritParams rmse #' #' @returns A numeric vector +#' @export mape <- function (obs, pred) { sum(abs((obs - pred))/obs)/length(obs) } @@ -83,6 +85,7 @@ mape <- function (obs, pred) { #' @inheritParams rmse #' #' @returns A numeric vector +#' @export mpe <- function (obs, pred) { sum((obs - pred)/obs)/length(obs) } From c366ba51f22840701aa4fad99a61573525d481aa Mon Sep 17 00:00:00 2001 From: Michael McCarthy <51542091+mccarthy-m-g@users.noreply.github.com> Date: Fri, 5 Jun 2026 12:56:55 -0700 Subject: [PATCH 4/5] fail gracefully when all fits fail --- R/run_eval_core.R | 6 ++++- tests/testthat/test-run_eval.R | 41 ++++++++++++++++++++++++++++++++++ 2 files changed, 46 insertions(+), 1 deletion(-) diff --git a/R/run_eval_core.R b/R/run_eval_core.R index 5741d33..48550ea 100644 --- a/R/run_eval_core.R +++ b/R/run_eval_core.R @@ -96,6 +96,9 @@ run_eval_core <- function( ) par_dummy <- as.data.frame(mod_upd$parameters) par_dummy[, 1:ncol(par_dummy)] <- NA + ## Also emit NA eta columns if a previous fit told us their names, so a + ## failed fit has the same shape as a successful one. + par_dummy[, eta_names] <- NA_real_ fit_pars <- dplyr::mutate(as.data.frame(par_dummy), id = obs_data$id[1]) } else { ## Data frame with predictive data @@ -169,12 +172,13 @@ run_eval_core <- function( ## named to parallel `map_ipred`. These differ from the iterative `eta_names` ## columns, which only use the data available up to each forecast. NA if the ## MAP fit failed. - map_eta_names <- paste0("map_", eta_names) map_etas <- if(!is.null(fit_map$fit$par)) { as.list(fit_map$fit$par[eta_names]) } else { as.list(rep(NA_real_, length(eta_names))) } + ## guard against the no-successful-fit case where `eta_names` is empty. + map_eta_names <- if(length(eta_names) > 0) paste0("map_", eta_names) else character(0) names(map_etas) <- map_eta_names ## pre-pend population predictions for the first observation diff --git a/tests/testthat/test-run_eval.R b/tests/testthat/test-run_eval.R index a0c5cad..b42d63e 100644 --- a/tests/testthat/test-run_eval.R +++ b/tests/testthat/test-run_eval.R @@ -88,6 +88,47 @@ test_that("Run also works when `model` argument just references the package", { # TODO: test outputs }) +test_that("run_eval() returns NA results with a warning when all fits fail", { + local_mipdeval_options(rlib_warning_verbosity = "default") + + # Force every MAP Bayesian fit to fail. get_map_estimates returns (rather than + # throws) an error object on failure, so we mock that behaviour. + local_mocked_bindings( + get_map_estimates = function(...) simpleError("forced failure for test"), + .package = "PKPDmap" + ) + + mod_obj <- parse_model("pkvancothomson") + expect_warning( + res <- run_eval( + model = mod_obj$model, + data = nm_vanco, + parameters = mod_obj$parameters, + omega = mod_obj$omega, + ruv = mod_obj$ruv, + fixed = mod_obj$fixed, + censor_covariates = FALSE, + ids = c(1, 2), + .vpc_options = vpc_options(skip = TRUE), + progress = FALSE, + verbose = FALSE + ), + regexp = "Errors were encountered in 10 out of 10 evaluated predictions" + ) + + ## Should still return a usable object rather than erroring out: + expect_equal( + names(res), + c("results", "mod_obj", "data", "sim", "stats_summ", "shrinkage", "bayesian_impact") + ) + expect_s3_class(res, "mipdeval_results") + + ## Predictions are all NA, but the result structure is intact: + expect_all_true(is.na(res$results$pred)) + expect_all_true(is.na(res$results$iter_ipred)) + expect_all_true(is.na(res$results$map_ipred)) +}) + test_that("Flattening of prior results in different predictions", { local_mipdeval_options() res <- run_eval( From efd88e45a0377caac6efe8ef156c25a7492a325f Mon Sep 17 00:00:00 2001 From: Michael McCarthy <51542091+mccarthy-m-g@users.noreply.github.com> Date: Fri, 5 Jun 2026 13:39:26 -0700 Subject: [PATCH 5/5] refactor failed fits warning to occur in `run_eval()` instead of `calculate_stats()` --- CLAUDE.md | 2 ++ R/calculate_stats.R | 17 +++++---------- R/check_failed_fits.R | 25 +++++++++++++++++++++ R/run_eval.R | 4 +++- man/calculate_stats.Rd | 11 +++++++++- man/check_failed_fits.Rd | 20 +++++++++++++++++ tests/testthat/test-calculate_stats.R | 14 ++++++++++++ tests/testthat/test-check_failed_fits.R | 29 +++++++++++++++++++++++++ 8 files changed, 109 insertions(+), 13 deletions(-) create mode 100644 R/check_failed_fits.R create mode 100644 man/check_failed_fits.Rd create mode 100644 tests/testthat/test-check_failed_fits.R diff --git a/CLAUDE.md b/CLAUDE.md index c23af75..19b1820 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -46,6 +46,8 @@ The main entry point is `run_eval()` (`R/run_eval.R`), which orchestrates: 4. **Statistics** (`calculate_stats.R`, `calculate_shrinkage.R`, `calculate_bayesian_impact.R`): Computes RMSE, NRMSE, MAPE, MPE, accuracy, shrinkage, and Bayesian impact metrics. +Failed fits (which surface as `NA` predictions) are detected and reported by `check_failed_fits.R`. `run_eval()` calls it once after the per-subject loop to emit a single warning; `calculate_stats()` also calls it (controllable via its `warn` argument, which `run_eval()` sets to `FALSE` to avoid a duplicate warning). + ### Output Structure `run_eval()` returns a list with class `"mipdeval_results"`: diff --git a/R/calculate_stats.R b/R/calculate_stats.R index 43f199e..d94dfba 100644 --- a/R/calculate_stats.R +++ b/R/calculate_stats.R @@ -7,6 +7,8 @@ #' providing an absolute or relative error margin. The cutoff is exclusive of #' the error margin. When `NULL` (the default), accuracy will not be #' calculated and will return `NA` instead. +#' @param warn should a warning be emitted when failed fits (NA predictions) are +#' detected? #' #' @returns tibble #' @@ -15,21 +17,14 @@ calculate_stats <- function( .res, rounding = 3, acc_error_abs = NULL, - acc_error_rel = NULL + acc_error_rel = NULL, + warn = TRUE ) { if(inherits(.res, "mipdeval_results")) { .res <- .res$results } - ## Check for errors during fits / predictions - errors <- dplyr::filter( - .res, - is.na(.data$pred) | - (is.na(.data$map_ipred) & !.data$apriori) | - is.na(.data$iter_ipred) - ) - if(nrow(errors) > 0) { - cli::cli_warn("Errors were encountered in {nrow(errors)} out of {nrow(.res)} evaluated predictions. The problems occurred in patient(s) {unique(errors$id)}.") - } + ## Warn about any failed fits / predictions (NA), unless the caller opts out. + if (isTRUE(warn)) check_failed_fits(.res) out <- .res |> tidyr::pivot_longer( cols = c("pred", "map_ipred", "iter_ipred"), names_to = "type" diff --git a/R/check_failed_fits.R b/R/check_failed_fits.R new file mode 100644 index 0000000..43e14ad --- /dev/null +++ b/R/check_failed_fits.R @@ -0,0 +1,25 @@ +#' Check for failed fits / predictions and (optionally) warn +#' +#' Detects predictions that came back as `NA`, which indicates the underlying +#' MAP Bayesian fit failed, and---when `warn = TRUE`---emits a warning +#' summarising how many predictions failed and in which subjects. +#' +#' @param .res output object (`mipdeval_results`) from [run_eval()], or a +#' `data.frame` with raw results. +#' +#' @returns invisibly, a `data.frame` of the rows with failed predictions. +check_failed_fits <- function(.res) { + if(inherits(.res, "mipdeval_results")) { + .res <- .res$results + } + errors <- dplyr::filter( + .res, + is.na(.data$pred) | + (is.na(.data$map_ipred) & !.data$apriori) | + is.na(.data$iter_ipred) + ) + if(nrow(errors) > 0) { + cli::cli_warn("Errors were encountered in {nrow(errors)} out of {nrow(.res)} evaluated predictions. The problems occurred in patient(s) {unique(errors$id)}.") + } + invisible(errors) +} diff --git a/R/run_eval.R b/R/run_eval.R index 24d85eb..b4bfc1f 100644 --- a/R/run_eval.R +++ b/R/run_eval.R @@ -236,12 +236,14 @@ run_eval <- function( # res is NULL when vpc_options(..., vpc_only = TRUE). if (!is.null(res)) { + check_failed_fits(res) if(verbose) cli::cli_progress_step("Calculating forecasting statistics") out$stats_summ <- calculate_stats( out, rounding = .stats_summ_options$rounding, acc_error_abs = .stats_summ_options$acc_error_abs, - acc_error_rel = .stats_summ_options$acc_error_rel + acc_error_rel = .stats_summ_options$acc_error_rel, + warn = FALSE # Avoid a duplicate warning from check_failed_fits() ) out$shrinkage <- calculate_shrinkage(out) out$bayesian_impact <- calculate_bayesian_impact(out) diff --git a/man/calculate_stats.Rd b/man/calculate_stats.Rd index cc08c98..77e2eb7 100644 --- a/man/calculate_stats.Rd +++ b/man/calculate_stats.Rd @@ -4,7 +4,13 @@ \alias{calculate_stats} \title{Calculate basic statistics, like RMSE, MPE, MAPE for forecasted data} \usage{ -calculate_stats(.res, rounding = 3, acc_error_abs = NULL, acc_error_rel = NULL) +calculate_stats( + .res, + rounding = 3, + acc_error_abs = NULL, + acc_error_rel = NULL, + warn = TRUE +) } \arguments{ \item{.res}{output object (\code{mipdeval_results}) from \code{run_eval()}, or @@ -16,6 +22,9 @@ calculate_stats(.res, rounding = 3, acc_error_abs = NULL, acc_error_rel = NULL) providing an absolute or relative error margin. The cutoff is exclusive of the error margin. When \code{NULL} (the default), accuracy will not be calculated and will return \code{NA} instead.} + +\item{warn}{should a warning be emitted when failed fits (NA predictions) are +detected?} } \value{ tibble diff --git a/man/check_failed_fits.Rd b/man/check_failed_fits.Rd new file mode 100644 index 0000000..f25b322 --- /dev/null +++ b/man/check_failed_fits.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/check_failed_fits.R +\name{check_failed_fits} +\alias{check_failed_fits} +\title{Check for failed fits / predictions and (optionally) warn} +\usage{ +check_failed_fits(.res) +} +\arguments{ +\item{.res}{output object (\code{mipdeval_results}) from \code{\link[=run_eval]{run_eval()}}, or a +\code{data.frame} with raw results.} +} +\value{ +invisibly, a \code{data.frame} of the rows with failed predictions. +} +\description{ +Detects predictions that came back as \code{NA}, which indicates the underlying +MAP Bayesian fit failed, and---when \code{warn = TRUE}---emits a warning +summarising how many predictions failed and in which subjects. +} diff --git a/tests/testthat/test-calculate_stats.R b/tests/testthat/test-calculate_stats.R index ac3ef1f..fcc0c4e 100644 --- a/tests/testthat/test-calculate_stats.R +++ b/tests/testthat/test-calculate_stats.R @@ -52,6 +52,20 @@ test_that("calculate_stats() warns when predictions contain errors (NAs)", { expect_all_false(is.na(out$accuracy)) }) +test_that("calculate_stats(warn = FALSE) does not warn about failed fits", { + res <- data.frame( + id = c(1, 1, 2, 2, 3, 3), + apriori = c(FALSE, TRUE, FALSE, TRUE, FALSE, TRUE), + dv = c(10, 12, NA, NA, 8, 14), + pred = c(9, 13, NA, NA, 7, 15), + map_ipred = c(10, 11, NA, NA, 9, 13), + iter_ipred = c(11, 12, NA, NA, 8, 14) + ) + expect_no_warning( + calculate_stats(res, acc_error_abs = 0.2, acc_error_rel = 0.05, warn = FALSE) + ) +}) + # stats_summ_options() -------------------------------------------------------- test_that("stats_summ_options() works", { actual <- stats_summ_options( diff --git a/tests/testthat/test-check_failed_fits.R b/tests/testthat/test-check_failed_fits.R new file mode 100644 index 0000000..dde05be --- /dev/null +++ b/tests/testthat/test-check_failed_fits.R @@ -0,0 +1,29 @@ +test_that("check_failed_fits() warns and returns failed rows", { + res <- data.frame( + id = c(1, 1, 2, 2, 3, 3), + apriori = c(FALSE, TRUE, FALSE, TRUE, FALSE, TRUE), + dv = c(10, 12, NA, NA, 8, 14), + pred = c(9, 13, NA, NA, 7, 15), + map_ipred = c(10, 11, NA, NA, 9, 13), + iter_ipred = c(11, 12, NA, NA, 8, 14) + ) + expect_warning( + errors <- check_failed_fits(res), + "Errors were encountered in 2 out of 6 evaluated predictions" + ) + expect_equal(nrow(errors), 2) + expect_equal(unique(errors$id), 2) +}) + +test_that("check_failed_fits() does not warn when there are no failures", { + res <- data.frame( + id = c(1, 1, 2, 2), + apriori = c(FALSE, TRUE, FALSE, TRUE), + dv = c(10, 12, 8, 14), + pred = c(9, 13, 7, 15), + map_ipred = c(10, 11, 9, 13), + iter_ipred = c(11, 12, 8, 14) + ) + expect_no_warning(errors <- check_failed_fits(res)) + expect_equal(nrow(errors), 0) +})