Fixes for use of models with IOV + grouping of samples + example busulfan (#5)

* misc updates + skeleton for first versions main functions

* updates

* add initial version + 1st test

* updates

* add comp with proseval

* doc updates + minor fixes

* remove windows CI for now

* activate threads option

* allow specifying model info only with library name

* fix issue

* add ids argument; add test for prior weight

* updates

* fix iov

* updates for grouping

* update parsing of proseval when data is grouped

* update busulfan dataset

* update grouping functions

* add test grouping

* misc fixes

* misc fixes

* misc fixes

* fixes from copilot review

* make sure dataset is reproducible

* Grouping suggestions (#6)

* make sure `nm_vanco` is an ungrouped tibble

* make sure `nm_busulfan` is an ungrouped tibble

* use tidyselect syntax, not data-masking

* tidy grouping functions + add tests

* rebuild README

* refactor `check_input_data()`, add tests

* tidy documentation

* add TODO comments

* tidy tests

* tidy code

* clarify abbrev

* bump version

---------

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

Author: roninsightrx

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: mipdeval
 Title: Run evaluation of iterative predictive performance based
   on historical datasets
-Version: 0.0.0.9000
+Version: 0.1.0
 Authors@R: c(
     person("Ron", "Keizer", , "ron@insight-rx.com", role = c("aut", "cre")),
     person("Michael", "McCarthy", , "michael.mccarthy@insight-rx.com", role = c("aut")),

NAMESPACE

@@ -3,15 +3,17 @@
 S3method(parse_model,PKPDsim)
 S3method(parse_model,character)
 S3method(parse_model,default)
+export(add_grouping_column)
 export(compare_psn_execute_results)
 export(compare_psn_proseval_results)
+export(group_by_dose)
+export(group_by_time)
 export(install_default_literature_model)
 export(new_ode_model)
 export(parse_psn_proseval_results)
 export(reldiff_psn_execute_results)
 export(reldiff_psn_proseval_results)
 export(run_eval)
-export(run_eval_core)
 importFrom(PKPDsim,install_default_literature_model)
 importFrom(PKPDsim,new_ode_model)
 importFrom(rlang,.data)

R/add_grouping_column.R

@@ -0,0 +1,106 @@
+#' Add grouping column using a function
+#'
+#' @inheritParams run_eval
+#'
+#' @param fun Function used to define groups. Either (1) a quoted or character
+#'   name referencing a function (e.g., [group_by_time()], [group_by_dose()], or a
+#'   custom function), or (2) an anonymous function; see examples.
+#' @param label column name of grouper column in dataset
+#' @param ... Additional arguments passed onto `fun`. See the respective
+#'   grouping functions for more details. E.g. for `group_by_time`, need a
+#'   `bins` argument.
+#'
+#' @details
+#'
+#' All functions supplied to `add_grouping_column()` must take at least the
+#' argument `data`, and return a numeric vector (or a vector that can be cast to
+#' numeric) of the same length as the number of rows in `data`, which will be
+#' used as the grouping column.
+#'
+#' @returns data.frame
+#' @seealso [group_by_time()], [group_by_dose()]
+#' @examples
+#' # group_by_dose:
+#' add_grouping_column(nm_busulfan, fun = group_by_dose, label = "group")
+#'
+#' # group_by_time:
+#' add_grouping_column(nm_busulfan, fun = "group_by_time", label = "group")
+#'
+#' # Anonymous function:
+#' add_grouping_column(
+#'   nm_busulfan,
+#'   fun = function(data) {
+#'     as.numeric(cut(
+#'       data$TIME, breaks = c(0, 24, 48, 72, 96, Inf), include.lowest = TRUE
+#'     ))
+#'   },
+#'   label = "group"
+#' )
+#'
+#' @export
+add_grouping_column <- function(
+  data,
+  fun = group_by_time,
+  label = "group",
+  ...
+) {
+  if(label %in% names(data)) {
+    cli::cli_alert_warning("Overwriting {.field {label}} column in {.arg data}.")
+  }
+  if(inherits(fun, "character")) {
+    fun <- get(fun)
+  }
+  data[[label]] <- fun(data, ...)
+  data
+}
+
+#' Group data by time using bin separators
+#'
+#' @inheritParams add_grouping_column
+#' @inheritParams run_eval
+#' @param bins vector of bin separators. Suggestion to keep the last bin
+#'   separator `Inf`, to ensure all points are included in a group.
+#'
+#' @returns a numeric vector of the same length as the number of rows in `data`.
+#'
+#' @examples
+#' group_by_time(nm_busulfan)
+#'
+#' @export
+group_by_time <- function(
+  data,
+  bins = c(0, 24, 48, 72, 96, Inf),
+  dictionary = list("TIME" = "TIME"),
+  ...
+) {
+  as.numeric(cut(data[[dictionary$TIME]], breaks = bins, include.lowest = TRUE))
+}
+
+#' Will create a separate group for each dose intervals that contains at least
+#' one sample
+#'
+#' @inheritParams add_grouping_column
+#' @inheritParams run_eval
+#'
+#' @returns a numeric vector of the same length as the number of rows in `data`.
+#'
+#' @examples
+#' group_by_dose(nm_busulfan)
+#'
+#' @export
+group_by_dose <- function(
+  data,
+  dictionary = list("ID" = "ID", "EVID" = "EVID", "TIME" = "TIME"),
+  ...
+) {
+  data |>
+    dplyr::group_by(dplyr::pick(dictionary$ID)) |>
+    dplyr::mutate(`_dose_idx` = cumsum(.data[[dictionary$EVID]])) |>
+    dplyr::group_by(dplyr::pick(dictionary$ID, "_dose_idx")) |>
+    dplyr::mutate(`_has_sample` = any(.data[[dictionary$EVID]] == 0)) |>
+    dplyr::group_by(dplyr::pick(dictionary$ID)) |>
+    dplyr::mutate(
+      `_group` = cumsum(.data[[dictionary$EVID]] == 1 & .data$`_has_sample`)
+    ) |>
+    dplyr::pull(.data$`_group`)
+}

R/check_input_data.R

@@ -8,13 +8,57 @@
 #' @returns a data.frame
 check_input_data <- function(data, dictionary) {
   if(!is.null(dictionary) && length(dictionary) > 0) {
-    if(!all(unlist(dictionary) %in% names(data))) {
-      cli::cli_abort("Not all columns listed in dictionary were found in input data.")
-    }
+    check_dictionary_columns(data, dictionary)
+    check_valid_dictionary(dictionary)
+    # TODO: We should move this into a different function like tidy_input_data().
+    # Checks and data manipulation should occur separately.
     cli::cli_alert_info("Renaming dataset columns based on dictionary.")
-    new_data <- dplyr::rename(data, dplyr::all_of(dictionary))
+    data <- data |>
+      dplyr::rename(
+        !!!rlang::set_names(
+          dictionary,
+          names(dictionary)
+        )
+      )
   } else {
     cli::cli_alert_info("No data `dictionary` provided, assuming common NONMEM column names.")
   }
+  check_required_cols(data)
   data
 }
+
+check_dictionary_columns <- function(data, dictionary) {
+  has_cols <- unlist(dictionary) %in% names(data)
+  if(all(has_cols)) return()
+  missing_cols <- unlist(dictionary)[!has_cols]
+  cli::cli_abort(c(
+           "Column names in {.arg dictionary} must exist in {.arg data}.",
+    "x" = "{.str {missing_cols}} {?is/are} not {?a/} column{?s} in {.arg data}."
+  ))
+}
+
+check_valid_dictionary <- function(dictionary) {
+  valid_names <- c("ID", "TIME", "EVID", "DV")
+  has_valid_names <- names(dictionary) %in% valid_names
+  if (all(has_valid_names)) return()
+  invalid_names <- names(dictionary)[!has_valid_names]
+  cli::cli_abort(c(
+    "{.arg dictionary} names must be any of the following: {.field {valid_names}}.",
+    "x" = "{.field {invalid_names}} {?is/are} not {?a/} valid name{?s}.",
+    "i" = "Did you make a typo?"
+  ))
+}
+
+check_required_cols <- function(data) {
+  required_cols <- c("ID", "TIME", "EVID", "DV")
+  has_required <- required_cols %in% names(data)
+  if (all(has_required)) return()
+  missing_cols <- required_cols[!has_required]
+  cli::cli_abort(c(
+          "{.arg data} must have the following columns: {.field {required_cols}}.",
+    "x" = "{.field {missing_cols}} {?is/are} missing.",
+    "i" = "Did you forget to rename a column?"
+  ))
+}
+
+

R/compare_psn_results.R

@@ -48,7 +48,10 @@ compare_psn_results <- function(mipdeval, psn, tol = 0.1, .apriori) {
   # rel diff, whether all observations were within tolerance
   dplyr::summarise(
     reldiff_psn_results(mipdeval, psn, .apriori),
-    within_tol = all(abs(.data$rel_diff) < tol)
+    mean_rel_diff = mean(.data$rel_diff),
+    sd_rel_diff = stats::sd(.data$rel_diff),
+    max_rel_diff = max(abs(.data$rel_diff)),
+    within_tol = all(.data$max_rel_diff < tol)
   )
 }
 

R/parse_input_data.R

@@ -8,7 +8,8 @@
 parse_input_data <- function(
     data,
     covariates,
-    ids
+    ids,
+    group = NULL
 ) {
 
   ## Filter IDs, if needed
@@ -26,19 +27,30 @@ parse_input_data <- function(
   out <- purrr::map(
     new_data,
     parse_nm_data,
-    covariates = covariates
+    covariates = covariates,
+    group = group
   )
 
   out
 }
 
-parse_nm_data <- function(data, covariates) {
+parse_nm_data <- function(data, covariates, group = NULL) {
   out <- list()
   out$covariates <- make_covariates_object(data, covariates)
   out$regimen <- PKPDsim::nm_to_regimen(data)
-  out$observations <- data |>
+  obs <- data |>
+    dplyr::filter(.data$EVID == 0)
+  if(!is.null(group)) {
+    obs[["_grouper"]] <- obs[[group]]
+  } else { # if no grouping, add a grouper that iterates over all observations
+    obs <- obs |>
+      dplyr::group_by(.data$ID) |>
+      dplyr::mutate("_grouper" = 1:length(.data$ID)) |>
+      dplyr::ungroup()
+  }
+  out$observations <- obs |>
     dplyr::filter(.data$EVID == 0) |>
-    dplyr::select("ID", "TIME", "DV", "CMT") |>
+    dplyr::select("ID", "TIME", "DV", "CMT", "_grouper") |>
     dplyr::rename(id = "ID", t = "TIME", y = "DV", cmt = "CMT")
   out
 }

R/parse_model.R

@@ -31,25 +31,43 @@ parse_model.character <- function(
   parameters = NULL,
   ruv = NULL,
   omega = NULL,
-  iov_bins = NULL
+  fixed = NULL,
+  iov = NULL
 ) {
   check_installed_model_library(model)
+
   withr::with_package(model, quietly = TRUE, {
     # Use model defaults if user has not supplied overrides:
     if (is.null(parameters)) parameters <- parameters()
     if (is.null(ruv)) ruv <- ruv()
     if (is.null(omega)) omega <- omega_matrix()
-    if (is.null(iov_bins)) iov_bins <- iov()
-    list(
+    if (is.null(fixed)) fixed <- fixed()
+    if (is.null(iov)) iov <- iov()
+    bins <- as.numeric(iov$bins)
+    out <- list(
       model = model(),
       parameters = parameters,
-      # covariates = PKPDsim::get_model_covariates(model()),
       ruv = ruv,
-      omega_matrix = omega,
-      iov_bins = iov_bins,
-      fixed = fixed()
+      omega = omega,
+      fixed = fixed(),
+      bins = bins,
+      kappa = NULL
     )
+    if(!is.null(iov)) {
+      iov_obj <- PKPDmap::create_iov_object(
+        cv = iov$cv,
+        omega = omega,
+        bins = bins,
+        parameters = parameters,
+        fixed = fixed,
+        ruv = ruv,
+        verbose = FALSE
+      )
+      iov_updates <- c("parameters", "omega", "fixed", "bins", "kappa")
+      out[iov_updates] <- iov_obj[iov_updates]
+    }
   })
+  out
 }
 
 #' @export
@@ -60,16 +78,32 @@ parse_model.PKPDsim <- function(
   parameters,
   ruv,
   omega,
-  iov_bins = NULL
+  fixed = NULL,
+  iov = NULL
 ) {
-  list(
+  out <- list(
     model = model,
     parameters = parameters,
     ruv = ruv,
-    omega_matrix = omega,
-    iov_bins = iov_bins,
-    fixed = attr(model, "fixed")
+    omega = omega,
+    fixed = fixed,
+    bins = numeric(0),
+    kappa = NULL
   )
+  if(!is.null(iov)) {
+    iov_obj <- PKPDmap::create_iov_object(
+      cv = iov$cv,
+      omega = omega,
+      bins = as.numeric(iov$bins),
+      parameters = parameters,
+      fixed = fixed,
+      ruv = ruv,
+      verbose = FALSE
+    )
+    iov_updates <- c("parameters", "omega", "fixed", "bins", "kappa")
+    out[iov_updates] <- iov_obj[iov_updates]
+  }
+  out
 }
 
 #' Check if PKPDsim model library is installed

R/parse_psn_proseval_results.R

@@ -1,22 +1,46 @@
 #' Parse PsN::proseval results.csv to filter out only the rows
-#' that we need (for prediction of next level)
+#' that we need (for prediction of next sample or group of samples)
 #'
 #' @param data results data.frame or path to csv file
+#' @param group optional. Group samples using the column name specified in
+#'   `group`
 #'
 #' @returns A data.frame
 #' @export
-parse_psn_proseval_results <- function(data) {
+parse_psn_proseval_results <- function(data, group = NULL) {
   if(inherits(data, "character")) {
     res <- read.csv(file = data)
   } else {
     res <- data
   }
-  res |>
+  res <- res |>
     dplyr::group_by(.data$ID) |>
-    dplyr::filter(.data$EVID == 2 & .data$DV > 0) |>
-    dplyr::filter(!duplicated(.data$OBS)) |>
-    dplyr::mutate(n = 1:length(.data$OBS)) |>
+    dplyr::filter(.data$EVID == 2 & .data$DV > 0)
+  if(is.null(group)) { # assume we just want to look at next observation
+    group <- "group"
+    res <- res |>
+      dplyr::group_by(.data$OBS, .add = TRUE) |>
+      dplyr::mutate(group = seq(.data$OBS[1] + 1, .data$OBS[1] + dplyr::n(), 1)) |>
+      dplyr::ungroup()
+  }
+  ## create a data.frame on which OBS iterations to take from proseval results
+  idx_tab <- res |>
+    dplyr::filter(.data$OBS == 1) |>
+    dplyr::group_by(.data$ID) |>
+    dplyr::mutate(idx = 1:dplyr::n()) |>
+    dplyr::filter(!duplicated(group)) |>
+    dplyr::ungroup() |>
+    dplyr::filter(group > 1) |>
+    dplyr::select("ID", "idx") |>
+    dplyr::group_by(.data$ID) |>
+    dplyr::summarise(idx = list(.data$idx))
+  dplyr::left_join(res, idx_tab) |>
+    dplyr::filter( # filter rows that match the index table
+      purrr::map2_lgl(.data$idx, .data$OBS, ~ .y %in% .x)
+    ) |>
+    dplyr::group_by(.data$ID, .data$OBS) |>
+    dplyr::filter(.data$group == .data$group[1]) |>
     dplyr::mutate(res = .data$DV - .data$IPRED) |>
-    dplyr::select("TIME", obs = "DV", pred = "IPRED", pop = "PRED", "res", "n") |>
+    dplyr::select("TIME", obs = "DV", pred = "IPRED", pop = "PRED", "res", !!group) |>
     dplyr::ungroup()
 }