# Guiding Oncology Dose-Escalation Trials

#### Andrew Bean, Sebastian Weber

#### 2020-05-07

## Introduction

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.

## Example use-case

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.

    kable(hist_combo2)

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.

    levels(hist_combo2$group_id)

    ## [1] "trial_A"  "trial_B"  "IIT"      "trial_AB"

### Fitting the model

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.

### Summary of prior specification

The function `prior_summary` provides a facility for printing, in a readable format, a summary of the prior specification.

    prior_summary(blrmfit) # not run here

### Summary of posterior

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.

### Data scenarios

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)

    })

Model inference when 1 DLT is observed in first cohort
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

Model inference when 2 DLTs are observed in first cohort
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

## Continuation of example

In the example of [1], at the time of completion of `trial_AB`, the complete historical and concurrent data are as follows.

    kable(codata_combo2)

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]")

## References

[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.

## Session Info

    sessionInfo()

    ## 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