---
title: "Guiding Oncology Dose-Escalation Trials"
author: "Andrew Bean, Sebastian Weber"
date: "`r Sys.Date()`"
output: html_vignette
---
<!--
%\VignetteEngine{knitr::rmarkdown}
%\VignetteIndexEntry{Guiding Oncology Dose-Escalation Trials}
-->

```{r, include=FALSE}
library(OncoBayes2)
library(knitr)
library(ggplot2)

theme_set(theme_bw())

knitr::opts_chunk$set(
  fig.width = 1.62*4,
  fig.height = 4
)
## setup up fast sampling when run on CRAN
is_CRAN <- !identical(Sys.getenv("NOT_CRAN"), "true")
## NOTE: for running this vignette locally, please uncomment the
## following line:
## is_CRAN <- TRUE
.user_mc_options <- list()
if (is_CRAN) {
  .user_mc_options <- options(OncoBayes2.MC.warmup=10, OncoBayes2.MC.iter=20, OncoBayes2.MC.chains=1)
}
```

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

```{r}
kable(hist_combo2)
```

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.

```{r}
levels(hist_combo2$group_id)
```

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

```{r, eval = FALSE}

## Load involved packages
library(RBesT)  ## defines logit function
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
)
```


```{r, include = FALSE}
library(RBesT)  ## defines logit function
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.

```{r, eval = FALSE}
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.

```{r}
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.

```{r}
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.

```{r}
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)
```

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.

```{r, eval = FALSE}
# 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)
  
})

```

```{r, include = FALSE}
# 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)
})

```


```{r, echo = FALSE}
kable(scenario_inference[["1 DLTs"]], digits = 3,
      caption = "Model inference when 1 DLT is observed in first cohort")
```

```{r, echo = FALSE}
kable(scenario_inference[["2 DLTs"]], digits = 3,
      caption = "Model inference when 2 DLTs are observed in first cohort")
```

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

```{r}
kable(codata_combo2)
```

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.

```{r, include = FALSE}
final_fit <- update(blrmfit, data = codata_combo2)

## we could equally obtain the same analysis with 
## final_fit <- update(blrmfit, newdata = codata_combo2[12:20,])

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)
```

```{r, eval = FALSE}
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)
```


```{r, fig.height = 1.05 * 4, fig.width=1.62 * 4}

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

```{r}
sessionInfo()
```


```{r, include=FALSE}
## restore previous global user options
options(.user_mc_options)
```
