Comparing dose-escalation designs by simulation
Source:vignettes/A710-SimulationComparison.Rmd
A710-SimulationComparison.RmdThe escalation package by Kristian Brock. Documentation
is hosted at https://brockk.github.io/escalation/
Introduction
This vignette focuses on the common task of comparing competing
dose-escalation designs by simulation. Before reading on, be sure you
have read the README
file for a general explanation of how to compose dose-finding
designs in escalation; and the Simulation vignette for a general
introduction to using simulation.
Comparing competing designs
In simulate_compare, we implement the method of Sweeting et al. (2024) to efficiently compare
dose-finding designs. The crux of this method is to ensure that the same
simulated patients are used for each competing design so that, for
instance, if a patient is given dose 2 by one design and experiences
toxicity, then that patient will also experience toxicity if given dose
2 or above by any other design. Ensuring consistency across simulated
iterates in this way reduces Monte Carlo error and allows much faster
identification of differences between designs.
For example, let us compare the behaviour of the 3+3 and CRM designs investigated above. We start by defining the competing designs in a list with convenient names:
library(escalation)
target <- 0.25
skeleton <- c(0.05, 0.1, 0.25, 0.4, 0.6)
designs <-
list(
"3+3" = get_three_plus_three(num_doses = 5),
"CRM" = get_dfcrm(skeleton = skeleton, target = target) %>%
stop_at_n(n = 12),
"Stopping CRM" = get_dfcrm(skeleton = skeleton, target = target) %>%
stop_at_n(n = 12) %>%
stop_when_too_toxic(dose = 1, tox_threshold = 0.35, confidence = 0.8)
)Here we will compare three designs: 3+3; CRM without a toxicity stopping rule; and an otherwise identical CRM design with a toxicity stopping rule. The names we provide will be reused.
For illustration we use only a small number of replicates:
num_sims <- 20
true_prob_tox <- c(0.12, 0.27, 0.44, 0.53, 0.57)
sims <- simulate_compare(
designs,
num_sims = num_sims,
true_prob_tox = true_prob_tox
)
#> Running 3+3
#> Running CRM
#> Running Stopping CRMWe provide a convenient function to quickly visualise how the probability of selecting each dose in each design evolved as the simulations progressed:
convergence_plot(sims)
We can see immediately, for instance, that the designs generally agree that dose 2 is the best MTD candidate, and that the CRM designs are much more likely to recommend dose 3. We can be more precise by formally contrasting the probability of selecting each dose for each pair of designs:
library(dplyr)
library(ggplot2)
as_tibble(sims) %>%
filter(n %% 5 == 0) %>%
ggplot(aes(x = n, y = delta)) +
geom_point(size = 0.4) +
geom_linerange(aes(ymin = delta_l, ymax = delta_u)) +
geom_hline(yintercept = 0, linetype = "dashed", col = "red") +
facet_grid(comparison ~ dose,
labeller = labeller(
.rows = label_both,
.cols = label_both)
)
The error bars here reflect 95% symmetric asymptotic normal
confidence intervals. Change the alpha = 0.05 parameter
when calling as_tibble(sims) to get confident intervals for
a different significance level. We see that, even with the very small
sample size of 50 simulated trials, the CRM designs are significantly
more likely to recommend dose 3 than the 3+3 design. In contrast, in
this scenario there is very little difference at all between the two CRM
variants.
Working with PatientSamples
The idea at the core of the method by Sweeting et al. (2024) for efficient comparison of competing designs is to use the same patients across different designs. This reduces Monte Carlo error by examining where the designs differ in their recommendations given identical inputs.
This is achieved in escalation using classes inheriting
from PatientSample. A single PatientSample
reflects one particular state of the world, where patient \(i\) would reliably experience a toxicity or
efficacy event if treated at a particular dose. The consistent
occurrence of events is managed using latent uniform random variables.
We show how to import and export these variables, and use them to infer
potential outcomes at a range of doses.
To illustrate some of the fine-grained features of working with patient samples, we reproduce the BOIN12 example in Sweeting et al. (2024):
num_doses <- 5
designs <- list(
"BOIN12 v1" = get_boin12(num_doses = num_doses,
phi_t = 0.35, phi_e = 0.25,
u2 = 40, u3 = 60,
c_t = 0.95, c_e = 0.9) %>%
stop_when_n_at_dose(n = 12, dose = 'any') %>%
stop_at_n(n = 36),
"BOIN12 v2" = get_boin12(num_doses = num_doses,
phi_t = 0.35, phi_e = 0.25,
u2 = 40, u3 = 60,
c_t = 0.85, c_e = 0.8
) %>%
stop_when_n_at_dose(n = 12, dose = 'any') %>%
stop_at_n(n = 36)
)Exporting the underlying latent variables
By setting return_patient_samples = TRUE in the call to
simulate_compare:
true_prob_tox <- c(0.05, 0.10, 0.15, 0.18, 0.45)
true_prob_eff <- c(0.40, 0.50, 0.52, 0.53, 0.53)
set.seed(2024)
x <- simulate_compare(
designs = designs,
num_sims = 50,
true_prob_tox = true_prob_tox,
true_prob_eff = true_prob_eff,
return_patient_samples = TRUE
)
#> Running BOIN12 v1
#> Running BOIN12 v2we ensure the generated patient samples are exported in the returned object:
ps <- attr(x, "patient_samples")There are 50 patient-samples because we performed 50 simulated iterates:
length(ps)
#> [1] 50Each iterate is associated with one patient-sample.
For instance, the uniformly-distributed random variables that determine the occurrence of toxicity events in the first simulated iterate are:
ps[[1]]$tox_u
#> [1] 0.69817312 0.70142031 0.30320913 0.44945646 0.11031108 0.13106918
#> [7] 0.29230717 0.81558326 0.02112692 0.96611989 0.64736749 0.91045726
#> [13] 0.65435328 0.06133627 0.45295471 0.19192290 0.27820403 0.42659125
#> [19] 0.88575333 0.30889712 0.65158215 0.03589462 0.22038706 0.24535641
#> [25] 0.11568014 0.55189023 0.62894470and the equivalent for efficacy events are:
ps[[1]]$eff_u
#> [1] 0.4570092149 0.4157109936 0.8765880703 0.8422317195 0.8632152062
#> [6] 0.2954267922 0.6152471751 0.0005296487 0.4411276847 0.0221378489
#> [11] 0.8728519802 0.6238860050 0.6504830588 0.7410900716 0.4894901910
#> [16] 0.4564330971 0.0749628954 0.3724637846 0.7188321317 0.2872392512
#> [21] 0.1984391382 0.6658069612 0.1914782405 0.4981285664 0.8587655805
#> [26] 0.9091297414 0.0394975643These values reflect the patient-specific propensities to toxicity and efficacy. For instance, we see that the first patient would always experience a toxicity if treated at a dose with associated true toxicity probability less than 0.6981731. Likewise, the same patient would experience efficacy if treated at a dose with associated true efficacy probability less than 0.4570092.
These latent variables can be exported using R’s many I/O functions and formats.
Calculating potential outcomes
We can calculate all potential outcomes for a list of patient-samples
by calling get_potential_outcomes and providing the true
event probabilities:
z <- get_potential_outcomes(
patient_samples = ps,
true_prob_tox = true_prob_tox,
true_prob_eff = true_prob_eff
)Note that were we working with tox-only dose-escalation designs like
CRM or mTPI (for example), we would omit the true_prob_eff
parameter.
For instance, the potential outcomes in the first simulated iterate are:
z[[1]]
#> $tox
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 0 0 0 0 0
#> [2,] 0 0 0 0 0
#> [3,] 0 0 0 0 1
#> [4,] 0 0 0 0 1
#> [5,] 0 0 1 1 1
#> [6,] 0 0 1 1 1
#> [7,] 0 0 0 0 1
#> [8,] 0 0 0 0 0
#> [9,] 1 1 1 1 1
#> [10,] 0 0 0 0 0
#> [11,] 0 0 0 0 0
#> [12,] 0 0 0 0 0
#> [13,] 0 0 0 0 0
#> [14,] 0 1 1 1 1
#> [15,] 0 0 0 0 0
#> [16,] 0 0 0 0 1
#> [17,] 0 0 0 0 1
#> [18,] 0 0 0 0 1
#> [19,] 0 0 0 0 0
#> [20,] 0 0 0 0 1
#> [21,] 0 0 0 0 0
#> [22,] 1 1 1 1 1
#> [23,] 0 0 0 0 1
#> [24,] 0 0 0 0 1
#> [25,] 0 0 1 1 1
#> [26,] 0 0 0 0 0
#> [27,] 0 0 0 0 0
#>
#> $eff
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 0 1 1 1 1
#> [2,] 0 1 1 1 1
#> [3,] 0 0 0 0 0
#> [4,] 0 0 0 0 0
#> [5,] 0 0 0 0 0
#> [6,] 1 1 1 1 1
#> [7,] 0 0 0 0 0
#> [8,] 1 1 1 1 1
#> [9,] 0 1 1 1 1
#> [10,] 1 1 1 1 1
#> [11,] 0 0 0 0 0
#> [12,] 0 0 0 0 0
#> [13,] 0 0 0 0 0
#> [14,] 0 0 0 0 0
#> [15,] 0 1 1 1 1
#> [16,] 0 1 1 1 1
#> [17,] 1 1 1 1 1
#> [18,] 1 1 1 1 1
#> [19,] 0 0 0 0 0
#> [20,] 1 1 1 1 1
#> [21,] 1 1 1 1 1
#> [22,] 0 0 0 0 0
#> [23,] 1 1 1 1 1
#> [24,] 0 1 1 1 1
#> [25,] 0 0 0 0 0
#> [26,] 0 0 0 0 0
#> [27,] 1 1 1 1 1We see that patient 1 would not have experienced a toxicity at any dose because their underlying toxicity propensity variable was quite high at 0.6981731. Gladly, however, the same patient would have experienced efficacy at any dose at or above dose-level 2.
Importing the underlying latent variables
We show above how to export the latent event variables to store of use them elsewhere. We can also import the latent variables, to mimic the patient population used in some external simulation, for instance.
To do so, we instantiate a list of PatientSample objects, with one for each simulated iterate. We will work with ten iterates in the interests of speed:
num_sims <- 10
ps <- lapply(1:num_sims, function(x) PatientSample$new())We then call the set_eff_and_tox function on each to set
the latent toxicity and efficacy propensities to the desired values. For
illustration, let us simply sample uniform random variables and use
those. In a more realistic example, we would set the latent variables to
the values we wished to import, i.e. those values that were used in some
external study or idealised values that we wish to use.
set.seed(2024)
for(i in seq_len(num_sims)) {
tox_u_new <- runif(n = 50)
eff_u_new <- runif(n = 50)
ps[[i]]$set_eff_and_tox(tox_u = tox_u_new, eff_u = eff_u_new)
}We have set the latent variable vectors to have length 50. This means we can work with up to 50 patients in each iterate. If the simulation tries to use a 51st patient, we will receive an error. When setting the latent variables in this way, ensure you use enough values to cover your maximum sample size.
To use the patient-samples we have created in simulated trials, we
specify the patient_samples parameter:
x1 <- simulate_compare(
designs = designs,
num_sims = length(ps),
true_prob_tox = true_prob_tox,
true_prob_eff = true_prob_eff,
patient_samples = ps
)
#> Running BOIN12 v1
#> Running BOIN12 v2Having specified the exact patient population in this way, if we run the simulation study a second time with the same patients:
x2 <- simulate_compare(
designs = designs,
num_sims = length(ps),
true_prob_tox = true_prob_tox,
true_prob_eff = true_prob_eff,
patient_samples = ps
)
#> Running BOIN12 v1
#> Running BOIN12 v2we will see that the decisions within design are (usually) identical within a design. E.g.
all(
recommended_dose(x1[["BOIN12 v1"]]) ==
recommended_dose(x2[["BOIN12 v1"]])
)
#> [1] TRUEThe first BOIN12 variant recommends the exact same dose in the two batches because the patients are the same.
For completeness, when might the decisions vary despite the same patients being used? In a design that has inherent randomness. E.g. Wages & Tait’s design features adaptive randomisation, meaning that, even with identical patients, different behaviours will be seen.
Correlated toxicity and efficacy outcomes
Finally, let us introduce the CorrelatedPatientSample to
sample correlated toxicity and efficacy events. To coerce a correlation
of 0.5 between the latent uniform tox and eff event
propensities, we run:
ps <- CorrelatedPatientSample$new(num_patients = 100, rho = 0.5)We can observe
cor(ps$tox_u, ps$eff_u)
#> [1] 0.4827035that the correlation is close to desired. Note however, that the correlation of the binary level events is not, in general, equal to the target value of rho, because they depend on the event probabilities:
Further refinements
Check out the Further refinements section of the Simulation vignette for further ways of
refining the behaviour of simulations, including changing the cohort
size and timing of patient arrivals, simulating the conclusion of
partly-observed trials, tweaking the immediate next dose, returning all
interim model fits, and working with large trials. The methods described
there apply to both simulate_trials and
simulate_compare.