crm12Comb: R Package of Phase I/II Adaptive Design for Drug Combinations based on CRM Design
Install and import the crm12Comb package
Summary of crm12Comb
The package crm12Comb is used to simulate the combined phase I and phase II adaptive dose finding design based on continual reassessment method (CRM) for drug combinations, and output the operating characteristics for design evaluations using different input scenarios/skeletons/prior distributions/link functions/other settings. The crm12Comb package contains eight functions with one main function SIM_phase_I_II() to realize the whole process of the design and other helper functions including get_ordering() to obtain the complete orderings for both toxicity and efficacy under drug combinations, priorSkeketons() to generate toxicity/efficacy skeletons, rBin2Corr() to generate correlated binary data for (toxicity, efficacy) outcomes, toxicity_est() (or efficacy_est()) to estimate toxicity (or efficacy) for one enrolled patient or cohort of patients, and randomization_phase() (or maximization_phase()) to perform adaptive randomization (or maximization) phase for selecting dose combination for next patient or cohort of patients allocations. The corresponding documentations can be viewed using the following statements.
### main function, please use this function to run simulations
?SIM_phase_I_II
### helper functions need to run before simulations (pre-defined settings)
?priorSkeletons
### helper functions included in the main function SIM_phase_I_II()
?rBin2Corr
?ger_ordering
?toxicity_est
?efficacy_est
?randomization_phase
?maximization_phaseLink functions used in CRM design
The crm12Comb package is currently able to perform various link functions based on corresponding prior distributions for both toxicity and efficacy probabilities. Let πT(di) be the toxicity probability for the ith dose combination di under each toxicity ordering m and toxicity skeleton 0 < p1m < p2m < ⋯ < pIm < 1, πT(di) = P[Yj = 1|di] ≈ Fm(di, β), m = 1, …, M
Empiric model Fm(di, β) = diexp (β), with normal on β or Fm(di, β) = diβ, with gamma prior on β.
Hyperbolic tangent model $$F_m(d_i,\beta) = \left(\frac{\tanh(d_i)+1}{2} \right)^{\beta},$$ with exponential prior on β.
One-parameter logistic model $$F_m(d_i,\beta) = \frac{1}{1+\exp(-a_0-\exp(\beta)d_i)},$$ with fixed a0 and normal prior on β or $$F_m(d_i,\beta) = \frac{1}{1+\exp(-a_0-\beta d_i)},$$ with fixed a0 and gamma prior on β.
Two-parameter logistic model $$F_m(d_i,\beta) = \frac{1}{1+\exp(-\alpha-\exp(\beta)d_i)},$$ with normal prior on α and β or $$F_m(d_i,\beta) = \frac{1}{1+\exp(-\alpha-\beta d_i)},$$ with gamma prior on α and β.
Similar link function and prior distributions can be derived for efficacy probabilities.
Exponential form was implemented on parameter β under the normal prior distribution to reveal the monotonic nature of both dose-toxicity curve and dose-efficacy curve, while no exponential form was added for gamma prior distribution because all the values are positive for gamma distribution.
Using crm12Comb
We provided two parts of examples in this section, one is for how to implement the functions in crm12Comb package and the other one is how to obtain simulation results given a list of input values.
Functions in crm12Comb
Here is an example of how to use the crm12Comb to perform simulation studies given a pre-specified scenarios of true (toxicity, efficacy) probabilities among 100 iterations, where the priorSkeletons() function is used to generate the toxicity skeleton and efficacy skeleton and SIM_phase_I_II() function is used to run the simulation for 100 iterations for this example.
The pre-defined true (toxicity, efficacy) probabilities are defined by two doses A and B with 4 levels for each dose, denoting by 1 to 4 as toxicity and efficacy from low to high.
| Scenario | ||||
|---|---|---|---|---|
| Doses of A | ||||
| Doses of B | 1 | 2 | 3 | 4 |
| 4 | (0.20,0.24) | (0.24,0.35) | (0.35,0.45) | (0.40,0.50) |
| 3 | (0.12,0.18) | (0.16,0.22) | (0.22,0.35) | (0.25,0.40) |
| 2 | (0.06,0.10) | (0.10,0.15) | (0.14,0.25) | (0.20,0.35) |
| 1 | (0.02,0.05) | (0.04,0.10) | (0.08,0.15) | (0.12,0.32) |
Toxicity skeleton and efficacy skeleton are further needed to be generated before conducting the simulations. Here we assume the empiric link function with normal prior distribution.
library(crm12Comb)
# generate skeletons
DLT_skeleton <- priorSkeletons(updelta=0.025, target=0.3, npos=10, ndose=16,
model = "empiric", prior = "normal", beta_mean=0)
print(paste0("DLT skeleton is: ", paste(round(DLT_skeleton,3), collapse=", ")))## [1] "DLT skeleton is: 0.015, 0.026, 0.042, 0.063, 0.09, 0.123, 0.161, 0.204, 0.251, 0.3, 0.351, 0.402, 0.452, 0.501, 0.548, 0.592"Efficacy_skeleton <- priorSkeletons(updelta=0.025, target=0.5, npos=10, ndose=16,
model = "empiric", prior = "normal", beta_mean=0)
print(paste0("Efficacy skeleton is: ", paste(round(Efficacy_skeleton,3), collapse=", ")))## [1] "Efficacy skeleton is: 0.079, 0.111, 0.149, 0.192, 0.24, 0.291, 0.343, 0.396, 0.449, 0.5, 0.549, 0.595, 0.638, 0.678, 0.714, 0.747"According to the above pre-specified scenario and skeletons, we can run the simulation by 100 times through phase I/II adaptive design for drug combinations based on CRM and obtain the 9 operating characteristics shown below.
scenario <- matrix(c(0.02, 0.05,
0.04, 0.10,
0.08, 0.15,
0.12, 0.32,
0.06, 0.10,
0.10, 0.15,
0.14, 0.25,
0.20, 0.35,
0.12, 0.18,
0.16, 0.22,
0.22, 0.35,
0.25, 0.40,
0.20, 0.24,
0.24, 0.35,
0.35, 0.45,
0.40, 0.50), ncol=2, byrow = TRUE)
# simulate 100 trials under the same model and prior distribution
simRes <- SIM_phase_I_II(nsim=100, Nmax=40, DoseComb=scenario, input_doseComb_forMat=c(4,4),
input_type_forMat="matrix", input_Nphase=20,
input_DLT_skeleton=DLT_skeleton,
input_efficacy_skeleton=Efficacy_skeleton,
input_DLT_thresh=0.3, input_efficacy_thresh=0.3,
input_cohortsize=1, input_corr=0,
input_early_stopping_safety_thresh=0.33,
input_early_stopping_futility_thresh=0.2,
input_model="empiric", input_para_prior="normal",
input_beta_mean=0, input_beta_sd=sqrt(1.34),
input_theta_mean=0, input_theta_sd=sqrt(1.34),
random_seed=123)
print(paste0("Probability of recommending safe/ineffective combinations as ODC is ", simRes$prob_safe))## [1] "Probability of recommending safe/ineffective combinations as ODC is 0.07"## [1] "Probability of recommending target combinations as ODC is 0.86"print(paste0("Probability of recommending overly toxic combinations as ODC is ", simRes$prob_toxic))## [1] "Probability of recommending overly toxic combinations as ODC is 0.06"## [1] "Mean # of patients enrolled is 39.94"## [1] "Proportion of patients allocated to true ODC(s) is 0.652"## [1] "Proportion stopped for safety is 0"## [1] "Proportion stopped for futility is 0.01"## [1] "Observed DLT rate is 0.184"## [1] "Observed response rate is 0.305"After obtaining the simulation results, the data for each patient allocation with toxicity and efficacy outcomes under each trial is stored sequentially in a list named “datALL” in simRes. This allows a comprehensive visualization of sequential patient enrollment with toxicity and efficacy outcomes and a summary of patient allocation by dose combinations can be plotted.
The first plot is for the sequential patient enrollment of the third simulated trial, where the red (green) half dot denoted patient with toxicity (efficacy) outcome.
The second plot is for patient allocations by dose combination of the third simulated trial shown by a histogram.
# generate plots of patient allocations by dose levels of the first trial
patient_allocation_plot(simRes$datALL[[3]])
In addition, the overall optimal dose combination (ODC) selection can be plotted by a histogram among total of 100 simulations for this example.
Simulations of list of input values
If there are lists of input values for several parameters, we can use the for loop with the SIM_phase_I_II() function to obtain multiple sets of results.
We provided the examples of the six scenarios given by Wages and Conaway (2014). We included several different commonly considered values of the parameters including (1) maximum sample size N = c(40, 50, 60), (2) pre-specified toxicity and efficacy skeletons Skeleton = c(1, 2), (3) number of patients for determination of randomization phase nR = c(10, 20, 30), and (4) Association parameter for efficacy and toxicity corr = c(0, −2.049, 0.814) were based on simulation results given by Thall and Cook (2004). The results are saved as “examples_results.RData” in the crm12Comb package. The result data can directly be used to obtain plots comparing among four different link functions using sample_plot() function.
Here are the six scenarios.
| True (toxicity, efficacy) probabilities | ||||
|---|---|---|---|---|
| Doses of | Doses of B | |||
| Scenario | A | 1 | 2 | 3 |
| 1 | 3 | (0.08, 0.15) | (0.10, 0.20) | (0.18, 0.40) |
| 2 | (0.04, 0.10) | (0.06, 0.16) | (0.08, 0.20) | |
| 1 | (0.02, 0.05) | (0.04, 0.10) | (0.06, 0.15) | |
| 2 | 3 | (0.16, 0.20) | (0.25, 0.35) | (0.35, 0.50) |
| 2 | (0.10, 0.10) | (0.14, 0.25) | (0.20, 0.40) | |
| 1 | (0.06, 0.05) | (0.08, 0.10) | (0.12, 0.20) | |
| 3 | 3 | (0.24, 0.40) | (0.33, 0.50) | (0.40, 0.60) |
| 2 | (0.16, 0.20) | (0.22, 0.40) | (0.35, 0.50) | |
| 1 | (0.08, 0.10) | (0.14, 0.25) | (0.20, 0.35) | |
| 4 | 3 | (0.33, 0.50) | (0.40, 0.60) | (0.55, 0.70) |
| 2 | (0.18, 0.35) | (0.25, 0.45) | (0.42, 0.55) | |
| 1 | (0.12, 0.20) | (0.20, 0.40) | (0.35, 0.50) | |
| 5 | 3 | (0.45, 0.55) | (0.55, 0.65) | (0.75, 0.75) |
| 2 | (0.20, 0.36) | (0.35, 0.49) | (0.40, 0.62) | |
| 1 | (0.15, 0.20) | (0.20, 0.35) | (0.25, 0.50) | |
| 6 | 3 | (0.65, 0.60) | (0.80, 0.65) | (0.85, 0.70) |
| 2 | (0.55, 0.55) | (0.70, 0.60) | (0.75, 0.65) | |
| 1 | (0.50, 0.50) | (0.55, 0.55) | (0.65, 0.60) |
The code below shows how to input the true probabilities in the form of matrix with scenario 1 as an example.
scenario1 <- matrix(c(0.02, 0.05,
0.04, 0.10,
0.06, 0.15,
0.04, 0.10,
0.06, 0.16,
0.08, 0.20,
0.08, 0.15,
0.10, 0.20,
0.18, 0.40), ncol=2, byrow = TRUE)The example code shows an approach to obtain the results with list of input values. Based on the scenarios, we first construct the total of 54 in consist of conditions 2 pre-specified toxicity and efficacy skeletons, 3 maximum sample sizes, 3 numbers of patients for determination of randomization phase and 3 association parameters for efficacy and toxicity.
orderings <- function(DLT1, DLT2, ORR1, ORR2){
input_Nphase <- c(10, 20, 30)
input_corr <- c(0, -2.049, 0.814)
input_N <- c(40, 50, 60)
DLTs <- list(DLT1, DLT2)
ORRs <- list(ORR1, ORR2)
conds <- list()
i <- 1
conds <- list()
i <- 1
for (s in 1:2){
for (n in 1:3){
for (np in 1:3){
for (c in 1:3){
conds[[i]] <- list(DLT=DLTs[[s]], ORR=ORRs[[s]], sklnum = s,
N=input_N[n], Nphase=input_Nphase[np],
corr=input_corr[c])
i <- i+1
}
}
}
}
return(conds)
}
SC <- list(scenario1, scenario2, scenario3, scenario4, scenario5, scenario6)The second step is to build the output data set with detailed names of 9 operating characteristics.
output <- data.frame(Scenario = double(), Skeleton = double(),
N = double(), nR = double(), corr = double(), safe = double(),
target = double(), toxic = double(), avgSS = double(),
prop_ODC = double(), stop_safety = double(), stop_futility = double(),
o_DLT = double(), o_ORR = double())
colnames(output) <- c("Scenario", "Skeleton", "N", "nR", "corr",
"Probability of recommending safe/ineffective combinations as ODC",
"Probability of recommending target combinations as ODC",
"Probability of recommending overly toxic combinations as ODC",
"Mean # of patients enrolled", "Proportion of patients allocated to true ODC(s)",
"Proportion stopped for safety", "Proportion stopped for futility",
"Observed DLT rate", "Observed response rate")The third step is to specify two sets of toxicity and efficacy skeletons so that we can plug the skeletons into the orderings() function to stored 54 conditions into a list named conds.
# empiric, normal prior
DLT_skeleton1 <- priorSkeletons(updelta=0.045, target=0.3, npos=5, ndose=9,
model = "empiric", prior = "normal", beta_mean=0)
DLT_skeleton2 <- priorSkeletons(updelta=0.06, target=0.3, npos=4, ndose=9,
model = "empiric", prior = "normal", beta_mean=0)
Efficacy_skeleton1 <- priorSkeletons(updelta=0.045, target=0.5, npos=5, ndose=9,
model = "empiric", prior = "normal", beta_mean=0)
Efficacy_skeleton2 <- priorSkeletons(updelta=0.06, target=0.5, npos=4, ndose=9,
model = "empiric", prior = "normal", beta_mean=0)
conds <- orderings(DLT1=DLT_skeleton1, DLT2=DLT_skeleton2,
ORR1=Efficacy_skeleton1, ORR2=Efficacy_skeleton2)The last step is to iterate over the 6 scenarios and 54 conditions each with 1000 simulations, storing the results into the output data set after each iteration.
for (s in 1:length(SC)){
for (c in 1:length(conds)){
print(paste0("Scenario=", s, ", skeleton=", conds[[c]]$sklnum,
", N=", conds[[c]]$N, ", nR=", conds[[c]]$Nphase,
", corr=", conds[[c]]$corr))
curr = SIM_phase_I_II(nsim=1000, Nmax=conds[[c]]$N, DoseComb=SC[[s]],
input_doseComb_forMat=c(3,3),
input_type_forMat="matrix",
input_Nphase=conds[[c]]$Nphase,
input_DLT_skeleton=conds[[c]]$DLT,
input_efficacy_skeleton=conds[[c]]$ORR,
input_DLT_thresh=0.3, input_efficacy_thresh=0.3,
input_cohortsize=1, input_corr=conds[[c]]$corr,
input_early_stopping_safety_thresh=0.33,
input_early_stopping_futility_thresh=0.2,
input_model="empiric", input_para_prior="normal",
input_beta_mean=0, input_beta_sd=sqrt(1.34),
input_theta_mean=0, input_theta_sd=sqrt(1.34),
random_seed=42)
currTmp = data.frame(s, conds[[c]]$sklnum, conds[[c]]$N, conds[[c]]$Nphase, conds[[c]]$corr,
curr$prob_safe, curr$prob_target, curr$prob_toxic, curr$mean_SS, curr$mean_ODC,
curr$prob_stop_safety, curr$prob_stop_futility, curr$mean_DLT, curr$mean_ORR)
output = rbind(output, currTmp)
}
}The process above is an example of using one link function “empiric” with normal prior. We also ran other link functions including “tanh” with exponential prior, “one-parameter logistic” with normal and gamma priors. We stored the corresponding results into a data set located at inst/extdata/example_results.RData in the crm12Comb package.
Here are two examples of plotting the result data by fixing three parameters and drawing the relationship between the specific operating characteristic and the non-fixed parameter. The first example showing the relationship between “Probability of ODC as target combinations” versus (3) “number of patients for determination of randomization phase” after fixing (1), (2) and (4).
sample_plot(examples_results, outcome = "prob_target",
outname = "Probability of ODC as target combinations",
N = 40, nR = NULL, Skeleton = 1, corr = 0)
The second example showing the relationship between “Proportion of patients allocated to true ODC(s)” versus (4) “Association parameter for efficacy and toxicity” after fixing (1), (2) and (3).
sample_plot(examples_results, outcome = "mean_ODC",
outname = "Proportion of patients allocated to true ODC(s)",
N = 60, nR = 30, Skeleton = 2, corr = NULL)
References
Murtaugh, P. A., & Fisher, L. D. (1990). Bivariate binary models of efficacy and toxicity in dose-ranging trials. Communications in Statistics-Theory and Methods, 19(6), 2003-2020. https://doi.org/10.1080/03610929008830305
O’Quigley, J., Pepe, M., & Fisher, L. (1990). Continual reassessment method: a practical design for phase 1 clinical trials in cancer. Biometrics, 33-48. https://doi.org/10.2307/2531628
Thall, P. F., & Cook, J. D. (2004). Dose‐finding based on efficacy–toxicity trade‐offs. Biometrics, 60(3), 684-693. https://doi.org/10.1111/j.0006-341X.2004.00218.x
Lee, S. M., & Cheung, Y. K. (2009). Model calibration in the continual reassessment method. Clinical Trials, 6(3), 227-238. https://doi.org/10.1177/1740774509105076
Wages, N. A., & Conaway, M. R. (2013). Specifications of a continual reassessment method design for phase I trials of combined drugs. Pharmaceutical statistics, 12(4), 217-224. https://doi.org/10.1002/pst.1575
Wages, N. A., & Conaway, M. R. (2014). Phase I/II adaptive design for drug combination oncology trials. Statistics in medicine, 33(12), 1990-2003. https://doi.org/10.1002/sim.6097