---
title: "Bayesian Optimal Interval Design"
output: 
  rmarkdown::html_vignette:
    df_print: tibble
bibliography: library.bib
vignette: >
  %\VignetteIndexEntry{Bayesian Optimal Interval Design}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7, fig.height = 5
)
```

# Introduction
The Bayesian Optimal Interval (BOIN) design was introduced by @liu_bayesian_2015.
It is one of a series of dose-finding trial designs that works by partitioning the probability of toxicity into a set of intervals.
These designs make dose-selection decisions that are determined by the interval in which the probability of toxicity for the current dose is believed to reside.


# Summary of the BOIN Design
BOIN seeks a dose with probability of toxicity close to some pre-specified target level, $p_T$.
Let the probability of toxicity at dose $i$ be $p_i$.
The entire range of possible values for $p_i$ can be broken up into the following intervals: 

* The underdosing interval (UI), defined as $(0, \lambda_{1i})$;
* The equivalence interval (EI), defined as $(\lambda_{1i}, \lambda_{2i})$;
* The overdosing interval (OI), defined as $(\lambda_{2i}, 1)$;

Let $\pi_{UI,i}, \pi_{EI,i}$ and $\pi_{OI,i}$ be the a-priori probabilities that the rate of toxicity associated with dose $i$ belongs to the intervals UI, EI and OI.
By definition, $\pi_{UI,i} + \pi_{EI,i} + \pi_{OI,i} = 1$.
The authors advocate $\pi_{UI,i} = \pi_{EI,i} = \pi_{OI,i} = \frac{1}{3}$.

Let $n_i$ be the number of patients that have been treated at dose $i$, yielding $x_i$ toxicity events.
The so-called _local_ BOIN variant (i.e. that advocated by @liu_bayesian_2015) defines: 

$$ \lambda_{1i} = \frac{\log{\left( \frac{1 - \phi_1}{1 - p_T} \right)} + \frac{1}{n_i} \log{\left( \frac{\pi_{UI,i}}{\pi_{EI,i}} \right)} }{ \log{\left( \frac{p_T (1 - \phi_1)}{\phi_1 (1 - p_T)} \right)} }$$

$$ \lambda_{2i} = \frac{\log{\left( \frac{1 - p_T}{1 - \phi_2} \right)} + \frac{1}{n_i} \log{\left( \frac{\pi_{EI,i}}{\pi_{OI,i}} \right)} }{ \log{\left( \frac{ \phi_2(1 - p_T)}{p_T (1 - \phi_2)} \right)} }$$

where $\phi_1$ and $\phi_2$ are model parameters.
The authors advocate $\phi_1 \in \left[ 0.5pT, 0.7pT \right]$ and $\phi_2 \in \left[ 1.3pT, 1.5pT \right]$.
As defaults, they recommend $\phi_1 = 0.6p_T$ and $\phi_2 = 1.4p_T$.

Having observed toxicity rate $\hat{p}_i = x_i / n_i$ at the current dose $i$, the logical action depends on the interval in which $\hat{p}_i$ resides. 
If $\hat{p}_i < \lambda_{1i}$, then the current dose is likely an underdose, so our desire should be to escalate dose to $i+1$.
In contrast, if $\hat{p}_i > \lambda_{2i}$, then the current dose is likely an overdose and we will want to de-escalate dose to $i-1$ for the next patient.
If $\lambda_{1i} < \hat{p}_i < \lambda_{2i}$, then the current dose is deemed sufficiently close to $p_T$ and we will want to stay at dose-level $i$.

The authors advocate a stopping rule to protect against repeated administration of a dose that is evidently excessively toxic.
The proposed rule is similar to that used in TPI and mTPI.
Using a $Beta(1, 1)$ prior, dose $i$ is deemed inadmissible for being excessively toxic if

$$ Pr(p_{i} > p_{T} | x_i, n_i) > \xi,$$

for a certainty threshold, $\xi$, with $\xi = 0.95$ being suggested.
If a dose is excluded by this rule, it should not be recommended by the model.
Irrespective the values of $\lambda_{1i}$ and $\lambda_{1i}$, the design will recommend to stay at dose $i$ rather than escalate to a dose previously identified as being inadmissible.
Furthermore, the design will advocate stopping if the lowest dose is inferred to be inadmissible.

See @liu_bayesian_2015 and @BOIN for full details.



# Implementation in `escalation`
To demonstrate the method, let us reproduce the dose selection sequence in a trial of five doses targeting $p_T = 0.3$, described in @BOIN. 

Opting to take the defaults, $\phi_1 = 0.6 p_T = 0.18$ and $\phi_2 = 1.4 p_T = 0.42$, we create an object to fit the model using:
```{r, message=FALSE}
library(escalation)

model <- get_boin(num_doses = 5, target = 0.3)
```

This is short-hand for 
```{r}
model <- get_boin(num_doses = 5, target = 0.3, p.saf = 0.18, p.tox = 0.42)
```

The text in the paper describes that outcomes '1NNN' were observed in the first cohort:
```{r}
fit <- model %>% fit('1NNN')
```

leading to advice to escalate:
```{r}
fit %>% recommended_dose()
```

The next cohort also saw three non-toxicity events, leading to advice:

```{r}
fit <- model %>% fit('1NNN 2NNN')
fit %>% recommended_dose()
```

In the third cohort, two patients had toxicity, leading to advice:

```{r}
fit <- model %>% fit('1NNN 2NNN 3NTT')
fit %>% recommended_dose()
```

The relatively low sample size at dose 3 means that the dose has not yet been rendered inadmissible:

```{r}
fit %>% dose_admissible()
```

## Dose paths
We can reproduce some of the advice above using dose-paths to exhaustively calculate all possible future model advice.
For example, after observing `1NNN` in the first cohort, we can reproduce the advice to escalate further after seeing `2NNN` and then de-escalate after seeing `3NTT` using:

```{r}
cohort_sizes <- c(3, 3)
paths <- model %>% get_dose_paths(cohort_sizes = cohort_sizes, 
                                  previous_outcomes = '1NNN', next_dose = 2)
graph_paths(paths)
```

Thus, we can trace the path `1NNN 2NNN 3NNT`, along with every other possible path in two cohorts of three after having observed `1NNN`.

For more information on working with dose-paths, refer to the dose-paths vignette.


## Simulation
@liu_bayesian_2015 present in their Table 4 some simulated operating characteristics.
We can use the `simulate_trials` function to reproduce the findings.

Their example concerns a clinical trial of six doses that targets 25% toxicity.
We specify the `model` object to reflect this.
They also elect to limit the trial to a sample size of $n=36$:

```{r}
model <- get_boin(num_doses = 6, target = 0.25) %>% 
  stop_at_n(n = 36)
```

Their scenario 1 investigates the true probability vector:
```{r}
true_prob_tox <- c(0.25, 0.35, 0.5, 0.6, 0.7, 0.8)
```

For the sake of speed, we will run just fifty iterations:

```{r}
num_sims <- 50
```

In real life, however, we would naturally run many thousands of iterations.

Running the simulation:

```{r}  
sims <- model %>% 
  simulate_trials(num_sims = num_sims, true_prob_tox = true_prob_tox)
```

we see that from this small sample size, the probability that each dose is recommended is similar to that reported:
```{r}
prob_recommend(sims)
```

and so is the expected number of patients treated at each dose level:
```{r}
colMeans(n_at_dose(sims))
```

For more information on running dose-finding simulations, refer to the simulation vignette.


# References
