Add a function to create VPC data

DESCRIPTION

@@ -28,11 +28,14 @@ Imports:
     rlang,
     tibble,
     tidyr,
+    vctrs,
     withr
 Suggests: 
     knitr,
     rmarkdown,
     testthat (>= 3.2.0),
+    vdiffr,
+    vpc
 Remotes: 
     InsightRX/PKPDmap,
     InsightRX/PKPDsim 

NAMESPACE

@@ -3,6 +3,7 @@
 S3method(parse_model,PKPDsim)
 S3method(parse_model,character)
 S3method(parse_model,default)
+S3method(plot,mipdeval_results)
 S3method(print,mipdeval_results)
 S3method(print,mipdeval_results_bayesian_impact)
 S3method(print,mipdeval_results_shrinkage)
@@ -21,6 +22,7 @@ export(parse_psn_proseval_results)
 export(reldiff_psn_execute_results)
 export(reldiff_psn_proseval_results)
 export(run_eval)
+export(vpc_options)
 importFrom(PKPDsim,install_default_literature_model)
 importFrom(PKPDsim,new_ode_model)
 importFrom(rlang,.data)

R/check_input_data.R

@@ -6,7 +6,8 @@
 #' @inheritParams run_eval
 #'
 #' @returns a data.frame
-check_input_data <- function(data, dictionary) {
+check_input_data <- function(data, dictionary, verbose = TRUE) {
+  if(verbose) cli::cli_progress_step("Checking integrity of input data")
   if(!is.null(dictionary) && length(dictionary) > 0) {
     check_dictionary_columns(data, dictionary)
     check_valid_dictionary(dictionary)
@@ -21,9 +22,10 @@ check_input_data <- function(data, dictionary) {
         )
       )
   } else {
-    cli::cli_alert_info("No data `dictionary` provided, assuming common NONMEM column names.")
+    cli::cli_alert_info("No data `dictionary` provided, assuming common NONMEM column names in input dataset.")
   }
   check_required_cols(data)
+  if(verbose) cli::cli_progress_done()
   data
 }
 

R/parse_input_data.R

