From 7270a9ac78212339d1200f9c99a69314cb3fb640 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Tue, 12 Sep 2023 15:25:10 +0200 Subject: [PATCH 01/31] add class with constructors --- NAMESPACE | 4 + R/Model-class.R | 91 +++++++++++++++++++ .../Model-class-LogisticLogNormalGrouped.R | 6 ++ man/LogisticLogNormalGrouped-class.Rd | 55 +++++++++++ 4 files changed, 156 insertions(+) create mode 100644 examples/Model-class-LogisticLogNormalGrouped.R create mode 100644 man/LogisticLogNormalGrouped-class.Rd diff --git a/NAMESPACE b/NAMESPACE index 50bcac1df..5fc20e077 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -43,6 +43,7 @@ export(.DefaultIncrementsRelativeParts) export(.DefaultLogisticKadane) export(.DefaultLogisticKadaneBetaGamma) export(.DefaultLogisticLogNormal) +export(.DefaultLogisticLogNormalGrouped) export(.DefaultLogisticLogNormalMixture) export(.DefaultLogisticLogNormalSub) export(.DefaultLogisticNormal) @@ -112,6 +113,7 @@ export(.LogisticIndepBeta) export(.LogisticKadane) export(.LogisticKadaneBetaGamma) export(.LogisticLogNormal) +export(.LogisticLogNormalGrouped) export(.LogisticLogNormalMixture) export(.LogisticLogNormalSub) export(.LogisticNormal) @@ -210,6 +212,7 @@ export(LogisticIndepBeta) export(LogisticKadane) export(LogisticKadaneBetaGamma) export(LogisticLogNormal) +export(LogisticLogNormalGrouped) export(LogisticLogNormalMixture) export(LogisticLogNormalSub) export(LogisticNormal) @@ -400,6 +403,7 @@ exportClasses(LogisticIndepBeta) exportClasses(LogisticKadane) exportClasses(LogisticKadaneBetaGamma) exportClasses(LogisticLogNormal) +exportClasses(LogisticLogNormalGrouped) exportClasses(LogisticLogNormalMixture) exportClasses(LogisticLogNormalSub) exportClasses(LogisticNormal) diff --git a/R/Model-class.R b/R/Model-class.R index 44a53d052..5b08aeb15 100644 --- a/R/Model-class.R +++ b/R/Model-class.R @@ -516,6 +516,97 @@ ProbitLogNormalRel <- function(mean, cov, ref_dose = 1) { ProbitLogNormalRel(mean = c(-0.85, 1), cov = matrix(c(1, -0.5, -0.5, 1), nrow = 2)) } +# LogisticLogNormalGrouped ---- + +## class ---- + +#' `LogisticLogNormalGrouped` +#' +#' @description `r lifecycle::badge("experimental")` +#' +#' [`LogisticLogNormalGrouped`] is the class for a logistic regression model +#' for both the mono and the combo arms of the simultaneous dose escalation +#' design. +#' +#' @details The continuous covariate is the natural logarithm of the dose \eqn{x} divided by +#' the reference dose \eqn{x*} as in [`LogisticLogNormal`]. In addition, +#' \eqn{I_c} is a binary indicator covariate which is 1 for the combo arm and 0 for the mono arm. +#' The model is then defined as: +#' \deqn{logit[p(x)] = (alpha0 + I_c * delta0) + (alpha1 + I_c * delta1) * log(x / x*),} +#' where \eqn{p(x)} is the probability of observing a DLT for a given dose \eqn{x}, +#' and `delta0` and `delta1` are the differences in the combo arm compared to the mono intercept +#' and slope parameters `alpha0` and `alpha1`. +#' The prior is defined as \deqn{(alpha0, log(delta0), log(alpha1), log(delta1)) ~ Normal(mean, cov).} +#' +#' @seealso [`ModelLogNormal`], [`LogisticLogNormal`]. +#' +#' @aliases LogisticLogNormalGrouped +#' @export +#' +.LogisticLogNormalGrouped <- setClass( + Class = "LogisticLogNormalGrouped", + contains = "ModelLogNormal" +) + +## constructor ---- + +#' @rdname LogisticLogNormalGrouped-class +#' +#' @inheritParams ModelLogNormal +#' +#' @export +#' @example examples/Model-class-LogisticLogNormalGrouped.R +#' +LogisticLogNormalGrouped <- function(mean, cov, ref_dose = 1) { + params <- ModelParamsNormal(mean, cov) + .LogisticLogNormalGrouped( + params = params, + ref_dose = positive_number(ref_dose), + priormodel = function() { + theta ~ dmnorm(mean, prec) + alpha0 <- theta[1] + delta0 <- exp(theta[2]) + alpha1 <- exp(theta[3]) + delta1 <- exp(theta[4]) + }, + datamodel = function() { + for (i in 1:nObs) { + logit(p[i]) <- (alpha0 + is_combo[i] * delta0) + + (alpha1 + is_combo[i] * delta1) * log(x[i] / ref_dose) + y[i] ~ dbern(p[i]) + } + }, + modelspecs = function(group, from_prior) { + ms <- list( + mean = params@mean, + prec = params@prec + ) + if (!from_prior) { + ms$ref_dose <- ref_dose + ms$is_combo <- as.integer(group == "combo") + } + ms + }, + init = function() { + list(theta = c(0, 1, 1, 1)) + }, + datanames = c("nObs", "y", "x"), + sample = c("alpha0", "delta0", "alpha1", "delta1") + ) +} + +## default constructor ---- + +#' @rdname LogisticLogNormalGrouped-class +#' @note Typically, end users will not use the `.DefaultLogisticLogNormalGrouped()` function. +#' @export +.DefaultLogisticLogNormalGrouped <- function() { + LogisticLogNormalGrouped( + mean = rep(0, 4), + cov = diag(rep(1, 4)), + ) +} + # LogisticKadane ---- ## class ---- diff --git a/examples/Model-class-LogisticLogNormalGrouped.R b/examples/Model-class-LogisticLogNormalGrouped.R new file mode 100644 index 000000000..acae7d244 --- /dev/null +++ b/examples/Model-class-LogisticLogNormalGrouped.R @@ -0,0 +1,6 @@ +my_model <- LogisticLogNormalGrouped( + mean = c(-0.85, 0, 1, 0), + cov = diag(1, 4), + ref_dose = 50 +) +my_model diff --git a/man/LogisticLogNormalGrouped-class.Rd b/man/LogisticLogNormalGrouped-class.Rd new file mode 100644 index 000000000..550b0d920 --- /dev/null +++ b/man/LogisticLogNormalGrouped-class.Rd @@ -0,0 +1,55 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Model-class.R +\docType{class} +\name{LogisticLogNormalGrouped-class} +\alias{LogisticLogNormalGrouped-class} +\alias{.LogisticLogNormalGrouped} +\alias{LogisticLogNormalGrouped} +\alias{.DefaultLogisticLogNormalGrouped} +\title{\code{LogisticLogNormalGrouped}} +\usage{ +LogisticLogNormalGrouped(mean, cov, ref_dose = 1) + +.DefaultLogisticLogNormalGrouped() +} +\arguments{ +\item{mean}{(\code{numeric})\cr the prior mean vector.} + +\item{cov}{(\code{matrix})\cr the prior covariance matrix. The precision matrix +\code{prec} is internally calculated as an inverse of \code{cov}.} + +\item{ref_dose}{(\code{number})\cr the reference dose \eqn{x*} (strictly positive +number).} +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}} + +\code{\link{LogisticLogNormalGrouped}} is the class for a logistic regression model +for both the mono and the combo arms of the simultaneous dose escalation +design. +} +\details{ +The continuous covariate is the natural logarithm of the dose \eqn{x} divided by +the reference dose \eqn{x*} as in \code{\link{LogisticLogNormal}}. In addition, +\eqn{I_c} is a binary indicator covariate which is 1 for the combo arm and 0 for the mono arm. +The model is then defined as: +\deqn{logit[p(x)] = (alpha0 + I_c * delta0) + (alpha1 + I_c * delta1) * log(x / x*),} +where \eqn{p(x)} is the probability of observing a DLT for a given dose \eqn{x}, +and \code{delta0} and \code{delta1} are the differences in the combo arm compared to the mono intercept +and slope parameters \code{alpha0} and \code{alpha1}. +The prior is defined as \deqn{(alpha0, log(delta0), log(alpha1), log(delta1)) ~ Normal(mean, cov).} +} +\note{ +Typically, end users will not use the \code{.DefaultLogisticLogNormalGrouped()} function. +} +\examples{ +my_model <- LogisticLogNormalGrouped( + mean = c(-0.85, 0, 1, 0), + cov = diag(1, 4), + ref_dose = 50 +) +my_model +} +\seealso{ +\code{\link{ModelLogNormal}}, \code{\link{LogisticLogNormal}}. +} From 74ab10125d651db3a094626ddd7800f1587790c6 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Tue, 12 Sep 2023 15:55:31 +0200 Subject: [PATCH 02/31] add tests for model class --- tests/testthat/_snaps/Model-class.md | 36 ++++++++++++++++++++++++ tests/testthat/test-Model-class.R | 42 ++++++++++++++++++++++++++++ 2 files changed, 78 insertions(+) diff --git a/tests/testthat/_snaps/Model-class.md b/tests/testthat/_snaps/Model-class.md index d50f57c6d..8f0ab4dfc 100644 --- a/tests/testthat/_snaps/Model-class.md +++ b/tests/testthat/_snaps/Model-class.md @@ -118,6 +118,42 @@ [1] 2.6848927 0.5973656 0.5178854 4.1900261 +# MCMC computes correct values for LogisticLogNormalGrouped model + + Code + result@data + Output + $alpha0 + [1] -1.848124 -1.848124 -2.195992 -2.195992 + + $alpha1 + [1] 0.08165851 0.08165851 0.31496464 0.31496464 + + $delta0 + [1] 0.1263061 0.1263061 0.1112928 0.1112928 + + $delta1 + [1] 0.06280840 0.06280840 0.09082656 0.09082656 + + +# MCMC computes correct values for LogisticLogNormalGrouped model and empty data + + Code + result@data + Output + $alpha0 + [1] -0.7258644 -2.3268647 1.2196258 0.9135351 + + $alpha1 + [1] 1.3236188 0.2700508 1.8853673 0.2820018 + + $delta0 + [1] 3.4298468 2.0426447 0.7105671 0.2355593 + + $delta1 + [1] 0.3087583 0.2956680 3.1296169 1.7350122 + + # MCMC computes correct values for LogisticKadane model Code diff --git a/tests/testthat/test-Model-class.R b/tests/testthat/test-Model-class.R index c47bec12a..ddba621bf 100644 --- a/tests/testthat/test-Model-class.R +++ b/tests/testthat/test-Model-class.R @@ -279,6 +279,48 @@ test_that("MCMC computes correct values for ProbitLogNormalRel model and empty d expect_snapshot(result@data) }) +# LogisticLogNormalGrouped ---- + +## constructor ---- + +test_that("LogisticLogNormalGrouped object can be created with user constructor", { + result <- expect_silent( + LogisticLogNormalGrouped( + mean = 1:4, + cov = diag(1:4, 4), + ref_dose = 2 + ) + ) + expect_valid(result, "LogisticLogNormalGrouped") +}) + +test_that(".DefaultLogisticLogNormalGrouped works as expected", { + expect_valid( + .DefaultLogisticLogNormalGrouped(), + "LogisticLogNormalGrouped" + ) +}) + +## mcmc ---- + +test_that("MCMC computes correct values for LogisticLogNormalGrouped model", { + data <- h_get_data_grouped() + model <- .DefaultLogisticLogNormalGrouped() + options <- h_get_mcmc_options() + + result <- mcmc(data = data, model = model, options = options) + expect_snapshot(result@data) +}) + +test_that("MCMC computes correct values for LogisticLogNormalGrouped model and empty data", { + data <- h_get_data_grouped(empty = TRUE) + model <- .DefaultLogisticLogNormalGrouped() + options <- h_get_mcmc_options() + + result <- mcmc(data = data, model = model, options = options) + expect_snapshot(result@data) +}) + # LogisticKadane ---- ## constructor ---- From e297b7742ffd007ad3181f1a6dfefea37abe6b4a Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Tue, 12 Sep 2023 20:23:31 +0200 Subject: [PATCH 03/31] add prob method --- R/Model-methods.R | 35 ++++++++++++++++++++++++++++ man/prob.Rd | 11 +++++++++ tests/testthat/test-Model-methods.R | 36 +++++++++++++++++++++++++++++ 3 files changed, 82 insertions(+) diff --git a/R/Model-methods.R b/R/Model-methods.R index 4a6bfa975..ae0259d2f 100644 --- a/R/Model-methods.R +++ b/R/Model-methods.R @@ -923,6 +923,41 @@ setMethod( } ) +## LogisticLogNormalGrouped ---- + +#' @describeIn prob method for [`LogisticLogNormalGrouped`] which needs `group` +#' argument in addition. +#' @param group (`character` or `factor`)\cr for [`LogisticLogNormalGrouped`], +#' indicating whether to calculate the probability for the `mono` or for +#' the `combo` arm. +#' @aliases prob-LogisticLogNormalGrouped +#' @export +#' +setMethod( + f = "prob", + signature = signature( + dose = "numeric", + model = "LogisticLogNormalGrouped", + samples = "Samples" + ), + definition = function(dose, model, samples, group) { + assert_numeric(dose, lower = 0L, any.missing = FALSE, min.len = 1L) + assert_subset(c("alpha0", "delta0", "alpha1", "delta1"), names(samples)) + assert_length(dose, len = size(samples)) + assert_multi_class(group, c("character", "factor")) + assert_subset(as.character(group), choices = c("mono", "combo")) + assert_length(group, len = size(samples)) + + alpha0 <- samples@data$alpha0 + delta0 <- samples@data$delta0 + alpha1 <- samples@data$alpha1 + delta1 <- samples@data$delta1 + ref_dose <- as.numeric(model@ref_dose) + is_combo <- as.integer(group == "combo") + plogis((alpha0 + is_combo * delta0) + (alpha1 + is_combo * delta1) * log(dose / ref_dose)) + } +) + ## LogisticKadane ---- #' @describeIn prob diff --git a/man/prob.Rd b/man/prob.Rd index 82684fcca..48450c936 100644 --- a/man/prob.Rd +++ b/man/prob.Rd @@ -12,6 +12,8 @@ \alias{prob-ProbitLogNormal} \alias{prob,numeric,ProbitLogNormalRel,Samples-method} \alias{prob-ProbitLogNormalRel} +\alias{prob,numeric,LogisticLogNormalGrouped,Samples-method} +\alias{prob-LogisticLogNormalGrouped} \alias{prob,numeric,LogisticKadane,Samples-method} \alias{prob-LogisticKadane} \alias{prob,numeric,LogisticKadaneBetaGamma,Samples-method} @@ -46,6 +48,8 @@ prob(dose, model, samples, ...) \S4method{prob}{numeric,ProbitLogNormalRel,Samples}(dose, model, samples) +\S4method{prob}{numeric,LogisticLogNormalGrouped,Samples}(dose, model, samples, group) + \S4method{prob}{numeric,LogisticKadane,Samples}(dose, model, samples) \S4method{prob}{numeric,LogisticKadaneBetaGamma,Samples}(dose, model, samples) @@ -80,6 +84,10 @@ dose escalation or pseudo DLE (dose-limiting events)/toxicity model.} used to compute toxicity probabilities. Can also be missing for some models.} \item{...}{model specific parameters when \code{samples} are not used.} + +\item{group}{(\code{character} or \code{factor})\cr for \code{\link{LogisticLogNormalGrouped}}, +indicating whether to calculate the probability for the \code{mono} or for +the \code{combo} arm.} } \value{ A \code{proportion} or \code{numeric} vector with the toxicity probabilities. @@ -115,6 +123,9 @@ correspond to the sampling index, i.e. the layout is then \item \code{prob(dose = numeric, model = ProbitLogNormalRel, samples = Samples)}: +\item \code{prob(dose = numeric, model = LogisticLogNormalGrouped, samples = Samples)}: method for \code{\link{LogisticLogNormalGrouped}} which needs \code{group} +argument in addition. + \item \code{prob(dose = numeric, model = LogisticKadane, samples = Samples)}: \item \code{prob(dose = numeric, model = LogisticKadaneBetaGamma, samples = Samples)}: diff --git a/tests/testthat/test-Model-methods.R b/tests/testthat/test-Model-methods.R index e02832338..f5b45dab1 100644 --- a/tests/testthat/test-Model-methods.R +++ b/tests/testthat/test-Model-methods.R @@ -1113,6 +1113,42 @@ test_that("prob-ProbitLogNormalRel throws the error when dose is not valid", { ) }) +## LogisticLogNormalGrouped ---- + +test_that("prob-LogisticLogNormalGrouped works as expected", { + model <- .DefaultLogisticLogNormalGrouped() + samples <- h_as_samples(list( + alpha0 = c(0, -1, 1, 2), + delta0 = c(0, 1, -1, 0), + alpha1 = c(0, 0.5, 1, -1), + delta1 = c(1, 0, -1, 2) + )) + + result <- prob(10, model, samples, group = "mono") + expect_equal(result, c(0.5, 0.5378, 0.9645, 0.4249), tolerance = 1e-4) +}) + +test_that("prob-LogisticLogNormalGrouped works as expected for scalar samples", { + model <- .DefaultLogisticLogNormalGrouped() + samples <- h_as_samples(list(alpha0 = 1, delta0 = -1, alpha1 = 1, delta1 = -0.5)) + + result <- prob(c(1, 30), model, samples, group = "combo") + expect_equal(result, c(0.5, 0.8456), tolerance = 1e-4) +}) + +test_that("prob-LogisticLogNormalGrouped works as expected for vectors", { + model <- .DefaultLogisticLogNormalGrouped() + samples <- h_as_samples(list( + alpha0 = c(1, 2), + delta0 = c(0.5, -0.5), + alpha1 = c(0, 1), + delta1 = c(1, 0.2) + )) + + result <- prob(c(1, 30), model, samples, group = c("mono", "combo")) + expect_equal(result, c(0.7311, 0.9962), tolerance = 1e-4) +}) + ## LogisticKadane ---- test_that("prob-LogisticKadane works as expected", { From c4c2342040c72390f695ae7ea7117a7159a57646 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Tue, 12 Sep 2023 20:31:32 +0200 Subject: [PATCH 04/31] progress --- R/Model-methods.R | 35 ++++++++++++++++++++++++++ tests/testthat/test-Model-methods.R | 39 +++++++++++++++++++++++++++++ 2 files changed, 74 insertions(+) diff --git a/R/Model-methods.R b/R/Model-methods.R index ae0259d2f..504a2b0e9 100644 --- a/R/Model-methods.R +++ b/R/Model-methods.R @@ -403,6 +403,41 @@ setMethod( } ) +## LogisticLogNormalGrouped ---- + +#' @describeIn dose method for [`LogisticLogNormalGrouped`] which needs `group` +#' argument in addition. +#' @param group (`character` or `factor`)\cr for [`LogisticLogNormalGrouped`], +#' indicating whether to calculate the dose for the `mono` or for +#' the `combo` arm. +#' @aliases dose-LogisticLogNormalGrouped +#' @export +#' +setMethod( + f = "dose", + signature = signature( + x = "numeric", + model = "LogisticLogNormalGrouped", + samples = "Samples" + ), + definition = function(x, model, samples, group) { + assert_probabilities(x) + assert_subset(c("alpha0", "delta0", "alpha1", "delta1"), names(samples)) + assert_length(x, len = size(samples)) + assert_multi_class(group, c("character", "factor")) + assert_subset(as.character(group), choices = c("mono", "combo")) + assert_length(group, len = size(samples)) + + alpha0 <- samples@data$alpha0 + delta0 <- samples@data$delta0 + alpha1 <- samples@data$alpha1 + delta1 <- samples@data$delta1 + ref_dose <- as.numeric(model@ref_dose) + is_combo <- as.integer(group == "combo") + exp((logit(x) - (alpha0 + is_combo * delta0)) / (alpha1 + is_combo * delta1)) * ref_dose + } +) + ## LogisticKadane ---- #' @describeIn dose compute the dose level reaching a specific target diff --git a/tests/testthat/test-Model-methods.R b/tests/testthat/test-Model-methods.R index f5b45dab1..d98fd2d6a 100644 --- a/tests/testthat/test-Model-methods.R +++ b/tests/testthat/test-Model-methods.R @@ -504,6 +504,45 @@ test_that("dose-ProbitLogNormalRel throws the error when x is not valid", { ) }) +## LogisticLogNormalGrouped ---- + +test_that("dose-LogisticLogNormalGrouped works as expected", { + model <- .DefaultLogisticLogNormalGrouped() + samples <- h_as_samples(list( + alpha0 = c(0.1, -1, 1, 2), + delta0 = c(0, 1, -1, 0), + alpha1 = c(0, 0.5, 1, -1), + delta1 = c(1, 0, -0.9, 2) + )) + + result_mono <- dose(0.5, model, samples, group = "mono") + result_combo <- dose(0.5, model, samples, group = "combo") + + expect_equal(result_mono, c(0, 7.3891, 0.3679, 7.3891), tolerance = 1e-4) + expect_equal(result_combo, c(0.9048, 1, 1, 0.1353), tolerance = 1e-4) +}) + +test_that("dose-LogisticLogNormalGrouped works as expected for scalar samples", { + model <- .DefaultLogisticLogNormalGrouped() + samples <- h_as_samples(list(alpha0 = 1, delta0 = -1, alpha1 = 1, delta1 = -0.5)) + + result <- dose(c(1, 30), model, samples, group = "combo") + expect_equal(result, c(0.5, 0.8456), tolerance = 1e-4) +}) + +test_that("dose-LogisticLogNormalGrouped works as expected for vectors", { + model <- .DefaultLogisticLogNormalGrouped() + samples <- h_as_samples(list( + alpha0 = c(1, 2), + delta0 = c(0.5, -0.5), + alpha1 = c(0, 1), + delta1 = c(1, 0.2) + )) + + result <- dose(c(1, 30), model, samples, group = c("mono", "combo")) + expect_equal(result, c(0.7311, 0.9962), tolerance = 1e-4) +}) + ## LogisticKadane ---- test_that("dose-LogisticKadane works as expected", { From 697bb7be65be4ae7f4a725508e8bfce7669563c6 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Wed, 13 Sep 2023 09:54:26 +0200 Subject: [PATCH 05/31] dose and prob function generics result now pass through add. args --- R/Model-methods.R | 8 +++--- man/dose.Rd | 11 ++++++++ tests/testthat/test-Model-methods.R | 40 ++++++++++++++++++++++++++--- 3 files changed, 51 insertions(+), 8 deletions(-) diff --git a/R/Model-methods.R b/R/Model-methods.R index 504a2b0e9..042575579 100644 --- a/R/Model-methods.R +++ b/R/Model-methods.R @@ -50,8 +50,8 @@ setMethod( data = model_params, options = McmcOptions(samples = NROW(model_params[[1]])) ) - function(x) { - dose(x = x, model = model, samples = samples) + function(x, ...) { + dose(x = x, model = model, samples = samples, ...) } } ) @@ -127,8 +127,8 @@ setMethod( data = model_params, options = McmcOptions(samples = NROW(model_params[[1]])) ) - function(dose) { - prob(dose = dose, model = model, samples = samples) + function(dose, ...) { + prob(dose = dose, model = model, samples = samples, ...) } } ) diff --git a/man/dose.Rd b/man/dose.Rd index 7c97c5d22..9ad8202fe 100644 --- a/man/dose.Rd +++ b/man/dose.Rd @@ -12,6 +12,8 @@ \alias{dose-ProbitLogNormal} \alias{dose,numeric,ProbitLogNormalRel,Samples-method} \alias{dose-ProbitLogNormalRel} +\alias{dose,numeric,LogisticLogNormalGrouped,Samples-method} +\alias{dose-LogisticLogNormalGrouped} \alias{dose,numeric,LogisticKadane,Samples-method} \alias{dose-LogisticKadane} \alias{dose,numeric,LogisticKadaneBetaGamma,Samples-method} @@ -50,6 +52,8 @@ dose(x, model, samples, ...) \S4method{dose}{numeric,ProbitLogNormalRel,Samples}(x, model, samples) +\S4method{dose}{numeric,LogisticLogNormalGrouped,Samples}(x, model, samples, group) + \S4method{dose}{numeric,LogisticKadane,Samples}(x, model, samples) \S4method{dose}{numeric,LogisticKadaneBetaGamma,Samples}(x, model, samples) @@ -88,6 +92,10 @@ as the sample.} used to compute the resulting doses. Can also be missing for some models.} \item{...}{model specific parameters when \code{samples} are not used.} + +\item{group}{(\code{character} or \code{factor})\cr for \code{\link{LogisticLogNormalGrouped}}, +indicating whether to calculate the dose for the \code{mono} or for +the \code{combo} arm.} } \value{ A \code{number} or \code{numeric} vector with the doses. @@ -133,6 +141,9 @@ probability of the occurrence of a DLE (\code{x}). \item \code{dose(x = numeric, model = ProbitLogNormalRel, samples = Samples)}: compute the dose level reaching a specific target probability of the occurrence of a DLE (\code{x}). +\item \code{dose(x = numeric, model = LogisticLogNormalGrouped, samples = Samples)}: method for \code{\link{LogisticLogNormalGrouped}} which needs \code{group} +argument in addition. + \item \code{dose(x = numeric, model = LogisticKadane, samples = Samples)}: compute the dose level reaching a specific target probability of the occurrence of a DLE (\code{x}). diff --git a/tests/testthat/test-Model-methods.R b/tests/testthat/test-Model-methods.R index d98fd2d6a..1930ff5d1 100644 --- a/tests/testthat/test-Model-methods.R +++ b/tests/testthat/test-Model-methods.R @@ -11,7 +11,7 @@ test_that("doseFunction-GeneralModel returns correct dose function", { dose_fun_env <- environment(dose_fun) expect_function(doseFunction, args = c("model", "..."), null.ok = FALSE) - expect_function(dose_fun, args = "x", nargs = 1, null.ok = FALSE) + expect_function(dose_fun, args = c("x", "..."), null.ok = FALSE) # Body of `dose_fun` must be a `dose` method with `x`, `model` and `samples` args. dose_fun_body <- as.list(body(dose_fun)[[2]]) @@ -52,7 +52,7 @@ test_that("doseFunction-GeneralModel returns correct dose function for matrix pa dose_fun_env <- environment(dose_fun) expect_function(doseFunction, args = c("model", "..."), null.ok = FALSE) - expect_function(dose_fun, args = "x", nargs = 1, null.ok = FALSE) + expect_function(dose_fun, args = c("x", "..."), null.ok = FALSE) # Body of `dose_fun` must be a `dose` method with `x`, `model` and `samples` args. dose_fun_body <- as.list(body(dose_fun))[[2]] @@ -119,6 +119,22 @@ test_that("doseFunction-ModelPseudo throws the error when no params are provided ) }) +## LogisticLogNormalGrouped ---- + +test_that("doseFunction-LogisticLogNormalGrouped works as expected", { + model <- .DefaultLogisticLogNormalGrouped() + + dose_fun <- expect_silent(doseFunction( + model, alpha0 = 1, delta0 = 0.5, alpha1 = 0.5, delta1 = -0.2 + )) + dose_fun <- h_covr_detrace(dose_fun) + + expect_function(dose_fun, args = c("x", "..."), null.ok = FALSE) + expect_error(dose_fun(1), "argument \"group\" is missing, with no default") + result <- expect_silent(dose_fun(0.5, group = "mono")) + expect_equal(result, 0.13534, tolerance = 1e-4) +}) + # probFunction ---- ## GeneralModel ---- @@ -132,7 +148,7 @@ test_that("probFunction-GeneralModel returns correct prob function", { prob_fun_env <- environment(prob_fun) expect_function(probFunction, args = c("model", "..."), null.ok = FALSE) - expect_function(prob_fun, args = "dose", nargs = 1, null.ok = FALSE) + expect_function(prob_fun, args = c("dose", "..."), null.ok = FALSE) # Body of `prob_fun` must be a `prob` method with `dose`, `model` and `samples` args. prob_fun_body <- as.list(body(prob_fun)[[2]]) @@ -173,7 +189,7 @@ test_that("probFunction-GeneralModel returns correct prob function for matrix pa prob_fun_env <- environment(prob_fun) expect_function(probFunction, args = c("model", "..."), null.ok = FALSE) - expect_function(prob_fun, args = "dose", nargs = 1, null.ok = FALSE) + expect_function(prob_fun, args = c("dose", "..."), null.ok = FALSE) # Body of `prob_fun` must be a `prob` method with `dose`, `model` and `samples` args. prob_fun_body <- as.list(body(prob_fun)[[2]]) @@ -239,6 +255,22 @@ test_that("probFunction-ModelTox throws the error when no params are provided", ) }) +## LogisticLogNormalGrouped ---- + +test_that("probFunction-LogisticLogNormalGrouped works as expected", { + model <- .DefaultLogisticLogNormalGrouped() + + prob_fun <- expect_silent(probFunction( + model, alpha0 = 1, delta0 = 0.5, alpha1 = 0.5, delta1 = -0.2 + )) + prob_fun <- h_covr_detrace(prob_fun) + + expect_function(prob_fun, args = c("dose", "..."), null.ok = FALSE) + expect_error(prob_fun(1), "argument \"group\" is missing, with no default") + result <- expect_silent(prob_fun(10, group = "mono")) + expect_equal(result, 0.8958, tolerance = 1e-4) +}) + # efficacyFunction ---- ## ModelEff ---- From b8717453ada5a683e713ad0bf3fee996d73bda55 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Wed, 13 Sep 2023 10:04:44 +0200 Subject: [PATCH 06/31] fit and plot method tests --- .../plot-samples-logisticlognormalgrouped.svg | 66 +++++++++++++++++++ tests/testthat/test-Samples-methods.R | 32 +++++++++ 2 files changed, 98 insertions(+) create mode 100644 tests/testthat/_snaps/Samples-methods/plot-samples-logisticlognormalgrouped.svg diff --git a/tests/testthat/_snaps/Samples-methods/plot-samples-logisticlognormalgrouped.svg b/tests/testthat/_snaps/Samples-methods/plot-samples-logisticlognormalgrouped.svg new file mode 100644 index 000000000..0ad3dfc68 --- /dev/null +++ b/tests/testthat/_snaps/Samples-methods/plot-samples-logisticlognormalgrouped.svg @@ -0,0 +1,66 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + +0 +25 +50 +75 +100 + + + + + + + + + +2.5 +5.0 +7.5 +10.0 +Dose level +Probability of DLT [%] + +Type + + + + +Estimate +95% Credible Interval +plot-Samples-LogisticLogNormalGrouped + + diff --git a/tests/testthat/test-Samples-methods.R b/tests/testthat/test-Samples-methods.R index e78424553..df22d5ec2 100644 --- a/tests/testthat/test-Samples-methods.R +++ b/tests/testthat/test-Samples-methods.R @@ -217,6 +217,23 @@ test_that("fit-Samples forwards additional arguments to prob inside", { expect_named(result, c("dose", "middle", "lower", "upper")) }) +## Samples-LogisticLogNormalGrouped ---- + +test_that("fit-Samples works specifically also for LogisticLogNormalGrouped", { + mcmcOptions <- McmcOptions(samples = 3) + samples <- Samples( + data = list(alpha0 = -1:1, delta0 = c(0, 1, -1), alpha1 = -1:1, delta1 = c(-1, 0, 2)), + options = mcmcOptions + ) + model <- .DefaultLogisticLogNormalGrouped() + emptyData <- Data(doseGrid = seq(10, 80, 10)) + + result <- expect_silent(fit(samples, model, emptyData, group = "combo")) + expect_data_frame(result) + expect_equal(nrow(result), length(emptyData@doseGrid)) + expect_named(result, c("dose", "middle", "lower", "upper")) +}) + ## Samples-DataModel ---- test_that("fit-Samples works correctly for dual models", { @@ -869,6 +886,21 @@ test_that("Check that plot-Samples-ModelTox works correctly", { vdiffr::expect_doppelganger("plot-Samples-ModelTox_showlegend-FALSE", actual1) }) +## Samples-LogisticLogNormalGrouped ---- + +test_that("plot-Samples works specifically also for LogisticLogNormalGrouped", { + mcmcOptions <- McmcOptions(samples = 3) + samples <- Samples( + data = list(alpha0 = -1:1, delta0 = c(0, 1, -1), alpha1 = -1:1, delta1 = c(-1, 0, 2)), + options = mcmcOptions + ) + model <- .DefaultLogisticLogNormalGrouped() + emptyData <- Data(doseGrid = seq(0.5, 10, by = 0.1)) + + result <- expect_silent(plot(samples, model, emptyData, group = "combo")) + vdiffr::expect_doppelganger("plot-Samples-LogisticLogNormalGrouped", result) +}) + ## Samples-DataDual ---- test_that("Check that plot-Samples-ModelEff fails gracefully with bad input", { From f2cb21fe05365df599408943434fa1b3a65adf2f Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Wed, 13 Sep 2023 10:07:17 +0200 Subject: [PATCH 07/31] remove dummy tests because now we have the actual ones --- tests/testthat/test-Samples-methods.R | 22 ---------------------- 1 file changed, 22 deletions(-) diff --git a/tests/testthat/test-Samples-methods.R b/tests/testthat/test-Samples-methods.R index df22d5ec2..427bb4faf 100644 --- a/tests/testthat/test-Samples-methods.R +++ b/tests/testthat/test-Samples-methods.R @@ -205,18 +205,6 @@ test_that("fit-Samples works correctly for tox-only models", { checkIt(seed = 789, lowerQuantile = 0.25, upperQuantile = 0.75) }) -test_that("fit-Samples forwards additional arguments to prob inside", { - mcmcOptions <- McmcOptions(samples = 3) - samples <- Samples(data = list(alpha0 = 1:3, alpha1 = 4:6), options = mcmcOptions) - model <- h_needs_extra_prob_model() - emptyData <- Data(doseGrid = seq(10, 80, 10)) - - result <- fit(samples, model, emptyData, extra_argument = "yes") - expect_data_frame(result) - expect_equal(nrow(result), length(emptyData@doseGrid)) - expect_named(result, c("dose", "middle", "lower", "upper")) -}) - ## Samples-LogisticLogNormalGrouped ---- test_that("fit-Samples works specifically also for LogisticLogNormalGrouped", { @@ -396,16 +384,6 @@ test_that("plot-Samples works correctly", { vdiffr::expect_doppelganger("plot-Samples_showLegend-FALSE", actual1) }) -test_that("plot-Samples forwards additional arguments to prob inside", { - mcmcOptions <- McmcOptions(samples = 3) - samples <- Samples(data = list(alpha0 = 1:3, alpha1 = 4:6), options = mcmcOptions) - model <- h_needs_extra_prob_model() - emptyData <- Data(doseGrid = seq(10, 80, 10)) - - result <- plot(samples, model, emptyData, extra_argument = "yes") - expect_list(result) -}) - test_that("plot-Samples-DualEndpoint fails gracefully with bad input", { data <- DataDual( x = c(0.1, 0.5, 1.5, 3, 6, 10, 10, 10, 20, 20, 20, 40, 40, 40, 50, 50, 50), From 75380124beefc64559a9079093590f061690f11a Mon Sep 17 00:00:00 2001 From: github-actions <41898282+github-actions[bot]@users.noreply.github.com> Date: Wed, 13 Sep 2023 08:11:57 +0000 Subject: [PATCH 08/31] [skip actions] Restyle files --- tests/testthat/test-Model-methods.R | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test-Model-methods.R b/tests/testthat/test-Model-methods.R index 1930ff5d1..0f9f27cf8 100644 --- a/tests/testthat/test-Model-methods.R +++ b/tests/testthat/test-Model-methods.R @@ -125,7 +125,8 @@ test_that("doseFunction-LogisticLogNormalGrouped works as expected", { model <- .DefaultLogisticLogNormalGrouped() dose_fun <- expect_silent(doseFunction( - model, alpha0 = 1, delta0 = 0.5, alpha1 = 0.5, delta1 = -0.2 + model, + alpha0 = 1, delta0 = 0.5, alpha1 = 0.5, delta1 = -0.2 )) dose_fun <- h_covr_detrace(dose_fun) @@ -261,7 +262,8 @@ test_that("probFunction-LogisticLogNormalGrouped works as expected", { model <- .DefaultLogisticLogNormalGrouped() prob_fun <- expect_silent(probFunction( - model, alpha0 = 1, delta0 = 0.5, alpha1 = 0.5, delta1 = -0.2 + model, + alpha0 = 1, delta0 = 0.5, alpha1 = 0.5, delta1 = -0.2 )) prob_fun <- h_covr_detrace(prob_fun) From d2eaa63a27d5c2b26471f276b4dda0d552c70d04 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Wed, 13 Sep 2023 10:14:56 +0200 Subject: [PATCH 09/31] remove indentation linter because it is not in defaults to avoid warn --- .lintr | 1 - 1 file changed, 1 deletion(-) diff --git a/.lintr b/.lintr index 0a19bc3fd..7049743a2 100644 --- a/.lintr +++ b/.lintr @@ -3,6 +3,5 @@ linters: linters_with_defaults( cyclocomp_linter = NULL, object_usage_linter = NULL, object_length_linter = NULL, - indentation_linter = NULL, object_name_linter = object_name_linter(c("CamelCase", "camelCase", "snake_case")) ) From 3a4fdf31d16efaaba6215c4d112492cb02ba03b0 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Wed, 13 Sep 2023 15:30:44 +0200 Subject: [PATCH 10/31] polishing --- R/crmPack-package.R | 5 ++++- _pkgdown.yaml | 1 + tests/testthat/test-Model-methods.R | 10 +++++----- 3 files changed, 10 insertions(+), 6 deletions(-) diff --git a/R/crmPack-package.R b/R/crmPack-package.R index 5a5385bf2..44fbf1ee0 100644 --- a/R/crmPack-package.R +++ b/R/crmPack-package.R @@ -81,7 +81,9 @@ globalVariables(c( "logit<-", "rho0", "alpha0", + "delta0", "alpha1", + "delta1", "inverse", "priorCov", "theta", @@ -135,7 +137,8 @@ globalVariables(c( "ref_dose", "comp", "X", - "skel_probs" + "skel_probs", + "is_combo" )) # nolint end diff --git a/_pkgdown.yaml b/_pkgdown.yaml index fdfdab0e1..91b3857a8 100644 --- a/_pkgdown.yaml +++ b/_pkgdown.yaml @@ -51,6 +51,7 @@ reference: - LogisticLogNormalSub - ProbitLogNormal - ProbitLogNormalRel + - LogisticLogNormalGrouped - LogisticKadane - LogisticKadaneBetaGamma - LogisticNormalMixture diff --git a/tests/testthat/test-Model-methods.R b/tests/testthat/test-Model-methods.R index 0f9f27cf8..7e4a5cc54 100644 --- a/tests/testthat/test-Model-methods.R +++ b/tests/testthat/test-Model-methods.R @@ -560,8 +560,8 @@ test_that("dose-LogisticLogNormalGrouped works as expected for scalar samples", model <- .DefaultLogisticLogNormalGrouped() samples <- h_as_samples(list(alpha0 = 1, delta0 = -1, alpha1 = 1, delta1 = -0.5)) - result <- dose(c(1, 30), model, samples, group = "combo") - expect_equal(result, c(0.5, 0.8456), tolerance = 1e-4) + result <- dose(c(0.2, 0.8), model, samples, group = "combo") + expect_equal(result, c(0.0625, 16), tolerance = 1e-4) }) test_that("dose-LogisticLogNormalGrouped works as expected for vectors", { @@ -569,12 +569,12 @@ test_that("dose-LogisticLogNormalGrouped works as expected for vectors", { samples <- h_as_samples(list( alpha0 = c(1, 2), delta0 = c(0.5, -0.5), - alpha1 = c(0, 1), + alpha1 = c(0.5, 1), delta1 = c(1, 0.2) )) - result <- dose(c(1, 30), model, samples, group = c("mono", "combo")) - expect_equal(result, c(0.7311, 0.9962), tolerance = 1e-4) + result <- dose(c(0.4, 0.8), model, samples, group = c("mono", "combo")) + expect_equal(result, c(0.0601, 0.9096), tolerance = 1e-4) }) ## LogisticKadane ---- From e43cbc05138b8c5739323796ba44d137cd2c3fff Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 14 Sep 2023 10:53:02 +0200 Subject: [PATCH 11/31] add DesignGrouped class --- NAMESPACE | 4 + R/Design-class.R | 82 ++++++++++++++++++ R/Design-validity.R | 16 ++++ _pkgdown.yaml | 1 + examples/Design-class-DesignGrouped.R | 55 ++++++++++++ man/DesignGrouped-class.Rd | 118 ++++++++++++++++++++++++++ man/v_design.Rd | 6 ++ tests/testthat/test-Design-class.R | 44 ++++++++++ 8 files changed, 326 insertions(+) create mode 100644 examples/Design-class-DesignGrouped.R create mode 100644 man/DesignGrouped-class.Rd diff --git a/NAMESPACE b/NAMESPACE index 5fc20e077..dc386d882 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -28,6 +28,7 @@ export(.DefaultCohortSizeParts) export(.DefaultCohortSizeRange) export(.DefaultDALogisticLogNormal) export(.DefaultDataGrouped) +export(.DefaultDesignGrouped) export(.DefaultDualEndpoint) export(.DefaultDualEndpointBeta) export(.DefaultDualEndpointEmax) @@ -86,6 +87,7 @@ export(.DefaultStoppingTargetBiomarker) export(.DefaultStoppingTargetProb) export(.DefaultTITELogisticLogNormal) export(.Design) +export(.DesignGrouped) export(.DualDesign) export(.DualEndpoint) export(.DualEndpointBeta) @@ -189,6 +191,7 @@ export(DataGrouped) export(DataMixture) export(DataParts) export(Design) +export(DesignGrouped) export(DualDesign) export(DualEndpoint) export(DualEndpointBeta) @@ -375,6 +378,7 @@ exportClasses(DataGrouped) exportClasses(DataMixture) exportClasses(DataParts) exportClasses(Design) +exportClasses(DesignGrouped) exportClasses(DualDesign) exportClasses(DualEndpoint) exportClasses(DualEndpointBeta) diff --git a/R/Design-class.R b/R/Design-class.R index a0687ccbc..3f6b0496d 100644 --- a/R/Design-class.R +++ b/R/Design-class.R @@ -529,3 +529,85 @@ 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 including the rules to be followed for the mono arm, the +#' `model` slot will be ignored since the joint model is used instead. +#' @slot combo (`Design`)\cr including the rules to be followed for the combo arm, the +#' `model` slot will be ignored since the joint model is used instead. +#' @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. +#' +#' @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-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/_pkgdown.yaml b/_pkgdown.yaml index 91b3857a8..2d24e74ef 100644 --- a/_pkgdown.yaml +++ b/_pkgdown.yaml @@ -128,6 +128,7 @@ reference: - DualDesign - TDsamplesDesign - TDDesign + - DesignGrouped - title: Internal Helper Functions contents: - h_doses_unique_per_cohort 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/man/DesignGrouped-class.Rd b/man/DesignGrouped-class.Rd new file mode 100644 index 000000000..40da3a1d1 --- /dev/null +++ b/man/DesignGrouped-class.Rd @@ -0,0 +1,118 @@ +% 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. +} +\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 including the rules to be followed for the mono arm, the +\code{model} slot will be ignored since the joint model is used instead.} + +\item{\code{combo}}{(\code{Design})\cr including the rules to be followed for the combo arm, the +\code{model} slot will be ignored since the joint model is used instead.} + +\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/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/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") +}) From 33071f7baff66fb96939c2baa89e35bb6d285a35 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 14 Sep 2023 11:02:48 +0200 Subject: [PATCH 12/31] also test validation fun --- tests/testthat/test-Design-validity.R | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) 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) +}) From 88366e2fd2ae440bfc3499f4343cd3bb2d82911e Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 14 Sep 2023 11:37:26 +0200 Subject: [PATCH 13/31] start simulate method --- R/Design-methods.R | 307 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 307 insertions(+) diff --git a/R/Design-methods.R b/R/Design-methods.R index 4828b4332..53d07c1be 100644 --- a/R/Design-methods.R +++ b/R/Design-methods.R @@ -4710,3 +4710,310 @@ 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 [setSeed()]. +#' @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` or `NULL`)\cr optional `data.frame` with arguments that work +#' for 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. In order to produce outcomes from the posterior predictive +#' distribution, e.g, pass an `object`that contains the data observed so +#' far, `truth` contains the `prob` function from the model in +#' `object`, and `args` contains posterior samples from the model. +#' @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 = NULL, + firstSeparate = FALSE, + mcmcOptions = McmcOptions(), + parallel = FALSE, + nCores = min(parallelly::availableCores(), 5), + ...) { + nsim <- safeInteger(nsim) + + stopifnot( + is.function(truth), + is.function(combo_truth), + is.scalar(nsim), + nsim > 0, + is.bool(parallel), + is.scalar(nCores), + nCores > 0 + ) + + args <- as.data.frame(args) + n_args <- max(nrow(args), 1L) + rng_state <- setSeed(seed) + sim_seeds <- sample.int(n = 2147483647, size = nsim) + + run_sim <- function(iter_sim) { + set.seed(sim_seeds[iter_sim]) + + this_args <- args[(iter_sim - 1) %% n_args + 1, , drop = FALSE] + this_mono_truth <- function(dose) { + do.call(truth, c(dose, this_args)) + } + this_combo_truth <- function(dose) { + do.call(combo_truth, c(dose, this_args)) + } + + # Start the simulated data with the provided one. + this_mono_data <- object@mono@data + this_combo_data <- object@combo@data + + # Shall we stop the trial? Separately for mono and combo. + stop_mono <- stop_combo <- FALSE + + # Are we in the first cohort? This is to support the staggering feature. + first_cohort <- TRUE + + # What are the next doses to be used? Initialize with starting doses. + if (object@mono@startingDose < object@combo@startingDose) { + warning("combo starting dose usually not higher than mono starting dose") + } + if (object@same_dose) { + this_mono_dose <- + this_combo_dose <- min(object@mono@startingDose, object@combo@startingDose) + } else { + this_mono_dose <- object@mono@startingDose + this_combo_dose <- object@combo@startingDose + } + + # Inside this loop we simulate the whole trial, until stopping. + while (!(stop_mono && stop_combo)) { + if (!stop_mono) { + this_mono_prob <- this_mono_truth(this_mono_dose) + this_mono_size <- size(object@mono@cohort_size, dose = this_mono_dose, data = this_mono_data) + this_mono_dlts <- if (firstSeparate && this_mono_size > 1) { + first_mono_dlts <- rbinom(n = 1, size = 1, prob = this_mono_prob) + if (first_mono_dlts == 0) { + c(first_mono_dlts, rbinom(n = this_mono_size - 1, size = 1, prob = this_mono_prob)) + } else { + first_mono_dlts + } + } else { + rbinom(n = this_mono_size, size = 1, prob = this_mono_prob) + } + this_mono_data <- update(object = this_mono_data, x = this_mono_dose, y = this_mono_dlts) + } + + if (!stop_combo && (!first_cohort || !object@first_cohort_mono_only)) { + this_combo_prob <- this_combo_truth(this_combo_dose) + this_combo_size <- size(object@combo@cohort_size, dose = this_combo_dose, data = this_combo_data) + this_combo_dlts <- if (firstSeparate && this_combo_size > 1) { + first_combo_dlts <- rbinom(n = 1, size = 1, prob = this_combo_prob) + if (first_combo_dlts == 0) { + c(first_combo_dlts, rbinom(n = this_combo_size - 1, size = 1, prob = this_combo_prob)) + } else { + first_combo_dlts + } + } else { + rbinom(n = this_combo_size, size = 1, prob = this_combo_prob) + } + this_combo_data <- update(object = this_combo_data, x = this_combo_dose, y = this_combo_dlts) + } + + if (first_cohort) { + first_cohort <- FALSE + } + + grouped_data <- groupData( + this_mono_data, + this_combo_data + ) + this_samples <- mcmc( + data = grouped_data, + model = object@model, + options = mcmcOptions + ) + + if (!stop_mono) { + mono_dose_limit <- maxDose( + object@mono@increments, + data = this_mono_data + ) + + this_mono_dose <- nextBest( + object@mono@nextBest, + doselimit = mono_dose_limit, + samples = this_samples, + model = object@model, + data = grouped_data, + group = factor("mono", levels = c("mono", "combo")) + )$value + + stop_mono <- stopTrial( + object@mono@stopping, + dose = this_mono_dose, + samples = this_samples, + model = object@model, + data = this_mono_data, + group = factor("mono", levels = c("mono", "combo")) + ) + + stop_mono_results <- h_unpack_stopit(stop_mono) + } + + if (!stop_combo) { + combo_dose_limit <- if (is.na(this_mono_dose)) { + 0 + } else { + combo_max_dose <- maxDose( + object@combo@increments, + data = this_combo_data + ) + min( + combo_max_dose, + this_mono_dose, + na.rm = TRUE + ) + } + + this_combo_dose <- nextBest( + object@combo@nextBest, + doselimit = combo_dose_limit, + samples = this_samples, + model = object@model, + data = grouped_data, + group = factor("combo", levels = c("mono", "combo")) + )$value + + stop_combo <- stopTrial( + object@combo@stopping, + dose = this_combo_dose, + samples = this_samples, + model = object@model, + data = this_combo_data, + group = factor("combo", levels = c("mono", "combo")) + ) + + stop_combo_results <- h_unpack_stopit(stop_combo) + + if (object@same_dose) { + this_mono_dose <- this_combo_dose <- min( + this_mono_dose, + this_combo_dose + ) + } + } + } + + fit_mono <- fit( + object = this_samples, + model = object@model, + data = grouped_data, + group = factor("mono", levels = c("mono", "combo")) + ) + fit_combo <- fit( + object = this_samples, + model = object@model, + data = grouped_data, + group = factor("combo", levels = c("mono", "combo")) + ) + + this_result <- list( + mono = list( + data = this_mono_data, + dose = this_mono_dose, + fit = subset(fit_mono, select = -dose), + stop = attr(stop_mono, "message"), + report_results = stop_mono_results + ), + combo = list( + data = this_combo_data, + dose = this_combo_dose, + fit = subset(fit_combo, select = -dose), + stop = attr(stop_combo, "message"), + report_results = stop_combo_results + ) + ) + + return(this_result) + } + + result_list <- getResultList( + fun = run_sim, + nsim = nsim, + vars = + c( + "simSeeds", + "args", + "nArgs", + "truth", + "combo_truth", + "firstSeparate", + "object", + "mcmcOptions" + ), + parallel = if (parallel) nCores else NULL + ) + + # Now we have a list with each element containing mono and combo, + # but we want it now the other way around, i.e. a list with 2 elements + # mono and combo and the iterations inside. + 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, "[[", "report_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 + ) + }) + } +) From 5931e405ebc9e84d6c6a643ebc06ce3f14837af4 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 14 Sep 2023 13:42:36 +0200 Subject: [PATCH 14/31] add example --- R/Design-methods.R | 2 +- .../Design-method-simulate-DesignGrouped.R | 74 ++++++++++ man/simulate-DesignGrouped-method.Rd | 139 ++++++++++++++++++ 3 files changed, 214 insertions(+), 1 deletion(-) create mode 100644 examples/Design-method-simulate-DesignGrouped.R create mode 100644 man/simulate-DesignGrouped-method.Rd diff --git a/R/Design-methods.R b/R/Design-methods.R index 53d07c1be..347232218 100644 --- a/R/Design-methods.R +++ b/R/Design-methods.R @@ -4858,7 +4858,7 @@ setMethod( first_cohort <- FALSE } - grouped_data <- groupData( + grouped_data <- h_group_data( this_mono_data, this_combo_data ) diff --git a/examples/Design-method-simulate-DesignGrouped.R b/examples/Design-method-simulate-DesignGrouped.R new file mode 100644 index 000000000..d278b1023 --- /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/simulate-DesignGrouped-method.Rd b/man/simulate-DesignGrouped-method.Rd new file mode 100644 index 000000000..0b7700767 --- /dev/null +++ b/man/simulate-DesignGrouped-method.Rd @@ -0,0 +1,139 @@ +% 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 = NULL, + 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[=setSeed]{setSeed()}}.} + +\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} or \code{NULL})\cr optional \code{data.frame} with arguments that work +for 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. In order to produce outcomes from the posterior predictive +distribution, e.g, pass an \code{object}that contains the data observed so +far, \code{truth} contains the \code{prob} function from the model in +\code{object}, and \code{args} contains posterior samples from the model.} + +\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]]) +} From 417c34f43c310ddbdb77eab71775df707f58ed10 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 14 Sep 2023 14:24:03 +0200 Subject: [PATCH 15/31] add tests --- tests/testthat/test-Design-methods.R | 145 +++++++++++++++++++++++++++ 1 file changed, 145 insertions(+) diff --git a/tests/testthat/test-Design-methods.R b/tests/testthat/test-Design-methods.R index 6478607ea..f5a7c01eb 100644 --- a/tests/testthat/test-Design-methods.R +++ b/tests/testthat/test-Design-methods.R @@ -134,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 ---- From ed686cf3e96cdd996b73d926bcbfeba4a28aefea Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 14 Sep 2023 17:43:56 +0200 Subject: [PATCH 16/31] forward ... from nextBest call to prob, dose, gain calls inside --- R/Rules-methods.R | 34 +++++++++++++++++----------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/R/Rules-methods.R b/R/Rules-methods.R index 8d5d32127..1133a6431 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 ) @@ -781,8 +781,8 @@ setMethod( # 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( From 4de798116aa895d26f3dfba880d713df0b6a5514 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 14 Sep 2023 17:51:56 +0200 Subject: [PATCH 17/31] same for stopTrial methods --- R/Rules-methods.R | 24 +++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/R/Rules-methods.R b/R/Rules-methods.R index 1133a6431..17fbaa007 100644 --- a/R/Rules-methods.R +++ b/R/Rules-methods.R @@ -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) { From 725edf9d21295616a590a190f455df1303b80bd9 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 14 Sep 2023 17:57:00 +0200 Subject: [PATCH 18/31] remove helper NeedsExtraProbModel as it is no longer needed --- tests/testthat/helper-model.R | 19 ------------------- tests/testthat/test-Rules-methods.R | 12 ++++++------ tests/testthat/test-helpers.R | 2 +- 3 files changed, 7 insertions(+), 26 deletions(-) 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-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, From 04e7f6a6263161a82bcef19f402793e8eee6d678 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 14 Sep 2023 18:04:13 +0200 Subject: [PATCH 19/31] add back indentation_linter because it is needed with 3.1 lintr pkg --- .lintr | 1 + 1 file changed, 1 insertion(+) 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")) From 76a05da38c0fbe084c3e827dfdedbd4b0fdb75aa Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 14 Sep 2023 18:07:04 +0200 Subject: [PATCH 20/31] add missing assertion --- R/Design-methods.R | 1 + 1 file changed, 1 insertion(+) diff --git a/R/Design-methods.R b/R/Design-methods.R index 347232218..9ec1c031c 100644 --- a/R/Design-methods.R +++ b/R/Design-methods.R @@ -4777,6 +4777,7 @@ setMethod( is.function(combo_truth), is.scalar(nsim), nsim > 0, + is.bool(firstSeparate), is.bool(parallel), is.scalar(nCores), nCores > 0 From ef01b23ae4d0b22f13ad7fca1aa4aa4d728af554 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 14 Sep 2023 18:08:15 +0200 Subject: [PATCH 21/31] add to pkgdown --- _pkgdown.yaml | 1 + 1 file changed, 1 insertion(+) diff --git a/_pkgdown.yaml b/_pkgdown.yaml index 59346e9d4..13dc8284e 100644 --- a/_pkgdown.yaml +++ b/_pkgdown.yaml @@ -454,6 +454,7 @@ reference: - simulate,RuleDesign-method - simulate,TDDesign-method - simulate,TDsamplesDesign-method + - simulate-DesignGrouped - summary,DualSimulations-method - summary,GeneralSimulations-method - summary,PseudoDualFlexiSimulations-method From 873bb0aed3c660a1602caabb44c0fbdc77a85891 Mon Sep 17 00:00:00 2001 From: github-actions <41898282+github-actions[bot]@users.noreply.github.com> Date: Thu, 14 Sep 2023 16:12:18 +0000 Subject: [PATCH 22/31] [skip actions] Restyle files --- examples/Design-method-simulate-DesignGrouped.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/Design-method-simulate-DesignGrouped.R b/examples/Design-method-simulate-DesignGrouped.R index d278b1023..f37b725f3 100644 --- a/examples/Design-method-simulate-DesignGrouped.R +++ b/examples/Design-method-simulate-DesignGrouped.R @@ -52,7 +52,7 @@ legend("topright", c("mono", "combo"), lty = c(1, 2), col = c(1, 2)) set.seed(123) my_sims <- simulate( my_design, - nsim = 4, # This should be at least 100 in actual applications. + nsim = 4, # This should be at least 100 in actual applications. seed = 123, truth = my_truth, combo_truth = my_combo_truth From b5f7e3452c2e7419e203247b659f0f5d8cf0cec2 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 14 Sep 2023 18:13:10 +0200 Subject: [PATCH 23/31] trigger checks From 3e958eb24a2462e642427dba9673c820335d16c7 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Thu, 14 Sep 2023 18:49:35 +0200 Subject: [PATCH 24/31] avoid in_sim being passed to dose --- R/Rules-methods.R | 2 +- man/nextBest.Rd | 2 +- man/simulate-DesignGrouped-method.Rd | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/R/Rules-methods.R b/R/Rules-methods.R index 17fbaa007..dd3af7367 100644 --- a/R/Rules-methods.R +++ b/R/Rules-methods.R @@ -777,7 +777,7 @@ 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). 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/simulate-DesignGrouped-method.Rd b/man/simulate-DesignGrouped-method.Rd index 0b7700767..dbe34c884 100644 --- a/man/simulate-DesignGrouped-method.Rd +++ b/man/simulate-DesignGrouped-method.Rd @@ -116,7 +116,7 @@ legend("topright", c("mono", "combo"), lty = c(1, 2), col = c(1, 2)) set.seed(123) my_sims <- simulate( my_design, - nsim = 4, # This should be at least 100 in actual applications. + nsim = 4, # This should be at least 100 in actual applications. seed = 123, truth = my_truth, combo_truth = my_combo_truth From 4a48ee308171912eca1c6b5537cc5a017d6ad77c Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Fri, 15 Sep 2023 15:07:23 +0200 Subject: [PATCH 25/31] address review comments --- DESCRIPTION | 1 + R/Design-class.R | 13 +- R/Design-methods.R | 269 +++++++-------------------- R/helpers_design.R | 25 +++ man/DesignGrouped-class.Rd | 14 +- man/simulate-DesignGrouped-method.Rd | 13 +- 6 files changed, 113 insertions(+), 222 deletions(-) create mode 100644 R/helpers_design.R 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/R/Design-class.R b/R/Design-class.R index 3f6b0496d..05db7632b 100644 --- a/R/Design-class.R +++ b/R/Design-class.R @@ -543,16 +543,21 @@ DADesign <- function(model, data, #' #' @slot model (`LogisticLogNormalGrouped`)\cr the model to be used, currently only one #' class is allowed. -#' @slot mono (`Design`)\cr including the rules to be followed for the mono arm, the -#' `model` slot will be ignored since the joint model is used instead. -#' @slot combo (`Design`)\cr including the rules to be followed for the combo arm, the -#' `model` slot will be ignored since the joint model is used instead. +#' @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 #' diff --git a/R/Design-methods.R b/R/Design-methods.R index 9ec1c031c..a9077907f 100644 --- a/R/Design-methods.R +++ b/R/Design-methods.R @@ -4728,13 +4728,10 @@ setMethod("simulate", #' 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` or `NULL`)\cr optional `data.frame` with arguments that work -#' for 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. In order to produce outcomes from the posterior predictive -#' distribution, e.g, pass an `object`that contains the data observed so -#' far, `truth` contains the `prob` function from the model in -#' `object`, and `args` contains posterior samples from the model. +#' @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. @@ -4764,240 +4761,100 @@ setMethod( seed = NULL, truth, combo_truth, - args = NULL, + 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) - stopifnot( - is.function(truth), - is.function(combo_truth), - is.scalar(nsim), - nsim > 0, - is.bool(firstSeparate), - is.bool(parallel), - is.scalar(nCores), - nCores > 0 - ) - - args <- as.data.frame(args) n_args <- max(nrow(args), 1L) rng_state <- setSeed(seed) sim_seeds <- sample.int(n = 2147483647, size = nsim) run_sim <- function(iter_sim) { set.seed(sim_seeds[iter_sim]) - - this_args <- args[(iter_sim - 1) %% n_args + 1, , drop = FALSE] - this_mono_truth <- function(dose) { - do.call(truth, c(dose, this_args)) - } - this_combo_truth <- function(dose) { - do.call(combo_truth, c(dose, this_args)) - } - + 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. - this_mono_data <- object@mono@data - this_combo_data <- object@combo@data - - # Shall we stop the trial? Separately for mono and combo. - stop_mono <- stop_combo <- FALSE - - # Are we in the first cohort? This is to support the staggering feature. - first_cohort <- TRUE - + 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@mono@startingDose < object@combo@startingDose) { - warning("combo starting dose usually not higher than mono starting dose") - } if (object@same_dose) { - this_mono_dose <- - this_combo_dose <- min(object@mono@startingDose, object@combo@startingDose) + current$mono$dose <- current$combo$dose <- min(object@mono@startingDose, object@combo@startingDose) } else { - this_mono_dose <- object@mono@startingDose - this_combo_dose <- object@combo@startingDose + current$mono$dose <- object@mono@startingDose + current$combo$dose <- object@combo@startingDose } - # Inside this loop we simulate the whole trial, until stopping. - while (!(stop_mono && stop_combo)) { - if (!stop_mono) { - this_mono_prob <- this_mono_truth(this_mono_dose) - this_mono_size <- size(object@mono@cohort_size, dose = this_mono_dose, data = this_mono_data) - this_mono_dlts <- if (firstSeparate && this_mono_size > 1) { - first_mono_dlts <- rbinom(n = 1, size = 1, prob = this_mono_prob) - if (first_mono_dlts == 0) { - c(first_mono_dlts, rbinom(n = this_mono_size - 1, size = 1, prob = this_mono_prob)) - } else { - first_mono_dlts - } - } else { - rbinom(n = this_mono_size, size = 1, prob = this_mono_prob) - } - this_mono_data <- update(object = this_mono_data, x = this_mono_dose, y = this_mono_dlts) + 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 (!stop_combo && (!first_cohort || !object@first_cohort_mono_only)) { - this_combo_prob <- this_combo_truth(this_combo_dose) - this_combo_size <- size(object@combo@cohort_size, dose = this_combo_dose, data = this_combo_data) - this_combo_dlts <- if (firstSeparate && this_combo_size > 1) { - first_combo_dlts <- rbinom(n = 1, size = 1, prob = this_combo_prob) - if (first_combo_dlts == 0) { - c(first_combo_dlts, rbinom(n = this_combo_size - 1, size = 1, prob = this_combo_prob)) - } else { - first_combo_dlts - } - } else { - rbinom(n = this_combo_size, size = 1, prob = this_combo_prob) - } - this_combo_data <- update(object = this_combo_data, x = this_combo_dose, y = this_combo_dlts) + 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 (first_cohort) { - first_cohort <- FALSE - } - - grouped_data <- h_group_data( - this_mono_data, - this_combo_data - ) - this_samples <- mcmc( - data = grouped_data, - model = object@model, - options = mcmcOptions - ) - + 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 (!stop_mono) { - mono_dose_limit <- maxDose( - object@mono@increments, - data = this_mono_data - ) - - this_mono_dose <- nextBest( - object@mono@nextBest, - doselimit = mono_dose_limit, - samples = this_samples, - model = object@model, - data = grouped_data, - group = factor("mono", levels = c("mono", "combo")) - )$value - - stop_mono <- stopTrial( - object@mono@stopping, - dose = this_mono_dose, - samples = this_samples, - model = object@model, - data = this_mono_data, - group = factor("mono", levels = c("mono", "combo")) - ) - - stop_mono_results <- h_unpack_stopit(stop_mono) + 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") |> + {\(x) x$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 (!stop_combo) { - combo_dose_limit <- if (is.na(this_mono_dose)) { + current$combo$limit <- if (is.na(current$mono$dose)) { 0 } else { - combo_max_dose <- maxDose( - object@combo@increments, - data = this_combo_data - ) - min( - combo_max_dose, - this_mono_dose, - na.rm = TRUE - ) - } - - this_combo_dose <- nextBest( - object@combo@nextBest, - doselimit = combo_dose_limit, - samples = this_samples, - model = object@model, - data = grouped_data, - group = factor("combo", levels = c("mono", "combo")) - )$value - - stop_combo <- stopTrial( - object@combo@stopping, - dose = this_combo_dose, - samples = this_samples, - model = object@model, - data = this_combo_data, - group = factor("combo", levels = c("mono", "combo")) - ) - - stop_combo_results <- h_unpack_stopit(stop_combo) - - if (object@same_dose) { - this_mono_dose <- this_combo_dose <- min( - this_mono_dose, - this_combo_dose - ) + 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") |> + {\(x) x$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 && !stop_mono && !stop_combo) { + current$mono$dose <- current$combo$dose <- min(current$mono$dose, current$combo$dose) } } - - fit_mono <- fit( - object = this_samples, - model = object@model, - data = grouped_data, - group = factor("mono", levels = c("mono", "combo")) - ) - fit_combo <- fit( - object = this_samples, - model = object@model, - data = grouped_data, - group = factor("combo", levels = c("mono", "combo")) + 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, + data = data, dose = dose, fit = subset(fit, select = -dose), stop = attr(stop, "message"), report_results = results ) - - this_result <- list( - mono = list( - data = this_mono_data, - dose = this_mono_dose, - fit = subset(fit_mono, select = -dose), - stop = attr(stop_mono, "message"), - report_results = stop_mono_results - ), - combo = list( - data = this_combo_data, - dose = this_combo_dose, - fit = subset(fit_combo, select = -dose), - stop = attr(stop_combo, "message"), - report_results = stop_combo_results - ) - ) - - return(this_result) } - - result_list <- getResultList( - fun = run_sim, - nsim = nsim, - vars = - c( - "simSeeds", - "args", - "nArgs", - "truth", - "combo_truth", - "firstSeparate", - "object", - "mcmcOptions" - ), - parallel = if (parallel) nCores else NULL - ) - - # Now we have a list with each element containing mono and combo, - # but we want it now the other way around, i.e. a list with 2 elements - # mono and combo and the iterations inside. + vars_needed <- c("simSeeds", "args", "nArgs", "truth", "combo_truth", "firstSeparate", "object", "mcmcOptions") + parallel_cores <- if (parallel) nCores else NULL + result_list <- getResultList(fun = run_sim, nsim = nsim, vars = vars_needed, parallel = parallel_cores) + # 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") diff --git a/R/helpers_design.R b/R/helpers_design.R new file mode 100644 index 000000000..67dc5f777 --- /dev/null +++ b/R/helpers_design.R @@ -0,0 +1,25 @@ +h_add_dlts <- function(data, + dose, + truth, + cohort_size, + first_separate) { + assert_class(data, "Data") + assert_number(dose) + assert_function(truth) + assert_count(cohort_size, positive = TRUE) + assert_flag(first_separate) + + prob <- truth(dose) + size <- size(object@mono@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(object = data, x = dose, y = dlts) +} diff --git a/man/DesignGrouped-class.Rd b/man/DesignGrouped-class.Rd index 40da3a1d1..d44ef56eb 100644 --- a/man/DesignGrouped-class.Rd +++ b/man/DesignGrouped-class.Rd @@ -36,17 +36,23 @@ DesignGrouped( \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 including the rules to be followed for the mono arm, the -\code{model} slot will be ignored since the joint model is used instead.} +\item{\code{mono}}{(\code{Design})\cr defines the dose escalation rules for the mono arm, see +details.} -\item{\code{combo}}{(\code{Design})\cr including the rules to be followed for the combo arm, the -\code{model} slot will be ignored since the joint model is used instead.} +\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 diff --git a/man/simulate-DesignGrouped-method.Rd b/man/simulate-DesignGrouped-method.Rd index dbe34c884..34fbe75d8 100644 --- a/man/simulate-DesignGrouped-method.Rd +++ b/man/simulate-DesignGrouped-method.Rd @@ -11,7 +11,7 @@ seed = NULL, truth, combo_truth, - args = NULL, + args = data.frame(), firstSeparate = FALSE, mcmcOptions = McmcOptions(), parallel = FALSE, @@ -32,13 +32,10 @@ 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} or \code{NULL})\cr optional \code{data.frame} with arguments that work -for 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. In order to produce outcomes from the posterior predictive -distribution, e.g, pass an \code{object}that contains the data observed so -far, \code{truth} contains the \code{prob} function from the model in -\code{object}, and \code{args} contains posterior samples from the model.} +\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 From a9b80bbca841aabae4a7dce64c579e7a9212a5d3 Mon Sep 17 00:00:00 2001 From: github-actions <41898282+github-actions[bot]@users.noreply.github.com> Date: Fri, 15 Sep 2023 13:11:56 +0000 Subject: [PATCH 26/31] [skip actions] Restyle files --- R/Design-methods.R | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/R/Design-methods.R b/R/Design-methods.R index a9077907f..fcedb3b58 100644 --- a/R/Design-methods.R +++ b/R/Design-methods.R @@ -4817,7 +4817,9 @@ setMethod( 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") |> - {\(x) x$value}() + { + \(x) x$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) @@ -4831,7 +4833,9 @@ setMethod( } current$combo$dose <- object@combo@nextBest |> nextBest(current$combo$limit, current$samples, object@model, current$grouped, group = "combo") |> - {\(x) x$value}() + { + \(x) x$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) From 4a4b897167b68d8237e70bd2329a2122eb9cc324 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Fri, 15 Sep 2023 15:13:13 +0200 Subject: [PATCH 27/31] few fixes --- R/Design-methods.R | 20 ++++++++------------ R/helpers_design.R | 4 ++-- 2 files changed, 10 insertions(+), 14 deletions(-) diff --git a/R/Design-methods.R b/R/Design-methods.R index fcedb3b58..c4deb8fc1 100644 --- a/R/Design-methods.R +++ b/R/Design-methods.R @@ -4813,18 +4813,16 @@ setMethod( 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 (!stop_mono) { + 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") |> - { - \(x) x$value - }() + 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 (!stop_combo) { + if (!current$combo$stop) { current$combo$limit <- if (is.na(current$mono$dose)) { 0 } else { @@ -4832,15 +4830,13 @@ setMethod( 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") |> - { - \(x) x$value - }() + 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 && !stop_mono && !stop_combo) { + if (object@same_dose && !current$mono$stop && !current$combo$stop) { current$mono$dose <- current$combo$dose <- min(current$mono$dose, current$combo$dose) } } @@ -4848,7 +4844,7 @@ setMethod( current$combo$fit <- fit(current$samples, object@model, current$grouped, group = "combo") lapply( X = current[c("mono", "combo")], FUN = with, - data = data, dose = dose, fit = subset(fit, select = -dose), stop = attr(stop, "message"), report_results = results + list(data = data, dose = dose, fit = subset(fit, select = -dose), stop = attr(stop, "message"), report_results = results) ) } vars_needed <- c("simSeeds", "args", "nArgs", "truth", "combo_truth", "firstSeparate", "object", "mcmcOptions") diff --git a/R/helpers_design.R b/R/helpers_design.R index 67dc5f777..6af449185 100644 --- a/R/helpers_design.R +++ b/R/helpers_design.R @@ -6,11 +6,11 @@ h_add_dlts <- function(data, assert_class(data, "Data") assert_number(dose) assert_function(truth) - assert_count(cohort_size, positive = TRUE) + assert_class(cohort_size, "CohortSize") assert_flag(first_separate) prob <- truth(dose) - size <- size(object@mono@cohort_size, dose = dose, data = data) + 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) { From 692761f50ba3b51a3f66ae2008f6d4e4e36a9aa3 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Fri, 15 Sep 2023 15:43:25 +0200 Subject: [PATCH 28/31] amend news, move helpers to separate file, adapt to main changes,tests --- NEWS.md | 2 +- R/Design-methods.R | 146 +-------------------------- R/helpers_design.R | 145 +++++++++++++++++++++++++- man/get_result_list.Rd | 6 +- man/h_add_dlts.Rd | 26 +++++ man/set_seed.Rd | 5 +- tests/testthat/test-Design-methods.R | 48 --------- tests/testthat/test-helpers_design.R | 100 ++++++++++++++++++ 8 files changed, 277 insertions(+), 201 deletions(-) create mode 100644 man/h_add_dlts.Rd create mode 100644 tests/testthat/test-helpers_design.R 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-methods.R b/R/Design-methods.R index 2ea7e588a..6ee5321ea 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 ## ============================================================ @@ -4775,7 +4634,7 @@ setMethod( assert_count(nCores, positive = TRUE) n_args <- max(nrow(args), 1L) - rng_state <- setSeed(seed) + rng_state <- set_seed(seed) sim_seeds <- sample.int(n = 2147483647, size = nsim) run_sim <- function(iter_sim) { @@ -4846,8 +4705,7 @@ setMethod( ) } vars_needed <- c("simSeeds", "args", "nArgs", "truth", "combo_truth", "firstSeparate", "object", "mcmcOptions") - parallel_cores <- if (parallel) nCores else NULL - result_list <- getResultList(fun = run_sim, nsim = nsim, vars = vars_needed, parallel = parallel_cores) + result_list <- get_result_list(fun = run_sim, nsim = nsim, vars = vars_needed, parallel = parallel, n_cores = 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"), diff --git a/R/helpers_design.R b/R/helpers_design.R index 6af449185..b30026b78 100644 --- a/R/helpers_design.R +++ b/R/helpers_design.R @@ -1,3 +1,146 @@ +#' 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, @@ -21,5 +164,5 @@ h_add_dlts <- function(data, } else { rbinom(n = size, size = 1, prob = prob) } - update(object = data, x = dose, y = dlts) + update(data, x = dose, y = dlts) } 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/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/tests/testthat/test-Design-methods.R b/tests/testthat/test-Design-methods.R index 1f2521829..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 ---- 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) +}) From daf83d2bdf7f99cf8e59752171fe69050534fd36 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Fri, 15 Sep 2023 15:46:47 +0200 Subject: [PATCH 29/31] fix line lengths --- R/Design-methods.R | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/R/Design-methods.R b/R/Design-methods.R index 6ee5321ea..0fc4bfd29 100644 --- a/R/Design-methods.R +++ b/R/Design-methods.R @@ -4701,11 +4701,14 @@ setMethod( 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"), report_results = results) + 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(fun = run_sim, nsim = nsim, vars = vars_needed, parallel = parallel, n_cores = nCores) + 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"), From 73b94f4f381d5a5c92cf95236e3540f8cafbf381 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Fri, 15 Sep 2023 20:37:05 +0200 Subject: [PATCH 30/31] fix typo --- R/Design-methods.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/Design-methods.R b/R/Design-methods.R index 0fc4bfd29..e5372c743 100644 --- a/R/Design-methods.R +++ b/R/Design-methods.R @@ -4720,7 +4720,7 @@ setMethod( 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, "[[", "report_results") + report_results <- lapply(this_list, "[[", "results") stop_report <- as.matrix(do.call(rbind, report_results)) Simulations( From ae051f44a457a42b0819af747ffb3dfe77925fa3 Mon Sep 17 00:00:00 2001 From: Daniel Sabanes Bove Date: Sat, 16 Sep 2023 21:34:57 +0200 Subject: [PATCH 31/31] fix warning and note --- R/Design-methods.R | 2 +- R/crmPack-package.R | 3 ++- man/simulate-DesignGrouped-method.Rd | 2 +- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/R/Design-methods.R b/R/Design-methods.R index e5372c743..7d07790b1 100644 --- a/R/Design-methods.R +++ b/R/Design-methods.R @@ -4580,7 +4580,7 @@ setMethod("simulate", #' #' @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 [setSeed()]. +#' @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`. 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/man/simulate-DesignGrouped-method.Rd b/man/simulate-DesignGrouped-method.Rd index 34fbe75d8..54c2c4880 100644 --- a/man/simulate-DesignGrouped-method.Rd +++ b/man/simulate-DesignGrouped-method.Rd @@ -24,7 +24,7 @@ \item{nsim}{(\code{number})\cr how many trials should be simulated.} -\item{seed}{(\code{RNGstate})\cr generated with \code{\link[=setSeed]{setSeed()}}.} +\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.