Suggestions for VPC simulation

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,7 +22,7 @@ export(parse_psn_proseval_results)
 export(reldiff_psn_execute_results)
 export(reldiff_psn_proseval_results)
 export(run_eval)
-export(run_vpc)
+export(vpc_options)
 importFrom(PKPDsim,install_default_literature_model)
 importFrom(PKPDsim,new_ode_model)
 importFrom(rlang,.data)

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/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,6 +59,7 @@ run_eval <- function(
   weight_prior = 1,
   censor_covariates = TRUE,
   incremental = FALSE,
+  .vpc_options = vpc_options(),
   threads = 1,
   progress = TRUE,
   verbose = TRUE
@@ -104,41 +107,50 @@ run_eval <- function(
 
   ## 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
-    )
+  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
+  )
+  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
+  )
+  # 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)
+  }
 
   ## 5. Return results
   out

R/run_vpc.R

@@ -1,139 +1,21 @@
-#' Run iterative predictive analysis, looping over each individual's data
-#' to create observation and simulation datasets for creating visual predictive
-#' checks (VPCs).
+#' Core function for creating visual predictive checks (VPCs). Runs `n_samples`
+#' simulations for a single subject
 #'
 #' @inheritParams run_eval
-#' @param n_samples number of iterations to perform for VPC simulations.
-#'
-#' @returns A named list of data frames with `obs` and `sim` data.
-#'
-#' @examples
-#' \dontrun{
-#' library(mipdeval)
-#' library(vpc)
-#' vpc_vanc <- run_vpc(
-#'   model = "pkvancothomson",
-#'   data = nm_vanco,
-#'   n_samples = 100,
-#'   progress = FALSE
-#' )
-#' vpc(
-#'   sim = vpc_vanc$sim,
-#'   obs = vpc_vanc$obs
-#' )
-#' }
-#'
-#' @export
-run_vpc <- function(
-    model,
-    data,
-    ids = NULL,
-    n_samples = 250,
-    parameters = NULL,
-    fixed = NULL,
-    omega = NULL,
-    iov = NULL,
-    ruv = NULL,
-    dictionary = list(),
-    group = NULL,
-    progress = TRUE,
-    threads = 1,
-    verbose = TRUE
-) {
-
-  ## 0. Gather model information in an object
-  mod_obj <- parse_model(
-    model,
-    parameters = parameters,
-    ruv = ruv,
-    omega = omega,
-    fixed = fixed,
-    iov = iov
-  )
-
-  ## 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
-    )
-
-  ## Select covariates
-  covariates <- attr(mod_obj$model, "covariates")
-
-  ## 2. parse into separate, individual-level datasets and parse into
-  ##    format convenient for PKPDsim/PKPDmap, joined into a list object:
-  data_parsed <- parse_input_data(
-    data = input_data,
-    covariates = covariates,
-    ids = ids,
-    group = group
-  )
-
-  ## 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() { }
-  }
-
-  ## 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) {
-    future::plan(future::multisession, workers = threads)
-    res <- furrr::future_map(
-      .x = data_parsed,
-      .f = run_vpc_core,
-      mod_obj = mod_obj,
-      n_samples = n_samples,
-      progress_function = p
-    )
-  } else {
-    res <- purrr::map(
-      .x = data_parsed,
-      .f = run_vpc_core,
-      mod_obj = mod_obj,
-      n_samples = n_samples,
-      progress_function = p
-    )
-  }
-  sim_data <- dplyr::bind_rows(res) |>
-    dplyr::select(ID = .data$id, TIME = .data$t, DV = .data$y, ITER = .data$iter) |>
-    dplyr::arrange(.data$ITER, .data$ID, .data$TIME)
-
-  out <- list(
-    obs = input_data |>
-      dplyr::filter(.data$EVID == 0),
-    sim = sim_data
-  )
-
-  out
-
-}
-
-#' Core function for `run_vpc`. Runs `n_samples` simulations for
-#' a single subject
-#'
-#' @inheritParams run_vpc
 #' @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
+    n_samples = 250,
+    seed = 123
 ) {
-
   progress_function()
 
-  obs_data <- data$observations
-  comb <- data.frame()
-
   ## Simulate using PKPDsim
   bins <- NULL
   if(length(mod_obj$bins) > 0) {
@@ -149,13 +31,47 @@ run_vpc_core <- function(
     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!
-    )
+      iter = .data$id,
+      id = data$id # data, not .data!
+    ) |>
+    dplyr::select(ID = "id", TIME = "t", DV = "y", ITER = "iter")
+}
 
-  dat
+#' 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")
 }