Add accuracy metric to results stats summary

NAMESPACE

@@ -8,6 +8,7 @@ S3method(print,mipdeval_results)
 S3method(print,mipdeval_results_bayesian_impact)
 S3method(print,mipdeval_results_shrinkage)
 S3method(print,mipdeval_results_stats_summ)
+export(accuracy)
 export(add_grouping_column)
 export(calculate_bayesian_impact)
 export(calculate_shrinkage)
@@ -17,13 +18,19 @@ export(compare_psn_proseval_results)
 export(group_by_dose)
 export(group_by_time)
 export(install_default_literature_model)
+export(is_accurate)
+export(is_accurate_abs)
+export(is_accurate_rel)
 export(new_ode_model)
 export(parse_psn_proseval_results)
 export(reldiff_psn_execute_results)
 export(reldiff_psn_proseval_results)
 export(run_eval)
+export(stats_summ_options)
 export(vpc_options)
 importFrom(PKPDsim,install_default_literature_model)
 importFrom(PKPDsim,new_ode_model)
 importFrom(rlang,.data)
+importFrom(rlang,caller_arg)
+importFrom(rlang,caller_env)
 importFrom(utils,read.csv)

NEWS.md

@@ -1,5 +1,7 @@
 # mipdeval (development version)
 
+* Added `accuracy()` statistic for calculating the proportion of predicted drug concentrations that fall within specified absolute and relative error margins. This statistic will only be calculated in `run_eval()` when the absolute and relative error margins have been specified in the `.stats_summ_options` argument (see `stats_summ_options()`). (#30)
+
 # mipdeval 0.1.0
 
 * Initial version.

R/calculate_stats.R

@@ -3,13 +3,19 @@
 #' @param res output object (`mipdeval_results`) from `run_eval()`, or
 #' `data.frame` with raw results.
 #' @param rounding number of decimals to round to.
+#' @param acc_error_abs,acc_error_rel For calculating [accuracy()]: Positive number
+#'   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.
 #'
 #' @returns tibble
 #'
 #' @export
 calculate_stats <- function(
     res,
-    rounding = 3
+    rounding = 3,
+    acc_error_abs = NULL,
+    acc_error_rel = NULL
 ) {
   if(inherits(res, "mipdeval_results")) {
     res <- res$results
@@ -23,10 +29,67 @@ calculate_stats <- function(
       rmse = rmse(.data$dv, .data$value),
       nrmse = nrmse(.data$dv, .data$value),
       mpe = mpe(.data$dv, .data$value),
-      mape = mape(.data$dv, .data$value)
+      mape = mape(.data$dv, .data$value),
+      accuracy = dplyr::if_else(
+        !is.null(acc_error_abs) & !is.null(acc_error_rel),
+        accuracy(.data$dv, .data$value, acc_error_abs, acc_error_rel),
+        NA
+      )
     ) |>
     dplyr::mutate(dplyr::across(dplyr::everything(), round, rounding)) |>
     dplyr::as_tibble()
   class(out) <- c("mipdeval_results_stats_summ", class(out))
   out
 }
+
+#' Options for summary statistics
+#'
+#' @param ... These dots are reserved for future extensibility and must be empty.
+#' @inheritParams calculate_stats
+#'
+#' @returns A list.
+#' @export
+stats_summ_options <- function(
+    ...,
+    rounding = 3,
+    acc_error_abs = NULL,
+    acc_error_rel = NULL
+) {
+  rlang::check_dots_empty()
+  check_required_accuracy(acc_error_abs, acc_error_rel)
+  out <- list(
+    rounding = vctrs::vec_assert(
+      rounding, ptype = numeric(), size = 1L, arg = "rounding"
+    ),
+    acc_error_abs = vec_assert_or_null(
+      acc_error_abs, ptype = numeric(), size = 1L, arg = "acc_error_abs"
+    ),
+    acc_error_rel = vec_assert_or_null(
+      acc_error_rel, ptype = numeric(), size = 1L, arg = "acc_error_rel"
+    )
+  )
+  structure(out, class = "mipdeval_stats_summ_options")
+}
+
+check_required_accuracy <- function(
+  acc_error_abs = NULL,
+  acc_error_rel = NULL,
+  call = caller_env()
+) {
+  acc_error_abs_null <- is.null(acc_error_abs) & !is.null(acc_error_rel)
+  acc_error_rel_null <- !is.null(acc_error_abs) & is.null(acc_error_rel)
+  if (!acc_error_abs_null & !acc_error_rel_null) return()
+  cli::cli_abort(
+    message = c(
+      paste0(
+        "{.arg acc_error_abs} and {.arg acc_error_rel} must both be specified ",
+        "to calculate accuracy."
+      ),
+      "x" = dplyr::case_when(
+        acc_error_abs_null ~ "{.arg acc_error_abs} is NULL.",
+        acc_error_rel_null ~ "{.arg acc_error_rel} is NULL."
+      )
+    ),
+    call = call
+  )
+}

R/mipdeval-package.R

@@ -2,7 +2,7 @@
 "_PACKAGE"
 
 ## usethis namespace: start
-#' @importFrom rlang .data
+#' @importFrom rlang .data caller_arg caller_env
 #' @importFrom utils read.csv
 ## usethis namespace: end
 NULL

R/misc.R

@@ -26,6 +26,25 @@ is_timevarying <- function(.data, .cols) {
   )
 }
 
+#' Assert an argument has known prototype and/or size or is NULL
+#'
+#' @inheritParams vctrs::vec_assert
+#'
+#' @returns Either throws an error or returns `x`, invisibly.
+vec_assert_or_null <- function(
+    x,
+    ptype = NULL,
+    size = NULL,
+    arg = caller_arg(x),
+    call = caller_env()
+) {
+  if (!is.null(x)) {
+    vctrs::vec_assert(x = x, ptype = ptype, size = size, arg = arg, call = call)
+  } else {
+    NULL
+  }
+}
+
 #' Root-mean-squared error
 #'
 #' @param obs observations vector
@@ -68,6 +87,63 @@ mpe <- function (obs, pred) {
   sum((obs - pred)/obs)/length(obs)
 }
 
+#' Accuracy
+#'
+#' Accuracy provides a measure of clinical suitability, defined by whether model
+#' predicted drug concentrations fall within an absolute OR relative error
+#' margin of the measured concentrations.
+#'
+#' @param obs Observations vector.
+#' @param pred Predictions vector.
+#' @param error_abs,error_rel Positive number providing an absolute or relative
+#'   error margin. The cutoff is exclusive of the error margin. Defaults to `0`,
+#'   meaning no predictions fall within the error margin.
+#'
+#' @returns For `is_accurate()`, `is_accurate_abs()`, and `is_accurate_rel()`: A
+#'   logical vector indicating whether or not each predicted drug concentration
+#'   was considered accurate according to the specified absolute or relative
+#'   error margin(s).
+#'
+#'   For `accuracy()`: A single value between 0 and 1 indicating the proportion
+#'   of predicted drug concentrations that fell within the specified absolute
+#'   and relative error margins.
+#'
+#' @examples
+#' # Does the predicted drug concentration fall within 0.5 mg/L error margin?
+#' is_accurate_abs(6, 5, error_abs = 0.5)
+#'
+#' # Does the predicted drug concentration fall within 25% error margin?
+#' is_accurate_rel(6, 5, error_rel = 0.25)
+#'
+#' # Does the predicted drug concentration fall within 0.5 mg/L OR 25%?
+#' is_accurate(6, 5, error_abs = 0.5, error_rel = 0.25)
+#'
+#' # What proportion of predicted drug concentrations fell within 0.5 mg/L OR 25%?
+#' accuracy(rnorm(10, 6), rnorm(10, 5), error_abs = 0.5, error_rel = 0.25)
+#'
+#' @export
+accuracy <- function(obs, pred, error_abs = 0, error_rel = 0) {
+  mean(is_accurate(obs, pred, error_abs, error_rel))
+}
+
+#' @rdname accuracy
+#' @export
+is_accurate <- function(obs, pred, error_abs = 0, error_rel = 0) {
+  is_accurate_abs(obs, pred, error_abs) | is_accurate_rel(obs, pred, error_rel)
+}
+
+#' @rdname accuracy
+#' @export
+is_accurate_abs <- function(obs, pred, error_abs = 0) {
+  abs(pred - obs) < error_abs
+}
+
+#' @rdname accuracy
+#' @export
+is_accurate_rel <- function(obs, pred, error_rel = 0) {
+  (pred/obs > 1 - error_rel) & (pred/obs < 1 + error_rel)
+}
+
 #' Weighted sum-of-squares of residuals
 #'
 #' @inheritParams rmse
@@ -83,3 +159,4 @@ ss <- function(obs, pred, w = NULL) {
   if(sum(w) == 0) return(NA)
   sum(w * (obs - pred)^2)
 }
+

R/run_eval.R

@@ -28,6 +28,8 @@
 #'   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 `FALSE`.
+#' @param .stats_summ_options Options for summary statistics. This must be the
+#'   result from a call to [stats_summ_options()].
 #' @param .vpc_options Options for VPC simulations. This must be the result from
 #'   a call to [vpc_options()].
 #' @param threads number of threads to divide computations on. Default is 1,
@@ -59,11 +61,15 @@ run_eval <- function(
   weight_prior = 1,
   censor_covariates = TRUE,
   incremental = FALSE,
+  .stats_summ_options = stats_summ_options(),
   .vpc_options = vpc_options(),
   threads = 1,
   progress = TRUE,
   verbose = TRUE
 ) {
+  # Evaluate options at the start so any errors are triggered immediately.
+  .stats_summ_options <- .stats_summ_options
+  .vpc_options <- .vpc_options
 
   if(progress) { # configure progressbars
     progressr::handlers(global = TRUE)
@@ -169,7 +175,12 @@ run_eval <- function(
   # res is NULL when vpc_options(..., vpc_only = TRUE).
   if (!is.null(res)) {
     if(verbose) cli::cli_progress_step("Calculating forecasting statistics")
-    out$stats_summ <- calculate_stats(out)
+    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
+    )
     out$shrinkage <- calculate_shrinkage(out)
     out$bayesian_impact <- calculate_bayesian_impact(out)
     cli::cli_progress_done()