From 7c6ed0e83a577f17f2563b9f60edb8b1c1716183 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Sat, 16 Sep 2023 22:31:29 +0200 Subject: [PATCH] 675: Add `DesignGrouped` (#677) Co-authored-by: github-actions <41898282+github-actions[bot]@users.noreply.github.com> --- .lintr | 1 + DESCRIPTION | 1 + NAMESPACE | 4 + NEWS.md | 2 +- R/Design-class.R | 87 +++++ R/Design-methods.R | 308 ++++++++++-------- R/Design-validity.R | 16 + R/Rules-methods.R | 60 ++-- R/crmPack-package.R | 3 +- R/helpers_design.R | 168 ++++++++++ _pkgdown.yaml | 2 + examples/Design-class-DesignGrouped.R | 55 ++++ .../Design-method-simulate-DesignGrouped.R | 74 +++++ man/DesignGrouped-class.Rd | 124 +++++++ man/get_result_list.Rd | 6 +- man/h_add_dlts.Rd | 26 ++ man/nextBest.Rd | 2 +- man/set_seed.Rd | 5 +- man/simulate-DesignGrouped-method.Rd | 136 ++++++++ man/v_design.Rd | 6 + tests/testthat/helper-model.R | 19 -- tests/testthat/test-Design-class.R | 44 +++ tests/testthat/test-Design-methods.R | 193 ++++++++--- tests/testthat/test-Design-validity.R | 21 ++ tests/testthat/test-Rules-methods.R | 12 +- tests/testthat/test-helpers.R | 2 +- tests/testthat/test-helpers_design.R | 100 ++++++ 27 files changed, 1225 insertions(+), 252 deletions(-) create mode 100644 R/helpers_design.R create mode 100644 examples/Design-class-DesignGrouped.R create mode 100644 examples/Design-method-simulate-DesignGrouped.R create mode 100644 man/DesignGrouped-class.Rd create mode 100644 man/h_add_dlts.Rd create mode 100644 man/simulate-DesignGrouped-method.Rd create mode 100644 tests/testthat/test-helpers_design.R diff --git a/.lintr b/.lintr index 7049743a2..e1cd584ba 100644 --- a/.lintr +++ b/.lintr @@ -1,6 +1,7 @@ linters: linters_with_defaults( line_length_linter = line_length_linter(120), cyclocomp_linter = NULL, + indentation_linter = NULL, object_usage_linter = NULL, object_length_linter = NULL, object_name_linter = object_name_linter(c("CamelCase", "camelCase", "snake_case")) diff --git a/DESCRIPTION b/DESCRIPTION index c2c85dc83..fcca6d997 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -110,6 +110,7 @@ Collate: 'checkmate.R' 'crmPack-package.R' 'helpers_covr.R' + 'helpers_design.R' 'helpers_model.R' 'logger.R' 'utils-pipe.R' diff --git a/NAMESPACE b/NAMESPACE index 22a2fa797..b6e7c10b3 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -29,6 +29,7 @@ export(.DefaultCohortSizeParts) export(.DefaultCohortSizeRange) export(.DefaultDALogisticLogNormal) export(.DefaultDataGrouped) +export(.DefaultDesignGrouped) export(.DefaultDualEndpoint) export(.DefaultDualEndpointBeta) export(.DefaultDualEndpointEmax) @@ -87,6 +88,7 @@ export(.DefaultStoppingTargetBiomarker) export(.DefaultStoppingTargetProb) export(.DefaultTITELogisticLogNormal) export(.Design) +export(.DesignGrouped) export(.DualDesign) export(.DualEndpoint) export(.DualEndpointBeta) @@ -191,6 +193,7 @@ export(DataMixture) export(DataOrdinal) export(DataParts) export(Design) +export(DesignGrouped) export(DualDesign) export(DualEndpoint) export(DualEndpointBeta) @@ -378,6 +381,7 @@ exportClasses(DataMixture) exportClasses(DataOrdinal) exportClasses(DataParts) exportClasses(Design) +exportClasses(DesignGrouped) exportClasses(DualDesign) exportClasses(DualEndpoint) exportClasses(DualEndpointBeta) diff --git a/NEWS.md b/NEWS.md index 25af45d59..67e80121a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,5 @@ # Version 1.0.9000.9133 -* Added new `DataGrouped` class to support simultaneous dose escalation with monotherapy and combination therapy. +* Added new `DataGrouped` and `DesignGrouped` classes with corresponding model `LogisticLogNormalGrouped` to support simultaneous dose escalation with monotherapy and combination therapy arms. * Created the `CrmPackClass` class as the ultimate ancestor of all other `crmPack` classes to allow identification of crmPack classes and simpler definition of generic methods. diff --git a/R/Design-class.R b/R/Design-class.R index a0687ccbc..05db7632b 100644 --- a/R/Design-class.R +++ b/R/Design-class.R @@ -529,3 +529,90 @@ DADesign <- function(model, data, safetyWindow = safetyWindow ) } + +# DesignGrouped ---- + +## class ---- + +#' `DesignGrouped` +#' +#' @description `r lifecycle::badge("experimental")` +#' +#' [`DesignGrouped`] combines two [`Design`] objects: one for the mono and one +#' for the combo arm of a joint dose escalation design. +#' +#' @slot model (`LogisticLogNormalGrouped`)\cr the model to be used, currently only one +#' class is allowed. +#' @slot mono (`Design`)\cr defines the dose escalation rules for the mono arm, see +#' details. +#' @slot combo (`Design`)\cr defines the dose escalation rules for the combo arm, see +#' details. +#' @slot first_cohort_mono_only (`flag`)\cr whether first test one mono agent cohort, and then +#' once its DLT data has been collected, we proceed from the second cohort onwards with +#' concurrent mono and combo cohorts. +#' @slot same_dose (`flag`)\cr whether the lower dose of the separately determined mono and combo +#' doses should be used as the next dose for both mono and combo. +#' +#' @details Note that the model slots inside the `mono` and `combo` parameters +#' are ignored (because we don't fit separate regression models for the mono and +#' combo arms). Instead, the `model` parameter is used to fit a joint regression +#' model for the mono and combo arms together. +#' +#' @aliases DesignGrouped +#' @export +#' +.DesignGrouped <- setClass( + Class = "DesignGrouped", + slots = c( + model = "LogisticLogNormalGrouped", + mono = "Design", + combo = "Design", + first_cohort_mono_only = "logical", + same_dose = "logical" + ), + prototype = prototype( + model = .DefaultLogisticLogNormalGrouped(), + mono = .Design(), + combo = .Design(), + first_cohort_mono_only = TRUE, + same_dose = TRUE + ), + validity = v_design_grouped, + contains = "CrmPackClass" +) + +## constructor ---- + +#' @rdname DesignGrouped-class +#' +#' @param model (`LogisticLogNormalGrouped`)\cr see slot definition. +#' @param mono (`Design`)\cr see slot definition. +#' @param combo (`Design`)\cr see slot definition. +#' @param first_cohort_mono_only (`flag`)\cr see slot definition. +#' @param same_dose (`flag`)\cr see slot definition. +#' @param ... not used. +#' +#' @export +#' @example examples/Design-class-DesignGrouped.R +#' +DesignGrouped <- function(model, + mono, + combo = mono, + first_cohort_mono_only = TRUE, + same_dose = TRUE, + ...) { + .DesignGrouped( + model = model, + mono = mono, + combo = combo, + first_cohort_mono_only = first_cohort_mono_only, + same_dose = same_dose + ) +} + +## default constructor ---- + +#' @rdname DesignGrouped-class +#' @note Typically, end-users will not use the `.DefaultDesignGrouped()` function. +#' @export +.DefaultDesignGrouped <- .DesignGrouped diff --git a/R/Design-methods.R b/R/Design-methods.R index 03e9b5be1..7d07790b1 100644 --- a/R/Design-methods.R +++ b/R/Design-methods.R @@ -7,147 +7,6 @@ #' @include mcmc.R NULL -# helper functions ---- - -## set_seed ---- - -#' Helper function to set and save the RNG seed -#' -#' @description `r lifecycle::badge("stable")` -#' -#' This code is basically copied from `stats:::simulate.lm`. -#' -#' @param seed an object specifying if and how the random number generator -#' should be initialized ("seeded"). Either `NULL` (default) or an -#' integer that will be used in a call to [set.seed()] before -#' simulating the response vectors. If set, the value is saved as the -#' `seed` slot of the returned object. The default, `NULL` will -#' not change the random generator state. -#' @return The integer vector containing the random number generate state will -#' be returned, in order to call this function with this input to reproduce -#' the obtained simulation results. -#' -#' @export -#' @keywords programming -set_seed <- function(seed = NULL) { - assert_number(seed, null.ok = TRUE) - - if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) { - runif(1) - } - - if (is.null(seed)) { - get(".Random.seed", envir = .GlobalEnv) - } else { - seed <- as.integer(seed) - r_seed <- get(".Random.seed", envir = .GlobalEnv) - # Make sure r_seed exists in parent frame. - assign(".r_seed", r_seed, envir = parent.frame()) - set.seed(seed) - # Here we need the r_seed in the parent.frame! - do.call( - "on.exit", - list(quote(assign(".Random.seed", .r_seed, envir = .GlobalEnv))), - envir = parent.frame() - ) - structure(seed, kind = as.list(RNGkind())) - } -} - -## get_result_list ---- - -#' Helper function to obtain simulation results list -#' -#' @description `r lifecycle::badge("stable")` -#' -#' The function `fun` can use variables that are visible to itself. -#' The names of these variables have to be given in the vector `vars`. -#' -#' @param fun (`function`)\cr the simulation function for a single iteration, which takes as -#' single parameter the iteration index. -#' @param nsim number of simulations to be conducted. -#' @param vars names of the variables. -#' @param parallel should the simulation runs be parallelized across the -#' clusters of the computer? -#' @param n_cores how many cores should be used for parallel computing? -#' @return The list with all simulation results (one iteration corresponds -#' to one list element). -#' -#' @importFrom parallel makeCluster -#' @importFrom parallelly availableCores -#' @keywords internal programming -get_result_list <- function( - fun, - nsim, - vars, - parallel, - n_cores) { - assert_flag(parallel) - assert_integerish(n_cores, lower = 1) - - if (!parallel) { - lapply( - X = seq_len(nsim), - FUN = fun - ) - } else { - # Process all simulations. - cores <- min( - safeInteger(n_cores), - parallelly::availableCores() - ) - - # Start the cluster. - cl <- parallel::makeCluster(cores) - - # Load the required R package. - parallel::clusterEvalQ(cl, { - library(crmPack) - NULL - }) - - # Export local variables from the caller environment. - # Note: parent.frame() is different from parent.env() which returns - # the environment where this function has been defined! - parallel::clusterExport( - cl = cl, - varlist = vars, - envir = parent.frame() - ) - - # Export all global variables. - parallel::clusterExport( - cl = cl, - varlist = ls(.GlobalEnv) - ) - - # Load user extensions from global options. - crmpack_extensions <- getOption("crmpack_extensions") - if (is.null(crmpack_extensions) != TRUE) { - tryCatch( - { - parallel::clusterCall(cl, crmpack_extensions) - }, - error = function(e) { - stop("Failed to export crmpack_extensions: ", e$message) - } - ) - } - - # Do the computations in parallel. - res <- parallel::parLapply( - cl = cl, - X = seq_len(nsim), - fun = fun - ) - - # Stop the cluster. - parallel::stopCluster(cl) - - res - } -} - # nolint start ## ============================================================ @@ -4708,3 +4567,170 @@ setMethod("simulate", ## -------------------------------------------------------------------------- # nolint end + +# simulate ---- + +## DesignGrouped ---- + +#' Simulate Method for the [`DesignGrouped`] Class +#' +#' @description `r lifecycle::badge("experimental")` +#' +#' A simulate method for [`DesignGrouped`] designs. +#' +#' @param object (`DesignGrouped`)\cr the design we want to simulate trials from. +#' @param nsim (`number`)\cr how many trials should be simulated. +#' @param seed (`RNGstate`)\cr generated with [set_seed()]. +#' @param truth (`function`)\cr a function which takes as input a dose (vector) and +#' returns the true probability (vector) for toxicity for the mono arm. +#' Additional arguments can be supplied in `args`. +#' @param combo_truth (`function`)\cr same as `truth` but for the combo arm. +#' @param args (`data.frame`)\cr optional `data.frame` with arguments that work +#' for both the `truth` and `combo_truth` functions. The column names correspond to +#' the argument names, the rows to the values of the arguments. The rows are +#' appropriately recycled in the `nsim` simulations. +#' @param firstSeparate (`flag`)\cr whether to enroll the first patient separately +#' from the rest of the cohort and close the cohort in case a DLT occurs in this +#' first patient. +#' @param mcmcOptions (`McmcOptions`)\cr MCMC options for each evaluation in the trial. +#' @param parallel (`flag`)\cr whether the simulation runs are parallelized across the +#' cores of the computer. +#' @param nCores (`number`)\cr how many cores should be used for parallel computing. +#' @param ... not used. +#' +#' @return A list of `mono` and `combo` simulation results as [`Simulations`] objects. +#' +#' @aliases simulate-DesignGrouped +#' @export +#' @example examples/Design-method-simulate-DesignGrouped.R +#' +setMethod( + "simulate", + signature = + signature( + object = "DesignGrouped", + nsim = "ANY", + seed = "ANY" + ), + def = + function(object, + nsim = 1L, + seed = NULL, + truth, + combo_truth, + args = data.frame(), + firstSeparate = FALSE, + mcmcOptions = McmcOptions(), + parallel = FALSE, + nCores = min(parallelly::availableCores(), 5), + ...) { + nsim <- safeInteger(nsim) + assert_function(truth) + assert_function(combo_truth) + assert_data_frame(args) + assert_count(nsim, positive = TRUE) + assert_flag(firstSeparate) + assert_flag(parallel) + assert_count(nCores, positive = TRUE) + + n_args <- max(nrow(args), 1L) + rng_state <- set_seed(seed) + sim_seeds <- sample.int(n = 2147483647, size = nsim) + + run_sim <- function(iter_sim) { + set.seed(sim_seeds[iter_sim]) + current <- list(mono = list(), combo = list()) + # Define true toxicity functions. + current$args <- args[(iter_sim - 1) %% n_args + 1, , drop = FALSE] + current$mono$truth <- function(dose) do.call(truth, c(dose, current$args)) + current$combo$truth <- function(dose) do.call(combo_truth, c(dose, current$args)) + # Start the simulated data with the provided one. + current$mono$data <- object@mono@data + current$combo$data <- object@combo@data + # We are in the first cohort and continue for mono and combo. + current$first <- TRUE + current$mono$stop <- current$combo$stop <- FALSE + # What are the next doses to be used? Initialize with starting doses. + if (object@same_dose) { + current$mono$dose <- current$combo$dose <- min(object@mono@startingDose, object@combo@startingDose) + } else { + current$mono$dose <- object@mono@startingDose + current$combo$dose <- object@combo@startingDose + } + # Inside this loop we simulate the whole trial, until stopping. + while (!(current$mono$stop && current$combo$stop)) { + if (!current$mono$stop) { + current$mono$data <- current$mono$data |> + h_add_dlts(current$mono$dose, current$mono$truth, object@mono@cohort_size, firstSeparate) + } + if (!current$combo$stop && (!current$first || !object@first_cohort_mono_only)) { + current$combo$data <- current$combo$data |> + h_add_dlts(current$combo$dose, current$combo$truth, object@combo@cohort_size, firstSeparate) + } + if (current$first) current$first <- FALSE + current$grouped <- h_group_data(current$mono$data, current$combo$data) + current$samples <- mcmc(current$grouped, object@model, mcmcOptions) + if (!current$mono$stop) { + current$mono$limit <- maxDose(object@mono@increments, data = current$mono$data) + current$mono$dose <- object@mono@nextBest |> + nextBest(current$mono$limit, current$samples, object@model, current$grouped, group = "mono") + current$mono$dose <- current$mono$dose$value + current$mono$stop <- object@mono@stopping |> + stopTrial(current$mono$dose, current$samples, object@model, current$mono$data, group = "mono") + current$mono$results <- h_unpack_stopit(current$mono$stop) + } + if (!current$combo$stop) { + current$combo$limit <- if (is.na(current$mono$dose)) { + 0 + } else { + maxDose(object@combo@increments, current$combo$data) |> + min(current$mono$dose, na.rm = TRUE) + } + current$combo$dose <- object@combo@nextBest |> + nextBest(current$combo$limit, current$samples, object@model, current$grouped, group = "combo") + current$combo$dose <- current$combo$dose$value + current$combo$stop <- object@combo@stopping |> + stopTrial(current$combo$dose, current$samples, object@model, current$combo$data, group = "combo") + current$combo$results <- h_unpack_stopit(current$combo$stop) + } + if (object@same_dose && !current$mono$stop && !current$combo$stop) { + current$mono$dose <- current$combo$dose <- min(current$mono$dose, current$combo$dose) + } + } + current$mono$fit <- fit(current$samples, object@model, current$grouped, group = "mono") + current$combo$fit <- fit(current$samples, object@model, current$grouped, group = "combo") + lapply( + X = current[c("mono", "combo")], FUN = with, + list( + data = data, dose = dose, fit = subset(fit, select = -dose), + stop = attr(stop, "message"), results = results + ) + ) + } + vars_needed <- c("simSeeds", "args", "nArgs", "truth", "combo_truth", "firstSeparate", "object", "mcmcOptions") + result_list <- get_result_list(run_sim, nsim, vars_needed, parallel, nCores) + # Now we have a list with each element containing mono and combo. Reorder this a bit: + result_list <- list( + mono = lapply(result_list, "[[", "mono"), + combo = lapply(result_list, "[[", "combo") + ) + # Put everything in a list with both mono and combo Simulations: + lapply(result_list, function(this_list) { + data_list <- lapply(this_list, "[[", "data") + recommended_doses <- as.numeric(sapply(this_list, "[[", "dose")) + fit_list <- lapply(this_list, "[[", "fit") + stop_reasons <- lapply(this_list, "[[", "stop") + report_results <- lapply(this_list, "[[", "results") + stop_report <- as.matrix(do.call(rbind, report_results)) + + Simulations( + data = data_list, + doses = recommended_doses, + fit = fit_list, + stop_reasons = stop_reasons, + stop_report = stop_report, + seed = rng_state + ) + }) + } +) diff --git a/R/Design-validity.R b/R/Design-validity.R index d20b0f4ca..e5e53be4c 100644 --- a/R/Design-validity.R +++ b/R/Design-validity.R @@ -27,3 +27,19 @@ v_rule_design <- function(object) { ) v$result() } + + +#' @describeIn v_design validates that the [`DesignGrouped`] object +#' contains valid flags. +v_design_grouped <- function(object) { + v <- Validate() + v$check( + test_flag(object@first_cohort_mono_only), + "first_cohort_mono_only must be a flag" + ) + v$check( + test_flag(object@same_dose), + "same_dose must be a flag" + ) + v$result() +} diff --git a/R/Rules-methods.R b/R/Rules-methods.R index 8d5d32127..dd3af7367 100644 --- a/R/Rules-methods.R +++ b/R/Rules-methods.R @@ -67,7 +67,7 @@ setMethod( ), definition = function(nextBest, doselimit = Inf, samples, model, data, ...) { # Generate the MTD samples and derive the next best dose. - dose_target_samples <- dose(x = nextBest@target, model, samples) + dose_target_samples <- dose(x = nextBest@target, model, samples, ...) dose_target <- nextBest@derive(dose_target_samples) # Round to the next possible grid point. @@ -292,7 +292,7 @@ setMethod("nextBest", ), definition = function(nextBest, doselimit = Inf, samples, model, data, ...) { # Matrix with samples from the dose-tox curve at the dose grid points. - prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples) + prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples, ...) # Compute probabilities to be in target and overdose tox interval. prob_underdosing <- colMeans(prob_samples < nextBest@target[1]) prob_target <- colMeans(h_in_range(prob_samples, nextBest@target)) @@ -553,7 +553,7 @@ setMethod( ), definition = function(nextBest, doselimit = Inf, samples, model, data, ...) { # Matrix with samples from the dose-tox curve at the dose grid points. - prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples) + prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples, ...) dlt_prob <- colMeans(prob_samples) # Determine the dose with the closest distance. @@ -651,7 +651,7 @@ setMethod( ), definition = function(nextBest, doselimit = Inf, samples, model, data, ...) { # Matrix with samples from the dose-tox curve at the dose grid points. - prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples) + prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples, ...) criterion <- colMeans(h_info_theory_dist(prob_samples, nextBest@target, nextBest@asymmetry)) @@ -697,8 +697,8 @@ setMethod( # Target dose estimates, i.e. the dose with probability of the occurrence of # a DLT that equals to the prob_target_drt or prob_target_eot. - dose_target_drt <- dose(x = prob_target_drt, model) - dose_target_eot <- dose(x = prob_target_eot, model) + dose_target_drt <- dose(x = prob_target_drt, model, ...) + dose_target_eot <- dose(x = prob_target_eot, model, ...) # Find the next best doses in the doseGrid. The next best dose is the dose # at level closest and below the target dose estimate. @@ -732,7 +732,7 @@ setMethod( prob_target_eot = prob_target_eot, dose_target_eot = dose_target_eot, data = data, - prob_dlt = prob(dose = data@doseGrid, model = model), + prob_dlt = prob(dose = data@doseGrid, model = model, ...), doselimit = doselimit, next_dose = next_dose_drt ) @@ -777,12 +777,12 @@ setMethod( model = "LogisticIndepBeta", data = "Data" ), - definition = function(nextBest, doselimit = Inf, samples, model, data, ...) { + definition = function(nextBest, doselimit = Inf, samples, model, data, in_sim, ...) { # Generate target dose samples, i.e. the doses with probability of the # occurrence of a DLT that equals to the nextBest@prob_target_drt # (or nextBest@prob_target_eot, respectively). - dose_target_drt_samples <- dose(x = nextBest@prob_target_drt, model, samples) - dose_target_eot_samples <- dose(x = nextBest@prob_target_eot, model, samples) + dose_target_drt_samples <- dose(x = nextBest@prob_target_drt, model, samples, ...) + dose_target_eot_samples <- dose(x = nextBest@prob_target_eot, model, samples, ...) # Derive the prior/posterior estimates based on two above samples. dose_target_drt <- nextBest@derive(dose_target_drt_samples) @@ -865,15 +865,15 @@ setMethod( # Target dose estimates, i.e. the dose with probability of the occurrence of # a DLT that equals to the prob_target_drt or prob_target_eot. - dose_target_drt <- dose(x = prob_target_drt, model) - dose_target_eot <- dose(x = prob_target_eot, model) + dose_target_drt <- dose(x = prob_target_drt, model, ...) + dose_target_eot <- dose(x = prob_target_eot, model, ...) # Find the dose which gives the maximum gain. dosegrid_range <- dose_grid_range(data) opt <- optim( par = dosegrid_range[1], fn = function(DOSE) { - -gain(DOSE, model_dle = model, model_eff = model_eff) + -gain(DOSE, model_dle = model, model_eff = model_eff, ...) }, method = "L-BFGS-B", lower = dosegrid_range[1], @@ -986,15 +986,15 @@ setMethod( # Generate target dose samples, i.e. the doses with probability of the # occurrence of a DLT that equals to the prob_target_drt or prob_target_eot. - dose_target_drt_samples <- dose(x = prob_target_drt, model, samples = samples) - dose_target_eot_samples <- dose(x = prob_target_eot, model, samples = samples) + dose_target_drt_samples <- dose(x = prob_target_drt, model, samples = samples, ...) + dose_target_eot_samples <- dose(x = prob_target_eot, model, samples = samples, ...) # Derive the prior/posterior estimates based on two above samples. dose_target_drt <- nextBest@derive(dose_target_drt_samples) dose_target_eot <- nextBest@derive(dose_target_eot_samples) # Gain samples. - gain_samples <- sapply(data@doseGrid, gain, model, samples, model_eff, samples_eff) + gain_samples <- sapply(data@doseGrid, gain, model, samples, model_eff, samples_eff, ...) # For every sample, get the dose (from the dose grid) that gives the maximum gain value. dose_lev_mg_samples <- apply(gain_samples, 1, which.max) dose_mg_samples <- data@doseGrid[dose_lev_mg_samples] @@ -1085,7 +1085,7 @@ setMethod( ), definition = function(nextBest, doselimit, samples, model, data, ...) { # Matrix with samples from the dose-tox curve at the dose grid points. - prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples) + prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples, ...) # Determine the maximum dose level with a toxicity probability below or # equal to the target and calculate how often a dose is selected as MTD @@ -1237,7 +1237,7 @@ setMethod( ), definition = function(nextBest, doselimit, samples, model, data, ...) { # Matrix with samples from the dose-tox curve at the dose grid points. - prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples) + prob_samples <- sapply(data@doseGrid, prob, model = model, samples = samples, ...) # Determine which dose level has the minimum distance to target. dose_min_mtd_dist <- apply( @@ -2284,7 +2284,8 @@ setMethod("stopTrial", mtdSamples <- dose( x = stopping@target, model, - samples + samples, + ... ) ## what is the absolute threshold? @@ -2354,7 +2355,8 @@ setMethod( mtd_samples <- dose( x = stopping@target, model, - samples + samples, + ... ) # CV of MTD expressed as percentage, derived based on MTD posterior samples. mtd_cv <- (mad(mtd_samples) / median(mtd_samples)) * 100 @@ -2459,7 +2461,7 @@ setMethod( definition = function(stopping, dose, samples, model, data, ...) { # Compute the target biomarker prob at this dose. # Get the biomarker level samples at the dose grid points. - biom_level_samples <- biomarker(xLevel = seq_len(data@nGrid), model, samples) + biom_level_samples <- biomarker(xLevel = seq_len(data@nGrid), model, samples, ...) # If target is relative to maximum. if (stopping@is_relative) { @@ -2910,7 +2912,8 @@ setMethod( dose_target_samples <- dose( x = stopping@prob_target, model = model, - samples = samples + samples = samples, + ... ) # 95% credibility interval. dose_target_ci <- quantile(dose_target_samples, probs = c(0.025, 0.975)) @@ -2955,7 +2958,7 @@ setMethod("stopTrial", assert_probability(stopping@prob_target) prob_target <- stopping@prob_target - dose_target_samples <- dose(x = prob_target, model = model) + dose_target_samples <- dose(x = prob_target, model = model, ...) ## Find the variance of the log of the dose_target_samples(eta) M1 <- matrix(c(-1 / (model@phi2), -(log(prob_target / (1 - prob_target)) - model@phi1) / (model@phi2)^2), 1, 2) M2 <- model@Pcov @@ -3024,7 +3027,8 @@ setMethod("stopTrial", TDtargetEndOfTrialSamples <- dose( x = prob_target, model = model, - samples = samples + samples = samples, + ... ) ## Find the TDtarget End of trial estimate TDtargetEndOfTrialEstimate <- TDderive(TDtargetEndOfTrialSamples) @@ -3047,7 +3051,8 @@ setMethod("stopTrial", model, samples, Effmodel, - Effsamples + Effsamples, + ... ) } @@ -3137,12 +3142,13 @@ setMethod("stopTrial", ## find the TDtarget End of Trial TDtargetEndOfTrial <- dose( x = prob_target, - model = model + model = model, + ... ) ## Find the dose with maximum gain value Gainfun <- function(DOSE) { - -gain(DOSE, model_dle = model, model_eff = Effmodel) + -gain(DOSE, model_dle = model, model_eff = Effmodel, ...) } # if(data@placebo) { diff --git a/R/crmPack-package.R b/R/crmPack-package.R index 44fbf1ee0..a07ec3797 100644 --- a/R/crmPack-package.R +++ b/R/crmPack-package.R @@ -138,7 +138,8 @@ globalVariables(c( "comp", "X", "skel_probs", - "is_combo" + "is_combo", + "results" )) # nolint end diff --git a/R/helpers_design.R b/R/helpers_design.R new file mode 100644 index 000000000..b30026b78 --- /dev/null +++ b/R/helpers_design.R @@ -0,0 +1,168 @@ +#' Helper Function to Set and Save the RNG Seed +#' +#' @description `r lifecycle::badge("stable")` +#' +#' This code is basically copied from `stats:::simulate.lm`. +#' +#' @param seed an object specifying if and how the random number generator +#' should be initialized ("seeded"). Either `NULL` (default) or an +#' integer that will be used in a call to [set.seed()] before +#' simulating the response vectors. If set, the value is saved as the +#' `seed` slot of the returned object. The default, `NULL` will +#' not change the random generator state. +#' @return The integer vector containing the random number generate state will +#' be returned, in order to call this function with this input to reproduce +#' the obtained simulation results. +#' +#' @export +set_seed <- function(seed = NULL) { + assert_number(seed, null.ok = TRUE) + + if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) { + runif(1) + } + + if (is.null(seed)) { + get(".Random.seed", envir = .GlobalEnv) + } else { + seed <- as.integer(seed) + r_seed <- get(".Random.seed", envir = .GlobalEnv) + # Make sure r_seed exists in parent frame. + assign(".r_seed", r_seed, envir = parent.frame()) + set.seed(seed) + # Here we need the r_seed in the parent.frame! + do.call( + "on.exit", + list(quote(assign(".Random.seed", .r_seed, envir = .GlobalEnv))), + envir = parent.frame() + ) + structure(seed, kind = as.list(RNGkind())) + } +} + +#' Helper Function to Obtain Simulation Results List +#' +#' The function `fun` can use variables that are visible to itself. +#' The names of these variables have to be given in the vector `vars`. +#' +#' @param fun (`function`)\cr the simulation function for a single iteration, which takes as +#' single parameter the iteration index. +#' @param nsim number of simulations to be conducted. +#' @param vars names of the variables. +#' @param parallel should the simulation runs be parallelized across the +#' clusters of the computer? +#' @param n_cores how many cores should be used for parallel computing? +#' @return The list with all simulation results (one iteration corresponds +#' to one list element). +#' +#' @importFrom parallel makeCluster +#' @importFrom parallelly availableCores +#' @keywords internal programming +get_result_list <- function( + fun, + nsim, + vars, + parallel, + n_cores) { + assert_flag(parallel) + assert_integerish(n_cores, lower = 1) + + if (!parallel) { + lapply( + X = seq_len(nsim), + FUN = fun + ) + } else { + # Process all simulations. + cores <- min( + safeInteger(n_cores), + parallelly::availableCores() + ) + + # Start the cluster. + cl <- parallel::makeCluster(cores) + + # Load the required R package. + parallel::clusterEvalQ(cl, { + library(crmPack) + NULL + }) + + # Export local variables from the caller environment. + # Note: parent.frame() is different from parent.env() which returns + # the environment where this function has been defined! + parallel::clusterExport( + cl = cl, + varlist = vars, + envir = parent.frame() + ) + + # Export all global variables. + parallel::clusterExport( + cl = cl, + varlist = ls(.GlobalEnv) + ) + + # Load user extensions from global options. + crmpack_extensions <- getOption("crmpack_extensions") + if (is.null(crmpack_extensions) != TRUE) { + tryCatch( + { + parallel::clusterCall(cl, crmpack_extensions) + }, + error = function(e) { + stop("Failed to export crmpack_extensions: ", e$message) + } + ) + } + + # Do the computations in parallel. + res <- parallel::parLapply( + cl = cl, + X = seq_len(nsim), + fun = fun + ) + + # Stop the cluster. + parallel::stopCluster(cl) + + res + } +} + +#' Helper Function to Add Randomly Generated DLTs During Simulations +#' +#' @param data (`Data`)\cr what data to start from. +#' @param dose (`number`)\cr current dose. +#' @param truth (`function`)\cr defines the true probability for a DLT at a dose. +#' @param cohort_size (`CohortSize`)\cr the cohort size rule to use. +#' @param first_separate (`flag`)\cr whether the first patient is enrolled separately. +#' +#' @return The updated `data`. +#' +#' @keywords internal +h_add_dlts <- function(data, + dose, + truth, + cohort_size, + first_separate) { + assert_class(data, "Data") + assert_number(dose) + assert_function(truth) + assert_class(cohort_size, "CohortSize") + assert_flag(first_separate) + + prob <- truth(dose) + size <- size(cohort_size, dose = dose, data = data) + dlts <- if (first_separate && size > 1) { + first_dlts <- rbinom(n = 1, size = 1, prob = prob) + if (first_dlts == 0) { + c(first_dlts, rbinom(n = size - 1, size = 1, prob = prob)) + } else { + first_dlts + } + } else { + rbinom(n = size, size = 1, prob = prob) + } + update(data, x = dose, y = dlts) +} diff --git a/_pkgdown.yaml b/_pkgdown.yaml index a8f2f343b..33dc7e0f6 100644 --- a/_pkgdown.yaml +++ b/_pkgdown.yaml @@ -129,6 +129,7 @@ reference: - DualDesign - TDsamplesDesign - TDDesign + - DesignGrouped - title: Internal Helper Functions contents: - h_blind_plot_data @@ -453,6 +454,7 @@ reference: - simulate,RuleDesign-method - simulate,TDDesign-method - simulate,TDsamplesDesign-method + - simulate-DesignGrouped - summary,DualSimulations-method - summary,GeneralSimulations-method - summary,PseudoDualFlexiSimulations-method diff --git a/examples/Design-class-DesignGrouped.R b/examples/Design-class-DesignGrouped.R new file mode 100644 index 000000000..d116306ca --- /dev/null +++ b/examples/Design-class-DesignGrouped.R @@ -0,0 +1,55 @@ +empty_data <- Data(doseGrid = c(1, 3, 5, 10, 15, 20, 25, 40, 50, 80, 100)) + +# Initialize the joint model. +my_model <- LogisticLogNormalGrouped( + mean = c(-0.85, 0, 1, 0), + cov = diag(1, 4), + ref_dose = 56 +) + +# Choose the rule for selecting the next dose. +my_next_best <- NextBestNCRM( + target = c(0.2, 0.35), + overdose = c(0.35, 1), + max_overdose_prob = 0.25 +) + +# Choose the rule for the cohort-size. +my_size1 <- CohortSizeRange( + intervals = c(0, 30), + cohort_size = c(1, 3) +) +my_size2 <- CohortSizeDLT( + intervals = c(0, 1), + cohort_size = c(1, 3) +) +my_size <- maxSize(my_size1, my_size2) + +# Choose the rule for stopping. +my_stopping1 <- StoppingMinCohorts(nCohorts = 3) +my_stopping2 <- StoppingTargetProb( + target = c(0.2, 0.35), + prob = 0.5 +) +my_stopping3 <- StoppingMinPatients(nPatients = 20) +my_stopping <- (my_stopping1 & my_stopping2) | my_stopping3 + +# Choose the rule for dose increments. +my_increments <- IncrementsRelative( + intervals = c(0, 20), + increments = c(1, 0.33) +) + +# Initialize the design. +design <- DesignGrouped( + model = my_model, + mono = Design( + model = .DefaultModelLogNormal(), # Ignored. + nextBest = my_next_best, + stopping = my_stopping, + increments = my_increments, + cohort_size = my_size, + data = empty_data, + startingDose = 3 + ) +) diff --git a/examples/Design-method-simulate-DesignGrouped.R b/examples/Design-method-simulate-DesignGrouped.R new file mode 100644 index 000000000..f37b725f3 --- /dev/null +++ b/examples/Design-method-simulate-DesignGrouped.R @@ -0,0 +1,74 @@ +# Assemble ingredients for our group design. +my_stopping <- StoppingTargetProb(target = c(0.2, 0.35), prob = 0.5) | + StoppingMinPatients(20) | + StoppingMissingDose() +my_increments <- IncrementsDoseLevels(levels = 3L) +my_next_best <- NextBestNCRM( + target = c(0.2, 0.3), + overdose = c(0.3, 1), + max_overdose_prob = 0.3 +) +my_cohort_size <- CohortSizeConst(3) +empty_data <- Data(doseGrid = c(0.1, 0.5, 1.5, 3, 6, seq(from = 10, to = 80, by = 2))) +my_model <- LogisticLogNormalGrouped( + mean = c(-4, -4, -4, -4), + cov = diag(rep(6, 4)), + ref_dose = 0.1 +) + +# Put together the design. Note that if we only specify the mono arm, +# then the combo arm is having the same settings. +my_design <- DesignGrouped( + model = my_model, + mono = Design( + model = my_model, + stopping = my_stopping, + increments = my_increments, + nextBest = my_next_best, + cohort_size = my_cohort_size, + data = empty_data, + startingDose = 0.1 + ), + first_cohort_mono_only = TRUE, + same_dose = TRUE +) + +# Set up a realistic simulation scenario. +my_truth <- function(x) plogis(-4 + 0.2 * log(x / 0.1)) +my_combo_truth <- function(x) plogis(-4 + 0.5 * log(x / 0.1)) +matplot( + x = empty_data@doseGrid, + y = cbind( + mono = my_truth(empty_data@doseGrid), + combo = my_combo_truth(empty_data@doseGrid) + ), + type = "l", + ylab = "true DLT prob", + xlab = "dose" +) +legend("topright", c("mono", "combo"), lty = c(1, 2), col = c(1, 2)) + +# Start the simulations. +set.seed(123) +my_sims <- simulate( + my_design, + nsim = 4, # This should be at least 100 in actual applications. + seed = 123, + truth = my_truth, + combo_truth = my_combo_truth +) + +# Looking at the summary of the simulations: +mono_sims_sum <- summary(my_sims$mono, truth = my_truth) +combo_sims_sum <- summary(my_sims$combo, truth = my_combo_truth) + +mono_sims_sum +combo_sims_sum + +plot(mono_sims_sum) +plot(combo_sims_sum) + +# Looking at specific simulated trials: +trial_index <- 1 +plot(my_sims$mono@data[[trial_index]]) +plot(my_sims$combo@data[[trial_index]]) diff --git a/man/DesignGrouped-class.Rd b/man/DesignGrouped-class.Rd new file mode 100644 index 000000000..d44ef56eb --- /dev/null +++ b/man/DesignGrouped-class.Rd @@ -0,0 +1,124 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Design-class.R +\docType{class} +\name{DesignGrouped-class} +\alias{DesignGrouped-class} +\alias{.DesignGrouped} +\alias{DesignGrouped} +\alias{.DefaultDesignGrouped} +\title{\code{DesignGrouped}} +\usage{ +DesignGrouped( + model, + mono, + combo = mono, + first_cohort_mono_only = TRUE, + same_dose = TRUE, + ... +) +} +\arguments{ +\item{model}{(\code{LogisticLogNormalGrouped})\cr see slot definition.} + +\item{mono}{(\code{Design})\cr see slot definition.} + +\item{combo}{(\code{Design})\cr see slot definition.} + +\item{first_cohort_mono_only}{(\code{flag})\cr see slot definition.} + +\item{same_dose}{(\code{flag})\cr see slot definition.} + +\item{...}{not used.} +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} + +\code{\link{DesignGrouped}} combines two \code{\link{Design}} objects: one for the mono and one +for the combo arm of a joint dose escalation design. +} +\details{ +Note that the model slots inside the \code{mono} and \code{combo} parameters +are ignored (because we don't fit separate regression models for the mono and +combo arms). Instead, the \code{model} parameter is used to fit a joint regression +model for the mono and combo arms together. +} +\section{Slots}{ + +\describe{ +\item{\code{model}}{(\code{LogisticLogNormalGrouped})\cr the model to be used, currently only one +class is allowed.} + +\item{\code{mono}}{(\code{Design})\cr defines the dose escalation rules for the mono arm, see +details.} + +\item{\code{combo}}{(\code{Design})\cr defines the dose escalation rules for the combo arm, see +details.} + +\item{\code{first_cohort_mono_only}}{(\code{flag})\cr whether first test one mono agent cohort, and then +once its DLT data has been collected, we proceed from the second cohort onwards with +concurrent mono and combo cohorts.} + +\item{\code{same_dose}}{(\code{flag})\cr whether the lower dose of the separately determined mono and combo +doses should be used as the next dose for both mono and combo.} +}} + +\note{ +Typically, end-users will not use the \code{.DefaultDesignGrouped()} function. +} +\examples{ +empty_data <- Data(doseGrid = c(1, 3, 5, 10, 15, 20, 25, 40, 50, 80, 100)) + +# Initialize the joint model. +my_model <- LogisticLogNormalGrouped( + mean = c(-0.85, 0, 1, 0), + cov = diag(1, 4), + ref_dose = 56 +) + +# Choose the rule for selecting the next dose. +my_next_best <- NextBestNCRM( + target = c(0.2, 0.35), + overdose = c(0.35, 1), + max_overdose_prob = 0.25 +) + +# Choose the rule for the cohort-size. +my_size1 <- CohortSizeRange( + intervals = c(0, 30), + cohort_size = c(1, 3) +) +my_size2 <- CohortSizeDLT( + intervals = c(0, 1), + cohort_size = c(1, 3) +) +my_size <- maxSize(my_size1, my_size2) + +# Choose the rule for stopping. +my_stopping1 <- StoppingMinCohorts(nCohorts = 3) +my_stopping2 <- StoppingTargetProb( + target = c(0.2, 0.35), + prob = 0.5 +) +my_stopping3 <- StoppingMinPatients(nPatients = 20) +my_stopping <- (my_stopping1 & my_stopping2) | my_stopping3 + +# Choose the rule for dose increments. +my_increments <- IncrementsRelative( + intervals = c(0, 20), + increments = c(1, 0.33) +) + +# Initialize the design. +design <- DesignGrouped( + model = my_model, + mono = Design( + model = .DefaultModelLogNormal(), # Ignored. + nextBest = my_next_best, + stopping = my_stopping, + increments = my_increments, + cohort_size = my_size, + data = empty_data, + startingDose = 3 + ) +) +} diff --git a/man/get_result_list.Rd b/man/get_result_list.Rd index da89b4edc..e47267ea0 100644 --- a/man/get_result_list.Rd +++ b/man/get_result_list.Rd @@ -1,8 +1,8 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/Design-methods.R +% Please edit documentation in R/helpers_design.R \name{get_result_list} \alias{get_result_list} -\title{Helper function to obtain simulation results list} +\title{Helper Function to Obtain Simulation Results List} \usage{ get_result_list(fun, nsim, vars, parallel, n_cores) } @@ -24,8 +24,6 @@ The list with all simulation results (one iteration corresponds to one list element). } \description{ -\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#stable}{\figure{lifecycle-stable.svg}{options: alt='[Stable]'}}}{\strong{[Stable]}} - The function \code{fun} can use variables that are visible to itself. The names of these variables have to be given in the vector \code{vars}. } diff --git a/man/h_add_dlts.Rd b/man/h_add_dlts.Rd new file mode 100644 index 000000000..a7809024e --- /dev/null +++ b/man/h_add_dlts.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/helpers_design.R +\name{h_add_dlts} +\alias{h_add_dlts} +\title{Helper Function to Add Randomly Generated DLTs During Simulations} +\usage{ +h_add_dlts(data, dose, truth, cohort_size, first_separate) +} +\arguments{ +\item{data}{(\code{Data})\cr what data to start from.} + +\item{dose}{(\code{number})\cr current dose.} + +\item{truth}{(\code{function})\cr defines the true probability for a DLT at a dose.} + +\item{cohort_size}{(\code{CohortSize})\cr the cohort size rule to use.} + +\item{first_separate}{(\code{flag})\cr whether the first patient is enrolled separately.} +} +\value{ +The updated \code{data}. +} +\description{ +Helper Function to Add Randomly Generated DLTs During Simulations +} +\keyword{internal} diff --git a/man/nextBest.Rd b/man/nextBest.Rd index fd438dbfc..b9b730b3d 100644 --- a/man/nextBest.Rd +++ b/man/nextBest.Rd @@ -52,7 +52,7 @@ nextBest(nextBest, doselimit, samples, model, data, ...) \S4method{nextBest}{NextBestTD,numeric,missing,LogisticIndepBeta,Data}(nextBest, doselimit = Inf, model, data, in_sim = FALSE, ...) -\S4method{nextBest}{NextBestTDsamples,numeric,Samples,LogisticIndepBeta,Data}(nextBest, doselimit = Inf, samples, model, data, ...) +\S4method{nextBest}{NextBestTDsamples,numeric,Samples,LogisticIndepBeta,Data}(nextBest, doselimit = Inf, samples, model, data, in_sim, ...) \S4method{nextBest}{NextBestMaxGain,numeric,missing,ModelTox,DataDual}( nextBest, diff --git a/man/set_seed.Rd b/man/set_seed.Rd index f63be09de..982d23bad 100644 --- a/man/set_seed.Rd +++ b/man/set_seed.Rd @@ -1,8 +1,8 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/Design-methods.R +% Please edit documentation in R/helpers_design.R \name{set_seed} \alias{set_seed} -\title{Helper function to set and save the RNG seed} +\title{Helper Function to Set and Save the RNG Seed} \usage{ set_seed(seed = NULL) } @@ -24,4 +24,3 @@ the obtained simulation results. This code is basically copied from \code{stats:::simulate.lm}. } -\keyword{programming} diff --git a/man/simulate-DesignGrouped-method.Rd b/man/simulate-DesignGrouped-method.Rd new file mode 100644 index 000000000..54c2c4880 --- /dev/null +++ b/man/simulate-DesignGrouped-method.Rd @@ -0,0 +1,136 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Design-methods.R +\name{simulate,DesignGrouped-method} +\alias{simulate,DesignGrouped-method} +\alias{simulate-DesignGrouped} +\title{Simulate Method for the \code{\link{DesignGrouped}} Class} +\usage{ +\S4method{simulate}{DesignGrouped}( + object, + nsim = 1L, + seed = NULL, + truth, + combo_truth, + args = data.frame(), + firstSeparate = FALSE, + mcmcOptions = McmcOptions(), + parallel = FALSE, + nCores = min(parallelly::availableCores(), 5), + ... +) +} +\arguments{ +\item{object}{(\code{DesignGrouped})\cr the design we want to simulate trials from.} + +\item{nsim}{(\code{number})\cr how many trials should be simulated.} + +\item{seed}{(\code{RNGstate})\cr generated with \code{\link[=set_seed]{set_seed()}}.} + +\item{truth}{(\code{function})\cr a function which takes as input a dose (vector) and +returns the true probability (vector) for toxicity for the mono arm. +Additional arguments can be supplied in \code{args}.} + +\item{combo_truth}{(\code{function})\cr same as \code{truth} but for the combo arm.} + +\item{args}{(\code{data.frame})\cr optional \code{data.frame} with arguments that work +for both the \code{truth} and \code{combo_truth} functions. The column names correspond to +the argument names, the rows to the values of the arguments. The rows are +appropriately recycled in the \code{nsim} simulations.} + +\item{firstSeparate}{(\code{flag})\cr whether to enroll the first patient separately +from the rest of the cohort and close the cohort in case a DLT occurs in this +first patient.} + +\item{mcmcOptions}{(\code{McmcOptions})\cr MCMC options for each evaluation in the trial.} + +\item{parallel}{(\code{flag})\cr whether the simulation runs are parallelized across the +cores of the computer.} + +\item{nCores}{(\code{number})\cr how many cores should be used for parallel computing.} + +\item{...}{not used.} +} +\value{ +A list of \code{mono} and \code{combo} simulation results as \code{\link{Simulations}} objects. +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} + +A simulate method for \code{\link{DesignGrouped}} designs. +} +\examples{ +# Assemble ingredients for our group design. +my_stopping <- StoppingTargetProb(target = c(0.2, 0.35), prob = 0.5) | + StoppingMinPatients(20) | + StoppingMissingDose() +my_increments <- IncrementsDoseLevels(levels = 3L) +my_next_best <- NextBestNCRM( + target = c(0.2, 0.3), + overdose = c(0.3, 1), + max_overdose_prob = 0.3 +) +my_cohort_size <- CohortSizeConst(3) +empty_data <- Data(doseGrid = c(0.1, 0.5, 1.5, 3, 6, seq(from = 10, to = 80, by = 2))) +my_model <- LogisticLogNormalGrouped( + mean = c(-4, -4, -4, -4), + cov = diag(rep(6, 4)), + ref_dose = 0.1 +) + +# Put together the design. Note that if we only specify the mono arm, +# then the combo arm is having the same settings. +my_design <- DesignGrouped( + model = my_model, + mono = Design( + model = my_model, + stopping = my_stopping, + increments = my_increments, + nextBest = my_next_best, + cohort_size = my_cohort_size, + data = empty_data, + startingDose = 0.1 + ), + first_cohort_mono_only = TRUE, + same_dose = TRUE +) + +# Set up a realistic simulation scenario. +my_truth <- function(x) plogis(-4 + 0.2 * log(x / 0.1)) +my_combo_truth <- function(x) plogis(-4 + 0.5 * log(x / 0.1)) +matplot( + x = empty_data@doseGrid, + y = cbind( + mono = my_truth(empty_data@doseGrid), + combo = my_combo_truth(empty_data@doseGrid) + ), + type = "l", + ylab = "true DLT prob", + xlab = "dose" +) +legend("topright", c("mono", "combo"), lty = c(1, 2), col = c(1, 2)) + +# Start the simulations. +set.seed(123) +my_sims <- simulate( + my_design, + nsim = 4, # This should be at least 100 in actual applications. + seed = 123, + truth = my_truth, + combo_truth = my_combo_truth +) + +# Looking at the summary of the simulations: +mono_sims_sum <- summary(my_sims$mono, truth = my_truth) +combo_sims_sum <- summary(my_sims$combo, truth = my_combo_truth) + +mono_sims_sum +combo_sims_sum + +plot(mono_sims_sum) +plot(combo_sims_sum) + +# Looking at specific simulated trials: +trial_index <- 1 +plot(my_sims$mono@data[[trial_index]]) +plot(my_sims$combo@data[[trial_index]]) +} diff --git a/man/v_design.Rd b/man/v_design.Rd index c77e39a5c..8f00076d7 100644 --- a/man/v_design.Rd +++ b/man/v_design.Rd @@ -3,9 +3,12 @@ \name{v_design} \alias{v_design} \alias{v_rule_design} +\alias{v_design_grouped} \title{Internal Helper Functions for Validation of \code{\link{RuleDesign}} Objects} \usage{ v_rule_design(object) + +v_design_grouped(object) } \arguments{ \item{object}{(\code{RuleDesign})\cr object to validate.} @@ -25,4 +28,7 @@ These functions are only used internally to validate the format of an input \item \code{v_rule_design()}: validates that the \code{\link{RuleDesign}} object contains valid \code{startingDose}. +\item \code{v_design_grouped()}: validates that the \code{\link{DesignGrouped}} object +contains valid flags. + }} diff --git a/tests/testthat/helper-model.R b/tests/testthat/helper-model.R index e1725a900..b56e23665 100644 --- a/tests/testthat/helper-model.R +++ b/tests/testthat/helper-model.R @@ -386,22 +386,3 @@ h_get_fractional_crm <- function() { sigma2 = 2 ) } - -.NeedsExtraProbModel <- setClass("NeedsExtraProbModel", contains = "LogisticKadane") -setMethod( - f = "prob", - signature = signature( - dose = "numeric", - model = "NeedsExtraProbModel", - samples = "Samples" - ), - definition = function(dose, model, samples, extra_argument) { - if (missing(extra_argument)) stop("we need extra_argument") - # We don't forward to LogisticKadane method here since that would - # not accept the extra argument. - rep(0.5, size(samples)) - } -) -h_needs_extra_prob_model <- function() { - .NeedsExtraProbModel(h_get_logistic_kadane()) -} diff --git a/tests/testthat/test-Design-class.R b/tests/testthat/test-Design-class.R index adf469a4b..ea21ead52 100644 --- a/tests/testthat/test-Design-class.R +++ b/tests/testthat/test-Design-class.R @@ -424,3 +424,47 @@ test_that("DADesign user constructor arguments names are as expected", { ordered = TRUE ) }) + +# DesignGrouped ---- + +test_that(".DesignGrouped works as expected", { + result <- .DesignGrouped() + + expect_true(inherits(result, "CrmPackClass")) + expect_valid(result, "DesignGrouped") +}) + +test_that("DesignGrouped works as expected", { + empty_data <- Data(doseGrid = 2:50) + model <- .DefaultLogisticLogNormalGrouped() + stopping <- h_stopping_target_prob() + increments <- h_increments_relative() + placebo_cohort_size <- CohortSizeConst(0L) + next_best <- h_next_best_ncrm() + cohort_size <- CohortSizeRange(intervals = c(0, 30), cohort_size = c(1, 3)) + + result <- expect_silent( + DesignGrouped( + model = model, + mono = Design( + model, + stopping, + increments, + nextBest = next_best, + cohort_size = cohort_size, + data = empty_data, + startingDose = 3 + ) + ) + ) + + expect_valid(result, "DesignGrouped") + expect_identical(result@mono, result@combo) + expect_true(result@first_cohort_mono_only) + expect_true(result@same_dose) +}) + +test_that(".DefaultDesignGrouped works as expected", { + result <- .DefaultDesignGrouped() + expect_valid(result, "DesignGrouped") +}) diff --git a/tests/testthat/test-Design-methods.R b/tests/testthat/test-Design-methods.R index 773875b4f..f5a7c01eb 100644 --- a/tests/testthat/test-Design-methods.R +++ b/tests/testthat/test-Design-methods.R @@ -1,51 +1,3 @@ -# helper functions ---- - -## set_seed ---- - -test_that("set_seed returns correct value if seed is a value", { - seed <- 1.909 - seed_int <- 1 - - RNGkind("default") - rng_state <- set_seed(seed) - attr(seed_int, "kind") <- list("Mersenne-Twister", "Inversion", "Rejection") - expect_equal(rng_state, seed_int) - - RNGkind("Super-Duper") - rng_state <- set_seed(seed) - attr(seed_int, "kind") <- list("Super-Duper", "Inversion", "Rejection") - expect_equal(rng_state, seed_int) - - RNGkind("default") -}) - -test_that("set_seed returns correct value if seed is NULL", { - seed <- NULL - - RNGkind("default") - rng_state <- set_seed(seed) - expect_equal(rng_state, .Random.seed) - - RNGkind("Super-Duper") - rng_state <- set_seed(seed) - expect_equal(rng_state, .Random.seed) - - RNGkind("default") -}) - -## get_result_list ---- - -test_that("get_result_list returns correct value", { - res <- get_result_list(mean, 2, NULL, FALSE, 5) - expect_equal(res, list(1, 2)) - - res <- get_result_list(length, 2, NULL, FALSE, 5) - expect_equal(res, list(1, 1)) - - expect_error(get_result_list(length, 2, NULL, 5, 5)) - expect_error(get_result_list(length, 2, NULL, FALSE, 0)) -}) - # simulate ---- ## NextBestInfTheory ---- @@ -182,6 +134,151 @@ test_that("stop_reasons can be NA with certain stopping rule settings", { expect_identical(result, expected) }) +## DesignGrouped ---- + +test_that("simulate for DesignGrouped works as expected", { + object <- DesignGrouped( + model = LogisticLogNormalGrouped(mean = rep(-1, 4), cov = diag(5, 4), ref_dose = 1), + mono = Design( + model = .LogisticNormal(), + nextBest = NextBestNCRM(target = c(0.3, 0.6), overdose = c(0.6, 1), max_overdose_prob = 0.7), + stopping = StoppingMinPatients(nPatients = 9), + increments = IncrementsDoseLevels(levels = 5), + cohort_size = CohortSizeConst(3), + data = Data(doseGrid = 1:100), + startingDose = 1 + ), + same_dose = TRUE, + first_cohort_mono_only = TRUE + ) + my_truth <- function(x) plogis(-4 + 0.2 * log(x / 0.1)) + my_combo_truth <- function(x) plogis(-4 + 0.5 * log(x / 0.1)) + + result <- expect_silent(simulate( + object, + nsim = 2, + seed = 123, + truth = my_truth, + combo_truth = my_combo_truth, + mcmcOptions = h_get_mcmc_options() + )) + + expect_list(result) + expect_names(names(result), identical.to = c("mono", "combo")) + expect_valid(result$mono, "Simulations") + expect_valid(result$combo, "Simulations") + + mono_trial <- result$mono@data[[2L]] + combo_trial <- result$combo@data[[2L]] + + # First cohort is only mono at the starting dose (lowest in dose grid). + expect_true(all(mono_trial@xLevel[1:3] == 1)) + + # We have the same dose for subsequent cohorts. + expect_true(all(mono_trial@xLevel[4:6] == combo_trial@xLevel[1:3])) + expect_true(all(mono_trial@xLevel[7:9] == combo_trial@xLevel[4:6])) +}) + +test_that("simulate for DesignGrouped works as expected with different doses, parallel first cohort", { + object <- DesignGrouped( + model = LogisticLogNormalGrouped(mean = rep(-1, 4), cov = diag(5, 4), ref_dose = 1), + mono = Design( + model = .LogisticNormal(), + nextBest = NextBestNCRM(target = c(0.3, 0.6), overdose = c(0.6, 1), max_overdose_prob = 0.7), + stopping = StoppingMinPatients(nPatients = 20), + increments = IncrementsDoseLevels(levels = 5), + cohort_size = CohortSizeConst(3), + data = Data(doseGrid = 1:100), + startingDose = 1 + ), + same_dose = FALSE, + first_cohort_mono_only = FALSE + ) + my_truth <- function(x) plogis(-4 + 0.2 * log(x / 0.1)) + my_combo_truth <- function(x) plogis(-4 + 0.9 * log(x / 0.1)) + + result <- expect_silent(simulate( + object, + nsim = 1, + seed = 123, + truth = my_truth, + combo_truth = my_combo_truth, + mcmcOptions = h_get_mcmc_options() + )) + + expect_list(result) + expect_names(names(result), identical.to = c("mono", "combo")) + expect_valid(result$mono, "Simulations") + expect_valid(result$combo, "Simulations") + + mono_trial <- result$mono@data[[1L]] + combo_trial <- result$combo@data[[1L]] + + # First cohort is joint at starting dose. + expect_true(all(mono_trial@xLevel[1:3] == combo_trial@xLevel[1:3])) + + # We have different doses in subsequent cohorts. + expect_false(all(mono_trial@xLevel[4:20] == combo_trial@xLevel[4:20])) +}) + +test_that("simulate for DesignGrouped works when first patient is dosed separately, different combo design", { + object <- DesignGrouped( + model = LogisticLogNormalGrouped(mean = rep(-1, 4), cov = diag(5, 4), ref_dose = 1), + mono = Design( + model = .LogisticNormal(), + nextBest = NextBestNCRM(target = c(0.3, 0.6), overdose = c(0.6, 1), max_overdose_prob = 0.7), + stopping = StoppingMinPatients(nPatients = 10), + increments = IncrementsDoseLevels(levels = 3), + cohort_size = CohortSizeConst(2), + data = Data(doseGrid = 1:100), + startingDose = 10 + ), + combo = Design( + model = .LogisticNormal(), + nextBest = NextBestNCRM(target = c(0.3, 0.6), overdose = c(0.6, 1), max_overdose_prob = 0.7), + stopping = StoppingMinPatients(nPatients = 20), + increments = IncrementsDoseLevels(levels = 5), + cohort_size = CohortSizeConst(3), + data = Data(doseGrid = 1:100), + startingDose = 1 + ), + same_dose = FALSE, + first_cohort_mono_only = FALSE + ) + my_truth <- function(x) plogis(-4 + 0.2 * log(x / 0.1)) + my_combo_truth <- function(x) plogis(-2 + 0.9 * log(x / 0.1)) + + result <- expect_silent(simulate( + object, + nsim = 1, + seed = 123, + truth = my_truth, + firstSeparate = TRUE, + combo_truth = my_combo_truth, + mcmcOptions = h_get_mcmc_options() + )) + + expect_list(result) + expect_names(names(result), identical.to = c("mono", "combo")) + expect_valid(result$mono, "Simulations") + expect_valid(result$combo, "Simulations") + + mono_trial <- result$mono@data[[1L]] + combo_trial <- result$combo@data[[1L]] + + # We expect at least one cohort with just one patient in the combo arm here + # because of the high toxicity. + expect_true(any(table(combo_trial@cohort) == 1)) + + # Check that we had the different cohort sizes between the two arms. + expect_true(max(table(mono_trial@cohort)) == 2) + expect_true(max(table(combo_trial@cohort)) == 3) + + # Check that we had different starting doses in the two arms. + expect_true(mono_trial@x[1] == 10) + expect_true(combo_trial@x[1] == 1) +}) + # examine ---- ## DADesign ---- diff --git a/tests/testthat/test-Design-validity.R b/tests/testthat/test-Design-validity.R index 5374ad254..651a1e3ba 100644 --- a/tests/testthat/test-Design-validity.R +++ b/tests/testthat/test-Design-validity.R @@ -61,3 +61,24 @@ test_that("v_rule_design returns message when startingDose is not on doseGrid", object@startingDose <- 6.5 expect_equal(v_rule_design(object), err_msg) }) + +## v_design_grouped ---- + +test_that("v_design_grouped passes for valid object", { + object <- .DesignGrouped() + expect_true(v_design_grouped(object)) +}) + +test_that("v_design_grouped messages wrong flag slots as expected", { + object <- .DesignGrouped() + + object@same_dose <- c(NA, TRUE) + object@first_cohort_mono_only <- logical() + + result <- v_design_grouped(object) + expected <- c( + "first_cohort_mono_only must be a flag", + "same_dose must be a flag" + ) + expect_identical(result, expected) +}) diff --git a/tests/testthat/test-Rules-methods.R b/tests/testthat/test-Rules-methods.R index 890dd0f06..57d504c4c 100644 --- a/tests/testthat/test-Rules-methods.R +++ b/tests/testthat/test-Rules-methods.R @@ -81,13 +81,13 @@ test_that("nextBest-NextBestNCRM returns expected values of the objects (no dose }) test_that("nextBest-NextBestNCRM can accept additional arguments and pass them to prob inside", { - my_data <- h_get_data() - my_model <- h_needs_extra_prob_model() + my_data <- h_get_data_grouped() + my_model <- .DefaultLogisticLogNormalGrouped() my_samples <- mcmc(my_data, my_model, h_get_mcmc_options(samples = 10, burnin = 10)) nb_ncrm <- NextBestNCRM( target = c(0.2, 0.35), overdose = c(0.35, 1), max_overdose_prob = 0.25 ) - result <- nextBest(nb_ncrm, Inf, my_samples, my_model, my_data, extra_argument = "foo") + result <- nextBest(nb_ncrm, Inf, my_samples, my_model, my_data, group = "mono") expect_identical(result$value, NA_real_) }) @@ -1857,8 +1857,8 @@ test_that("StoppingTargetProb works correctly when above threshold", { }) test_that("stopTrial-StoppingTargetProb can accept additional arguments and pass them to prob", { - my_data <- h_get_data() - my_model <- h_needs_extra_prob_model() + my_data <- h_get_data_grouped() + my_model <- .DefaultLogisticLogNormalGrouped() my_samples <- mcmc(my_data, my_model, h_get_mcmc_options(samples = 10, burnin = 10)) stopping <- StoppingTargetProb(target = c(0.1, 0.4), prob = 0.3) result <- stopTrial( @@ -1867,7 +1867,7 @@ test_that("stopTrial-StoppingTargetProb can accept additional arguments and pass samples = my_samples, model = my_model, data = my_data, - extra_argument = "bla" + group = "combo" ) expect_false(result) }) diff --git a/tests/testthat/test-helpers.R b/tests/testthat/test-helpers.R index ad687b3ba..90e0be660 100644 --- a/tests/testthat/test-helpers.R +++ b/tests/testthat/test-helpers.R @@ -473,7 +473,7 @@ test_that("h_find_interval works as expected for custom replacement", { test_that("default constructors exist for all subclasses of GeneralModel", { allModelSubclasses <- names(getClassDef("GeneralModel")@subclasses) # Exceptions. - classesNotToTest <- c("DualEndpoint", "NeedsExtraProbModel") + classesNotToTest <- "DualEndpoint" classesToTest <- setdiff(allModelSubclasses, classesNotToTest) lapply( classesToTest, diff --git a/tests/testthat/test-helpers_design.R b/tests/testthat/test-helpers_design.R new file mode 100644 index 000000000..2c8c18641 --- /dev/null +++ b/tests/testthat/test-helpers_design.R @@ -0,0 +1,100 @@ +# set_seed ---- + +test_that("set_seed returns correct value if seed is a value", { + seed <- 1.909 + seed_int <- 1 + + RNGkind("default") + rng_state <- set_seed(seed) + attr(seed_int, "kind") <- list("Mersenne-Twister", "Inversion", "Rejection") + expect_equal(rng_state, seed_int) + + RNGkind("Super-Duper") + rng_state <- set_seed(seed) + attr(seed_int, "kind") <- list("Super-Duper", "Inversion", "Rejection") + expect_equal(rng_state, seed_int) + + RNGkind("default") +}) + +test_that("set_seed returns correct value if seed is NULL", { + seed <- NULL + + RNGkind("default") + rng_state <- set_seed(seed) + expect_equal(rng_state, .Random.seed) + + RNGkind("Super-Duper") + rng_state <- set_seed(seed) + expect_equal(rng_state, .Random.seed) + + RNGkind("default") +}) + +# get_result_list ---- + +test_that("get_result_list returns correct value", { + res <- get_result_list(mean, 2, NULL, FALSE, 5) + expect_equal(res, list(1, 2)) + + res <- get_result_list(length, 2, NULL, FALSE, 5) + expect_equal(res, list(1, 1)) + + expect_error(get_result_list(length, 2, NULL, 5, 5)) + expect_error(get_result_list(length, 2, NULL, FALSE, 0)) +}) + +# h_add_dlts ---- + +test_that("h_add_dlts works as expected", { + data <- h_get_data() + cohort_size <- CohortSizeConst(3) + + set.seed(123) + result <- expect_silent(h_add_dlts( + data = data, + dose = data@doseGrid[3], + truth = plogis, + cohort_size = cohort_size, + first_separate = FALSE + )) + expect_valid(result, "Data") + expect_equal(tail(result@x, 3), rep(data@doseGrid[3], 3)) + expect_true(data@nObs + 3 == result@nObs) +}) + +test_that("h_add_dlts works as expected when first separate patient has a DLT", { + data <- h_get_data() + cohort_size <- CohortSizeConst(3) + + set.seed(123) + result <- expect_silent(h_add_dlts( + data = data, + dose = data@doseGrid[3], + truth = function(dose, ...) 1, # Make sure the first patient has a DLT. + cohort_size = cohort_size, + first_separate = TRUE + )) + expect_valid(result, "Data") + expect_true(tail(result@y, 1) == 1) + expect_equal(tail(result@x, 1), data@doseGrid[3]) + expect_true(data@nObs + 1 == result@nObs) +}) + +test_that("h_add_dlts works as expected when first separate patient does not have a DLT", { + data <- h_get_data() + cohort_size <- CohortSizeConst(3) + + set.seed(123) + result <- expect_silent(h_add_dlts( + data = data, + dose = data@doseGrid[3], + truth = function(dose, ...) 0, # Make sure the first patient does not have a DLT. + cohort_size = cohort_size, + first_separate = TRUE + )) + expect_valid(result, "Data") + expect_equal(tail(result@y, 3), rep(0, 3)) + expect_equal(tail(result@x, 3), rep(data@doseGrid[3], 3)) + expect_true(data@nObs + 3 == result@nObs) +})