@@ -9,9 +9,12 @@ parse_input_data <- function(
     data,
     covariates,
     ids,
-    group = NULL
+    group = NULL,
+    verbose = TRUE
 ) {
 
+  if(verbose) cli::cli_progress_step("Parsing input data for analysis")
+
   ## Filter IDs, if needed
   if(!is.null(ids)) {
     cli::cli_alert_info("Running analysis on subset of data (n = {length(ids)})")
@@ -31,6 +34,8 @@ parse_input_data <- function(
     group = group
   )
 
+  if(verbose) cli::cli_progress_done()
+
   out
 }
 
@@ -52,6 +57,7 @@ parse_nm_data <- function(data, covariates, group = NULL) {
     dplyr::filter(.data$EVID == 0) |>
     dplyr::select("ID", "TIME", "DV", "CMT", "_grouper") |>
     dplyr::rename(id = "ID", t = "TIME", y = "DV", cmt = "CMT")
+  out$id <- obs$ID[1]
   out
 }
 

R/plot.R

@@ -0,0 +1,31 @@
+#' Plot method for a `run_eval()` object
+#'
+#' @description
+#' One plot `type` is currently available:
+#'
+#' - `type = "vpc"` plots visual predictive checks (VPC) with [vpc::vpc()].
+#'
+#' @param x An object.
+#' @param type Character string, indicating the type of plot. Options are `"vpc"`.
+#' @param ... Arguments passed to or from other methods.
+#'
+#' @returns A plot.
+#' @export
+plot.mipdeval_results <- function(x, type = "vpc", ...) {
+  type <- rlang::arg_match(type)
+  plot_fun <- switch (
+    type,
+    vpc = plot_vpc
+  )
+  plot_fun(x, ...)
+}
+
+plot_vpc <- function(res, ...) {
+  rlang::check_installed("vpc", reason = "for VPC plotting.")
+  rlang::check_dots_used()
+  vpc::vpc(
+    sim = res$sim,
+    obs = dplyr::filter(res$data, .data$EVID == 0),
+    ...
+  )
+}

R/read_input_data.R

@@ -3,7 +3,8 @@
 #' @inheritParams run_eval
 #'
 #' @returns data.frame
-read_input_data <- function(data) {
+read_input_data <- function(data, verbose = TRUE) {
+  if(verbose) cli::cli_progress_step("Reading input data")
   if(inherits(data, "character")) {
     if(!file.exists(data)) {
       cli::cli_abort("Sorry, filename supplied to `data` does not exist.")
@@ -14,5 +15,6 @@ read_input_data <- function(data) {
   } else {
     cli::cli_abort("Sorry, object supplied to `data` argument should be either a filename or a data.frame.")
   }
+  if(verbose) cli::cli_progress_done()
   input_data
 }

R/run.R

@@ -0,0 +1,17 @@
+# wrapper for purrr::map() that supports parallel processing via
+# furrr::future_map(), as well as optional skipping.
+run <- function(.x, .f, ..., .threads = 1, .skip = FALSE) {
+  if (isTRUE(.skip)) {
+    return()
+  } else if (.threads > 1) {
+    # TODO: consider using purrr::in_parallel() in the future when it's stable.
+    future::plan(future::multisession, workers = .threads)
+    res <- purrr::list_rbind(furrr::future_map(
+      .x = .x, .f = .f, ...
+    ))
+  } else {
+    res <- purrr::list_rbind(purrr::map(
+      .x = .x, .f = .f, ...
+    ))
+  }
+}

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 .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,
 #'   i.e. no parallel execution
 #' @param ruv residual error variability magnitude, specified as list.
@@ -57,11 +59,17 @@ run_eval <- function(
   weight_prior = 1,
   censor_covariates = TRUE,
   incremental = FALSE,
+  .vpc_options = vpc_options(),
   threads = 1,
   progress = TRUE,
   verbose = TRUE
 ) {
 
+  if(progress) { # configure progressbars
+    progressr::handlers(global = TRUE)
+    progressr::handlers("cli")
+  }
+
   ## 0. Gather model information in an object
   mod_obj <- parse_model(
     model,
@@ -73,10 +81,10 @@ run_eval <- function(
   )
 
   ## 1. read NONMEM data from file or data.frame. Do some simple checks
-  if(verbose) cli::cli_progress_step("Reading and parsing input data")
-  input_data <- read_input_data(data) |>
+  input_data <- read_input_data(data, verbose = verbose) |>
     check_input_data(
-      dictionary = dictionary
+      dictionary = dictionary,
+      verbose = verbose
     )
 
   ## Select covariates
@@ -88,57 +96,83 @@ run_eval <- function(
     data = input_data,
     covariates = covariates,
     ids = ids,
-    group = group
+    group = group,
+    verbose = verbose
   )
 
   ## Set up progress bars
   ## Need to loop this through the progressr package
   ## because furrr currently doesn't support progressbars with cli.
-  if(progress) {
-    progressr::handlers(global = TRUE)
-    progressr::handlers("cli")
-    p <- progressr::progressor(along = data_parsed)
-  } else {
-    p <- function() { }
-  }
+  p <- if(progress) { progressr::progressor(along = data_parsed) } else { \(x) {} }
 
   ## 3. run the core function on each individual-level dataset in the list
-  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)
-    res <- furrr::future_map(
-      .x = data_parsed,
-      .f = run_eval_core,
-      mod_obj = mod_obj,
-      censor_covariates = censor_covariates,
-      weight_prior = weight_prior,
-      progress_function = p
-    )
-  } else {
-    res <- purrr::map(
-      .x = data_parsed,
-      .f = run_eval_core,
-      mod_obj = mod_obj,
-      censor_covariates = censor_covariates,
-      weight_prior = weight_prior,
-      progress_function = p
+  if(verbose) {
+    cli::cli_alert_info("Running forecasts for subjects in dataset")
+    cli::cli_progress_step(
+      msg = NULL,
+      msg_done = "Forecasting analysis done."
     )
   }
+  res <- run(
+    .x = data_parsed,
+    .f = run_eval_core,
+    mod_obj = mod_obj,
+    censor_covariates = censor_covariates,
+    weight_prior = weight_prior,
+    progress_function = p,
+    .threads = threads,
+    .skip = .vpc_options$vpc_only
+  )
+  cli::cli_progress_done()
+
+  if(verbose) {
+    if(.vpc_options$skip) {
+      cli::cli_alert_info("Skipping simulations for VPC / NPDE")
+    } else {
+      cli::cli_alert_info("Running simulations for VPC / NPDE")
+      cli::cli_progress_step(
+        msg = NULL,
+        msg_done = "Simulations for VPC / NPDE done"
+      )
+    }
+  }
+  p <- if(progress) { progressr::progressor(along = data_parsed) } else { \(x) {} }
+  res_vpc <- run(
+    .x = data_parsed,
+    .f = run_vpc_core,
+    mod_obj = mod_obj,
+    n_samples = .vpc_options$n_samples,
+    seed = .vpc_options$seed,
+    progress_function = p,
+    .threads = threads,
+    .skip = .vpc_options$skip
+  )
+  cli::cli_progress_done()
+  # NULL is returned if VPC is skipped.
+  if (!is.null(res_vpc)) {
+    res_vpc <- dplyr::arrange(res_vpc, .data$ITER, .data$ID, .data$TIME)
+  }
 
   ## 4. Combine results and basic stats into return object
-  if(verbose) cli::cli_progress_step("Calculating forecasting statistics")
-  res_tbl <- dplyr::bind_rows(res) |>
-    tibble::as_tibble()
   out <- list(
-    results = res_tbl,
+    results = res,
     mod_obj = mod_obj,
-    data = data
+    data = input_data,
+    sim = res_vpc,
+    stats_summ = NULL,
+    shrinkage = NULL,
+    bayesian_impact = NULL
   )
   class(out) <- c("mipdeval_results", "list")
-  out$stats_summ <- calculate_stats(out)
-  out$shrinkage <- calculate_shrinkage(out)
-  out$bayesian_impact <- calculate_bayesian_impact(out)
+
+  # 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$shrinkage <- calculate_shrinkage(out)
+    out$bayesian_impact <- calculate_bayesian_impact(out)
+    cli::cli_progress_done()
+  }
 
   ## 5. Return results
   out

R/run_vpc.R

@@ -0,0 +1,77 @@
+#' Core function for creating visual predictive checks (VPCs). Runs `n_samples`
+#' simulations for a single subject
+#'
+#' @inheritParams run_eval
+#' @inheritParams run_eval_core
+#' @param n_samples Number of iterations to perform for VPC simulations.
+#' @param seed Seed for random number generation.
+#'
+#' @returns data.frame
+run_vpc_core <- function(
+    data,
+    mod_obj,
+    progress_function,
+    n_samples = 250,
+    seed = 123
+) {
+  progress_function()
+
+  ## Simulate using PKPDsim
+  bins <- NULL
+  if(length(mod_obj$bins) > 0) {
+    bins <- mod_obj$bins
+  }
+  dat <- PKPDsim::sim(
+    ode = mod_obj$model,
+    parameters = mod_obj$parameters,
+    omega = mod_obj$omega,
+    res_var = mod_obj$ruv,
+    t_obs = data$observations$t,
+    covariates = data$covariates,
+    regimen = data$regimen,
+    iov_bins = bins,
+    n_ind = n_samples,
+    seed = seed,
+    only_obs = TRUE,
+    verbose = FALSE
+  )
+
+  dat |>
+    dplyr::mutate(
+      iter = .data$id,
+      id = data$id # data, not .data!
+    ) |>
+    dplyr::select(ID = "id", TIME = "t", DV = "y", ITER = "iter")
+}
+
+#' Options for VPC simulations
+#'
+#' @param ... These dots are reserved for future extensibility and must be empty.
+#' @inheritParams run_vpc_core
+#' @param skip Logical. Skip VPC simulations?
+#' @param vpc_only Logical. Only return VPC simulations?
+#'
+#' @returns A list.
+#' @export
+vpc_options <- function(
+  ...,
+  n_samples = 250,
+  seed = 123,
+  skip = FALSE,
+  vpc_only = FALSE
+) {
+  rlang::check_dots_empty()
+  out <- list(
+    n_samples = vctrs::vec_assert(
+      n_samples, ptype = numeric(), size = 1L, arg = "n_samples"
+    ),
+    seed = seed,
+    skip = vctrs::vec_assert(
+      skip, ptype = logical(), size = 1L, arg = "skip"
+    ),
+    vpc_only = vctrs::vec_assert(
+      vpc_only, ptype = logical(), size = 1L, arg = "vpc_only"
+    )
+  )
+  structure(out, class = "mipdeval_vpc_options")
+}