The OncoBayes2
package provides flexible functions for Bayesian meta-analytic modeling of the incidence of Dose Limiting Toxicities (DLTs) by dose level, under treatment regimes involving any number of combination partners. Such models may be used to support dose-escalation decisions and estimation of the Maximum Tolerated Dose (MTD) in adaptive Bayesian dose-escalation designs.
The package supports incorporation of historical data through a Meta-Analytic-Combined (MAC) framework [1], which stratifies these heterogeneous sources of information through a hierarchical model. Additionally, it allows the use of EXchangeable/Non-EXchangeable (EX/NEX) priors to manage the amount of information-sharing across subgroups of historical and/or concurrent sources of data.
Consider the application described in Section 3.2 of [1], in which the risk of DLT is to be studied as a function of dose for two drugs, drug A and drug B. Historical information on the toxicity profiles of these two drugs is available from single agent trials trial_A
and trial_B
. The historical data for this example is available in an internal data set.
group_id | drug_A | drug_B | num_patients | num_toxicities | cohort_time |
---|---|---|---|---|---|
trial_A | 3.0 | 0.0 | 3 | 0 | 0 |
trial_A | 4.5 | 0.0 | 3 | 0 | 0 |
trial_A | 6.0 | 0.0 | 6 | 0 | 0 |
trial_A | 8.0 | 0.0 | 3 | 2 | 0 |
trial_B | 0.0 | 33.3 | 3 | 0 | 0 |
trial_B | 0.0 | 50.0 | 3 | 0 | 0 |
trial_B | 0.0 | 100.0 | 4 | 0 | 0 |
trial_B | 0.0 | 200.0 | 9 | 0 | 0 |
trial_B | 0.0 | 400.0 | 15 | 0 | 0 |
trial_B | 0.0 | 800.0 | 20 | 2 | 0 |
trial_B | 0.0 | 1120.0 | 17 | 4 | 0 |
The objective is to aid dosing and dose-escalation decisions in a future trial, trial_AB
, in which the drugs will be combined. Additionally, another investigator-initiated trial IIT
will study the same combination concurrently. Note that these as-yet-unobserved sources of data are included in the input data as unobserved factor levels. This mechanism allows us to specify a joint meta-analytic prior for all four sources of historical and concurrent data.
## [1] "trial_A" "trial_B" "IIT" "trial_AB"
To fit the hierarchical model described in [4], one makes a call to the function blrm_exnex
, as below; although we slightly deviate from the model in [4] by allowing an EXchangeable/Non-EXchangeable prior for the drug components.
## Load involved packages
library(dplyr) ## for mutate
library(tidyr) ## defines expand_grid
## Design parameters ---------------------
dref <- c(3, 960)
num_comp <- 2 # two investigational drugs
num_inter <- 1 # one drug-drug interaction needs to be modeled
num_groups <- nlevels(hist_combo2$group_id) # four groups of data
num_strata <- 1 # no stratification needed
## abbreviate labels (wait for prior_summary)
options(OncoBayes2.abbreviate.min=8)
## disables abbreviations (default)
#options(OncoBayes2.abbreviate.min=0)
## STRONGLY RECOMMENDED => uses 4 cores
options(mc.cores=4)
## Model fit -----------------------------
blrmfit <- blrm_exnex(
cbind(num_toxicities, num_patients - num_toxicities) ~
1 + I(log(drug_A / dref[1])) |
1 + I(log(drug_B / dref[2])) |
0 + I(drug_A/dref[1] * drug_B/dref[2]) |
group_id,
data = hist_combo2,
prior_EX_mu_mean_comp = matrix(
c(logit(0.1), 0, # hyper-mean of (intercept, log-slope) for drug A
logit(0.1), 0), # hyper-mean of (intercept, log-slope) for drug B
nrow = num_comp,
ncol = 2,
byrow = TRUE
),
prior_EX_mu_sd_comp = matrix(
c(3.33, 1, # hyper-sd of mean mu for (intercept, log-slope) for drug B
3.33, 1), # hyper-sd of mean mu for (intercept, log-slope) for drug B
nrow = num_comp,
ncol = 2,
byrow = TRUE
),
prior_EX_tau_mean_comp = matrix(
c(log(0.25), log(0.125),
log(0.25), log(0.125)),
nrow = num_comp,
ncol = 2,
byrow = TRUE
),
prior_EX_tau_sd_comp = matrix(
c(log(4) / 1.96, log(4) / 1.96,
log(4) / 1.96, log(4) / 1.96),
nrow = num_comp,
ncol = 2,
byrow = TRUE
),
prior_EX_mu_mean_inter = 0,
prior_EX_mu_sd_inter = 1.121,
prior_EX_tau_mean_inter = matrix(log(0.125), nrow = num_inter, ncol = num_strata),
prior_EX_tau_sd_inter = matrix(log(4) / 1.96, nrow = num_inter, ncol = num_strata),
prior_is_EXNEX_comp = rep(TRUE, num_comp),
prior_is_EXNEX_inter = rep(FALSE, num_inter),
## historical data is 100% EX
## new data only 80% EX
prior_EX_prob_comp = matrix(c(1.0, 1.0,
1.0, 1.0,
0.8, 0.8,
0.8, 0.8),
nrow = num_groups,
ncol = num_comp,
byrow=TRUE),
prior_EX_prob_inter = matrix(1, nrow = num_groups, ncol = num_inter),
prior_tau_dist = 1
)
The function blrm_exnex
returns an object from which numerical and graphical posterior summaries can be extracted using OncoBayes2 functions. We recommend making use of the methods described below.
The function prior_summary
provides a facility for printing, in a readable format, a summary of the prior specification.
The main target of inference is generally the probability of DLT at a selection of possible doses. In order to obtain this inference, one needs to specify the covariate levels of interest (which need not be present in the observed data).
In this case, we are interested in predicitons for trial_AB
, with the possible combination doses of drugs A and B below.
newdata <- expand_grid(
group_id = c("trial_AB"),
drug_A = c(0, 3, 4.5, 6),
drug_B = c(0, 400, 600, 800)
)
newdata$group_id <- factor(newdata$group_id, levels(hist_combo2$group_id))
Note it is important that the factor levels associated with newdata$group_id
be consistent with the levels of the grouping factor used when calling blrm_exnex
. In case you use expand.grid
, then the stringsAsFactors = FALSE
argument must be used to ensures that the variable will initially be treated as a character. We recommend using the tidyr
command expand_grid
which avoids casting strings to factors in an uncontrolled manner. The next line converts the group_id
column to a factor with the correct set of four levels. Alternativley, we can create the group_id
column directly as a factor
column within newdata
. The code below results in the same newdata
in a more compact form.
newdata <- expand_grid(
group_id = factor(c("trial_AB"), levels(hist_combo2$group_id)),
drug_A = c(0, 3, 4.5, 6),
drug_B = c(0, 400, 600, 800)
)
Posterior summary statistics for the DLT rates at these provisional doses can be extracted from the blrmfit
object using the summary
method.
summ_stats <- summary(blrmfit,
newdata = newdata,
prob = 0.95,
interval_prob = c(0, 0.16, 0.33, 1))
kable(cbind(newdata, summ_stats), digits = 3)
group_id | drug_A | drug_B | mean | sd | 2.5% | 50% | 97.5% | (0,0.16] | (0.16,0.33] | (0.33,1] |
---|---|---|---|---|---|---|---|---|---|---|
trial_AB | 0.0 | 0 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 1.000 | 0.000 | 0.000 |
trial_AB | 0.0 | 400 | 0.069 | 0.157 | 0.000 | 0.022 | 0.697 | 0.918 | 0.030 | 0.051 |
trial_AB | 0.0 | 600 | 0.100 | 0.169 | 0.001 | 0.051 | 0.779 | 0.880 | 0.056 | 0.064 |
trial_AB | 0.0 | 800 | 0.142 | 0.178 | 0.002 | 0.094 | 0.829 | 0.775 | 0.142 | 0.084 |
trial_AB | 3.0 | 0 | 0.101 | 0.185 | 0.000 | 0.035 | 0.814 | 0.848 | 0.076 | 0.075 |
trial_AB | 3.0 | 400 | 0.167 | 0.227 | 0.004 | 0.078 | 0.939 | 0.713 | 0.146 | 0.140 |
trial_AB | 3.0 | 600 | 0.202 | 0.236 | 0.009 | 0.109 | 0.952 | 0.625 | 0.187 | 0.188 |
trial_AB | 3.0 | 800 | 0.251 | 0.250 | 0.012 | 0.160 | 0.961 | 0.500 | 0.238 | 0.262 |
trial_AB | 4.5 | 0 | 0.147 | 0.207 | 0.002 | 0.073 | 0.913 | 0.746 | 0.150 | 0.104 |
trial_AB | 4.5 | 400 | 0.216 | 0.245 | 0.009 | 0.120 | 0.965 | 0.596 | 0.202 | 0.202 |
trial_AB | 4.5 | 600 | 0.258 | 0.263 | 0.010 | 0.154 | 0.974 | 0.510 | 0.208 | 0.282 |
trial_AB | 4.5 | 800 | 0.310 | 0.288 | 0.009 | 0.204 | 0.982 | 0.437 | 0.193 | 0.370 |
trial_AB | 6.0 | 0 | 0.207 | 0.227 | 0.004 | 0.128 | 0.957 | 0.580 | 0.244 | 0.176 |
trial_AB | 6.0 | 400 | 0.281 | 0.271 | 0.011 | 0.179 | 0.980 | 0.462 | 0.222 | 0.317 |
trial_AB | 6.0 | 600 | 0.325 | 0.298 | 0.008 | 0.217 | 0.985 | 0.421 | 0.190 | 0.389 |
trial_AB | 6.0 | 800 | 0.374 | 0.328 | 0.006 | 0.270 | 0.992 | 0.391 | 0.156 | 0.453 |
Such summaries may be used to assess which combination doses have acceptable risk of toxicity. For example, according the escalation with overdose control (EWOC) design criteria [3], any doses for which the last column does not exceed 25% are eligible for enrollment.
Before trial_AB
is initiated, it may be of interest to test how this model responds in various scenarios for the initial combination cohort(s).
This can be done easily in OncoBayes2 by updating the initial model fit with additional rows of hypothetical data. In the code below, we explore 3 possible outcomes for an initial cohort enrolled at 3 mg Drug A + 400 mg Drug B, and review the model’s inference at adjacent doses. Note that the update
method for accepts an add_data
argument which appends the additional data to the existing data of the previous analysis.
# set up two scenarios at the starting dose level
# store them as data frames in a named list
scenarios <- expand_grid(
group_id = factor("trial_AB", levels(hist_combo2$group_id)),
drug_A = 3,
drug_B = 400,
num_patients = 3,
num_toxicities = 1:2
) %>% split(1:2) %>% setNames(paste(1:2, "DLTs"))
candidate_doses <- expand_grid(
group_id = factor("trial_AB", levels(hist_combo2$group_id)),
drug_A = c(3, 4.5),
drug_B = 400
)
scenario_inference <- lapply(scenarios, function(scenario_newdata){
# refit the model with each scenario's additional data
scenario_fit <- update(blrmfit, add_data = scenario_newdata)
# summarize posterior at candidate doses
scenario_summ <- summary(scenario_fit,
newdata = candidate_doses,
interval_prob = c(0, 0.16, 0.33, 1))
cbind(candidate_doses, scenario_summ)
})
group_id | drug_A | drug_B | mean | sd | 2.5% | 50% | 97.5% | (0,0.16] | (0.16,0.33] | (0.33,1] |
---|---|---|---|---|---|---|---|---|---|---|
trial_AB | 3.0 | 400 | 0.190 | 0.158 | 0.021 | 0.143 | 0.629 | 0.554 | 0.293 | 0.153 |
trial_AB | 4.5 | 400 | 0.264 | 0.199 | 0.030 | 0.206 | 0.783 | 0.382 | 0.324 | 0.294 |
group_id | drug_A | drug_B | mean | sd | 2.5% | 50% | 97.5% | (0,0.16] | (0.16,0.33] | (0.33,1] |
---|---|---|---|---|---|---|---|---|---|---|
trial_AB | 3.0 | 400 | 0.44 | 0.259 | 0.071 | 0.398 | 0.933 | 0.160 | 0.256 | 0.584 |
trial_AB | 4.5 | 400 | 0.52 | 0.265 | 0.090 | 0.508 | 0.970 | 0.094 | 0.200 | 0.706 |
In the example of [1], at the time of completion of trial_AB
, the complete historical and concurrent data are as follows.
group_id | drug_A | drug_B | num_patients | num_toxicities | cohort_time |
---|---|---|---|---|---|
trial_A | 3.0 | 0.0 | 3 | 0 | 0 |
trial_A | 4.5 | 0.0 | 3 | 0 | 0 |
trial_A | 6.0 | 0.0 | 6 | 0 | 0 |
trial_A | 8.0 | 0.0 | 3 | 2 | 0 |
trial_A | 4.5 | 0.0 | 3 | 0 | 1 |
trial_A | 6.0 | 0.0 | 5 | 0 | 1 |
trial_B | 0.0 | 33.3 | 3 | 0 | 0 |
trial_B | 0.0 | 50.0 | 3 | 0 | 0 |
trial_B | 0.0 | 100.0 | 4 | 0 | 0 |
trial_B | 0.0 | 200.0 | 9 | 0 | 0 |
trial_B | 0.0 | 400.0 | 15 | 0 | 0 |
trial_B | 0.0 | 800.0 | 20 | 2 | 0 |
trial_B | 0.0 | 1120.0 | 17 | 4 | 0 |
IIT | 3.0 | 400.0 | 3 | 0 | 1 |
IIT | 3.0 | 800.0 | 7 | 5 | 1 |
IIT | 4.5 | 400.0 | 3 | 0 | 1 |
IIT | 6.0 | 400.0 | 6 | 0 | 1 |
IIT | 6.0 | 600.0 | 3 | 2 | 1 |
trial_AB | 3.0 | 400.0 | 3 | 0 | 1 |
trial_AB | 3.0 | 800.0 | 6 | 2 | 1 |
trial_AB | 4.5 | 600.0 | 10 | 2 | 1 |
trial_AB | 6.0 | 400.0 | 10 | 3 | 1 |
Numerous toxicities were observed in the concurrent IIT
study. Through the MAC framework, these data can influence the model summaries for trial_AB
. Note that we use the update
function differently than before, since we specify the entire data-set now we use the data
argument.
final_fit <- update(blrmfit, data = codata_combo2)
summ <- summary(final_fit, newdata, prob = c(0.5, 0.95), interval_prob = c(0,0.33,1))
final_summ_stats <- cbind(newdata, summ) %>%
mutate(EWOC=1*`(0.33,1]`<=0.25)
ggplot(final_summ_stats,
aes(x=factor(drug_B), colour=EWOC)) +
facet_wrap(~drug_A, labeller=label_both) +
scale_y_continuous(breaks=c(0, 0.16, 0.33, 0.4, 0.6, 0.8, 1.0)) +
coord_cartesian(ylim=c(0,0.8)) +
geom_hline(yintercept = c(0.16, 0.33),
linetype = "dotted") +
geom_pointrange(aes(y=`50%`, ymin=`2.5%`, ymax=`97.5%`)) +
geom_linerange(aes(ymin=`25%`, ymax=`75%`), size=1.5) +
ggtitle("DLT Probability", "Shown is the median (dot), 50% CrI (thick line) and 95% CrI (thin line)") +
ylab(NULL) +
xlab("Dose Drug B [mg]")
[1] Neuenschwander, B., Roychoudhury, S., & Schmidli, H. (2016). On the use of co-data in clinical trials. Statistics in Biopharmaceutical Research, 8(3), 345-354.
[2] Neuenschwander, B., Wandel, S., Roychoudhury, S., & Bailey, S. (2016). Robust exchangeability designs for early phase clinical trials with multiple strata. Pharmaceutical statistics, 15(2), 123-134.
[3] Neuenschwander, B., Branson, M., & Gsponer, T. (2008). Critical aspects of the Bayesian approach to phase I cancer trials. Statistics in medicine, 27(13), 2420-2439.
[4] Neuenschwander, B., Matano, A., Tang, Z., Roychoudhury, S., Wandel, S. Bailey, Stuart. (2014). A Bayesian Industry Approach to Phase I Combination Trials in Oncology. In Statistical methods in drug combination studies (Vol. 69). CRC Press.
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.6 LTS
##
## Matrix products: default
## BLAS: /usr/lib/libblas/libblas.so.3.6.0
## LAPACK: /usr/lib/lapack/liblapack.so.3.6.0
##
## locale:
## [1] C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] tidyr_1.0.0 dplyr_0.8.3 ggplot2_3.2.1 knitr_1.25
## [5] OncoBayes2_0.6-5 Rcpp_1.0.2
##
## loaded via a namespace (and not attached):
## [1] rstan_2.19.3 tidyselect_0.2.5 xfun_0.10
## [4] purrr_0.3.3 colorspace_1.4-1 vctrs_0.2.0
## [7] htmltools_0.4.0 stats4_3.6.1 loo_2.1.0
## [10] yaml_2.2.0 rlang_0.4.0 pkgbuild_1.0.6
## [13] pillar_1.4.2 glue_1.3.1 withr_2.1.2
## [16] matrixStats_0.55.0 lifecycle_0.1.0 plyr_1.8.4
## [19] stringr_1.4.0 munsell_0.5.0 gtable_0.3.0
## [22] codetools_0.2-16 evaluate_0.14 inline_0.3.15
## [25] callr_3.3.2 ps_1.3.0 parallel_3.6.1
## [28] bayesplot_1.7.0 rstantools_2.0.0 highr_0.8
## [31] backports_1.1.5 scales_1.0.0 checkmate_1.9.4
## [34] StanHeaders_2.19.2 abind_1.4-5 gridExtra_2.3
## [37] digest_0.6.21 stringi_1.4.3 processx_3.4.1
## [40] grid_3.6.1 cli_1.1.0 tools_3.6.1
## [43] magrittr_1.5 lazyeval_0.2.2 tibble_2.1.3
## [46] Formula_1.2-3 crayon_1.3.4 pkgconfig_2.0.3
## [49] zeallot_0.1.0 ellipsis_0.3.0 prettyunits_1.0.2
## [52] ggridges_0.5.1 assertthat_0.2.1 rmarkdown_1.16
## [55] R6_2.4.0 compiler_3.6.1