---
title: "Trial simulation with binary toxicity measure"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Trial simulation with binary toxicity measure}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
```

## Overview 

For illustrative purposes, we start by showing how to simulate a single trial, followed by simulations for repeated trials (100 sims).
See the README file and the full [publication](https://www.tandfonline.com/doi/abs/10.1080/19466315.2018.1462727) for greater overview of the design framework. 

The examples we will be working with are the following: suppose we have a study with 5 dose levels to test, each of which has an associated "true" toxicity level and "true" mean efficacy (used for data generation). Additionally, we specify an acceptable and unacceptable rate of dose-limiting toxicities (DLTs). 
In stage 1, equal cohort sizes are assigned to each dose level (default is 3). In stage 2, we also assign equal cohort sizes to only acceptable doses (default changes to 1).
We employ a continuous efficacy outcome and a binary toxicity endpoint. The design assumes a monotonically increasing dose-toxicity relationship (i.e, as the dose increases, so does toxicity), but no monotonicity is imposed for the dose-efficacy relationship (accommodates increasing, flat, plateau trends).

In scenario 1, we will work with a linear increasing dose-efficacy relationship; in scenario 2, we show simulations for a quadratic dose-efficacy relationship. 

## Scenario 1: Monotone dose-efficacy

```{r}
library(iAdapt)

# Number of pre-specified dose levels
dose <- 5 

# Vector of true toxicities associated with each dose
dose.tox <- c(0.05, 0.10, 0.20, 0.35, 0.45)       

# Acceptable (p_yes) and unacceptable (p_no) DLT rates used for establishing safety
p_no <- 0.40                                     
p_yes <- 0.15    

# Likelihood-ratio (LR) threshold
K <- 2                                          

# Cohort size used in stage 1
coh.size <- 3 

# Vector of true mean efficacies per dose (here mean T-cell persistence per dose (%))
m <- c(5, 15, 40, 65, 80)   # MUST BE THE SAME LENGTH AS dose.tox                  

# Efficacy (equal) variance per dose
v <- rep(0.01, 5) 

# Total sample size (stages 1&2)                            
N <- 25                                        

# Stopping rule: if dose 1 is the only safe dose, allocate up to 9 pts before ending the trial to collect more information
stop.rule <- 9   
```

Note that, as we've specified them, we have monotone increasing dose-toxicity and dose-efficacy curves, as shown below. If the unacceptable dose-limiting toxicity (DLT) rate is 0.40, (horizontal dotted line, figure a), then our target/optimal dose is dose 4 (green point) - the safe dose with the highest efficacy. 


```{r, echo=FALSE, fig.width=6, fig.height=3}
par(mfrow = c(1,2))
plot(x = 1:5, y = dose.tox, ylim = c(0, p_no + 0.1),
     type = 'b',
     xlab = "Dose level", ylab = "True toxicity rate",
     main = "(a) Dose-toxicity")
lines(x = 0:6, y = rep(p_no, 7), lty = 3) # threshold toxicity value
points(x = 4, y = dose.tox[4], pch = 15, col = "green") # highlight best dose

plot(x = 1:5, y = m, ylim = c(0, 100),
     type = 'b',
     xlab = "Dose level", ylab = "True efficacy (T-cell % persistence)",
     main = "(b) Dose-efficacy")
points(x = 4, y = m[4], pch = 15, col = "green") # highlight best dose

dev.off()
```

### One Trial Simulation
#### Stage 1: Establish the safety profile for all initial doses. 

Stage 1 establishes the safety profiles of the predefined doses. 

Function to generate and tabulate toxicities per dose level:
```{r}
set.seed(3)
tox.profile(dose = dose, dose.tox = dose.tox, p1 = p_no, p2 = p_yes, K = K, coh.size = coh.size)
```

The first column indicates the cohort number and third columns gives the dose assignment for all specified doses (in this case, 5). The second column gives the number of DLTs observed at that dose. The fourth column gives the likelihood ratio calculated from the observed DLTs. 
A dose is considered acceptably safe if the LR > 1/K and unacceptably safe if LR <= 1/K.

Now let's see which doses the design selects as being acceptably safe using K=2. 

Function to select only the acceptable toxic doses: 
```{r}
set.seed(3)
safe.dose(dose = dose, dose.tox = dose.tox, p1 = p_no, p2 = p_yes, K = K, coh.size = coh.size) 
```
 
We can see that the design selects only doses 1 through 4; \$alloc.safe gives the dose assignment (first column) and number of DLTs (second column). \$alloc.total gives the dose assignment for all enrolled patients (\$n1 gives the total sample size used in stage 1, in this case 15 subjects), where we see that 3 patients were assigned to each dose as specified by `coh.size`.

Stage 1 is mainly used to establish safety, but efficacy outcomes are also collected for each dose.
 
Function to generate efficacy outcomes (here T-cell percent persistence) for each dose:
```{r}
set.seed(3)
eff.stg1(dose = dose, dose.tox = dose.tox, p1 = p_no, p2 = p_yes, K = K, coh.size = coh.size, m, v, nbb = 100)
```

\$Y.safe and \$d.safe give the efficacy values and dose allocations for all subjects enrolled at acceptably safe doses; \$tox.safe gives the number of DLTs for each dose level; \$Y.alloc and \$d.alloc gives the efficacy values and dose allocations for all subjects enrolled in stage 1 (safe and unsafe doses). Notice that the \$Y.safe and \$d.safe are subsets of \$Y.alloc and \$d.alloc.

### Stage 2: Adaptive randomization based on efficacy outcomes.

If 2 or more doses are considered acceptable after stage 1, the remaining patients are randomized to these open doses until the total sample size N is reached. If only dose 1 is acceptable after stage 1, allocate up to 9 patients (`stop.rule = 9`). Toxicity is still being monitored (in the 'background') throughout stage 2, so acceptable doses (declared in stage 1) can still be discarded based on observed DLTs. The discarded dose and all levels above it cannot be revisited.

Function to fit a linear regression for the continuous efficacy outcomes, compute the randomization probabilities per dose and allocate the next subject to an acceptable safe dose that has the highest randomization probability:
```{r}
set.seed(3)
rand.stg2(dose, dose.tox, p_no, p_yes, K, coh.size, m, v, N, stop.rule = stop.rule, cohort = 1, samedose = TRUE, nbb = 100) 
```

Notice that after stage 1 (dose 5, efficacy value 85), patients were randomized only to dose 4 or lower (acceptably safe).
The complete vectors of dose allocations and efficacy outcomes for the entire trial (N=25) can be used to compute the operating characteristics of the design in repeated simulations.
 
### 100 Trials Simulations 

To simulate this trial 100 times, we can run the following:
```{r, results='hide'}
sims = 1e2

set.seed(3)
simulations = sim.trials(numsims = sims, dose, dose.tox, p1 = p_no, p2 = p_yes, K, coh.size, m, v, N, stop.rule = stop.rule, cohort = 1, samedose = TRUE, nbb = 100)
```

\$safe.d indicates whether each dose (column) was declared safe in stage 1 (1 = yes, 0 = no) for each trial (row).  

```{r}
head(simulations$safe.d)
head(simulations$sim.Y)
head(simulations$sim.d)
```

\$sim.Y gives the observed outcomes, where each column corresponds to 1 patient (maximum of N columns), and each row is a simulated trial. Correspondingly, \$sim.d gives the dose allocation for each patient (column) in each trial (row). 

To see the proportion of times we've designated each dose as safe in stage 1, we can simply take the column totals of the \$safe.d matrix:
```{r}
colSums(simulations$safe.d) / sims
```

Simulation results can be summarized. For each dose level, the inter-quartile range (25th percentile, median, 75th percentile) for the percent of subjects treated and observed efficacy are given in tables.

```{r}
sim.summary(simulations)
```

 
## Scenario 2: Non-monotone dose-efficacy

Here, all parameters are the same as in scenario 1, except for the true mean efficacies. Instead, they take on a parabolic nature.

```{r}
# Vector of true mean efficacies per dose (here mean percent persistence per dose)
m <- c(15, 35, 80, 60, 40)   # MUST BE THE SAME LENGTH AS dose.tox            
```

Now we have a non-monotone dose-efficacy curves, as shown below. If we have an acceptable dose-limiting toxicity (DLT) rate of 0.40, (horizontal dotted line, figure a), then our target/optimal dose would be dose 4 as we did in the previous scenario (red point). However, note that according to the true efficacies of each dose, the lesser dose number 3 is actually best (green point). Assuming a monotone relationship, we would be incorrectly pushing forward a more toxic and less effective dose level. 

```{r, echo=F, fig.width=6, fig.height=3}
par(mfrow = c(1,2))
plot(x = 1:5, y = dose.tox, ylim = c(0, p_no + 0.1),
     type = 'b',
     xlab = "Dose level", ylab = "True toxicity rate",
     main = "(a) Dose-toxicity")
lines(x = 0:6, y = rep(p_no, 7), lty = 3) # threshold toxicity value
points(x = 3, y = dose.tox[3], pch = 15, col = "green") # highlight most effecive dose
points(x = 4, y = dose.tox[4], pch = 15, col = "red") # highlight MTD

plot(x = 1:5, y = m, ylim = c(0, 100),
     type = 'b',
     xlab = "Dose level", ylab = "True efficacy (T-cell % persistence)",
     main = "(b) Dose-efficacy")
points(x = 3, y = m[3], pch = 15, col = "green") # highlight best dose
points(x = 4, y = m[4], pch = 15, col = "red") # highlight MTD

dev.off()
```

### One Trial Simulation
#### Stage 1: Establish the safety profile for all initial doses

Stage 1 establishes the safety profiles of the predefined doses. All output can be interpreted in the same way as in scenario 1. However, note that in this case we identify dose level 3 as the last acceptably safe dose.


Function to generate and tabulate toxicities per dose level.
```{r}
set.seed(1)
tox.profile(dose = dose, dose.tox = dose.tox, p1 = p_no, p2 = p_yes, K = K, coh.size = coh.size)
```
Dose level 4 had a LR=0.20 < 1/2; it was considered unsafe, so allocation stopped at this level (no dose 5 was assigned).

Function to select only the acceptable safe doses 
```{r}
set.seed(1)
safe.dose(dose = dose, dose.tox = dose.tox, p1 = p_no, p2 = p_yes, K = K, coh.size = coh.size) 
```
 
Function to generate efficacy outcomes (here percent persistence) for each dose
```{r}
set.seed(1)
eff.stg1(dose = dose, dose.tox = dose.tox, p1 = p_no, p2 = p_yes, K = K, coh.size = coh.size, m, v, nbb = 100)
```

 

#### Stage 2: Adaptive randomization based on efficacy outcomes.

If 2 or more doses are considered acceptable after stage 1, the remaining patients are randomized to these open doses until the total sample size N is reached. If only dose 1 is acceptable after stage 1, allocate up to 9 patients (stop.rule=9)	

Function to fit a linear regression for the continuous efficacy outcomes, compute the randomization probabilities per dose and allocate the next subject to an acceptable safe dose that has the highest randomization probability.

```{r}
set.seed(1)
rand.stg2(dose, dose.tox, p_no, p_yes, K, coh.size, m, v, N, stop.rule = stop.rule, cohort = 1, samedose = TRUE, nbb = 100) 
```

Notice that after stage 1 (dose 4, efficacy value 79), patients were randomized to only to dose 3 orlower (acceptably safe).
The complete vectors of dose allocations and efficacy outcomes for the entire trial (N=25) can be used to compute the operating characteristics of the design in repeated simulations.
The complete vectors of dose allocations and efficacy outcomes can be used to compute the operating characteristics of the design in repeated simulations.
 
 
 To simulate this trial 100 times, we can run the following:
```{r, results='hide'}
sims = 1e2

set.seed(1)
simulations = sim.trials(numsims = sims, dose, dose.tox, p1 = p_no, p2 = p_yes, K, coh.size, m, v, N, stop.rule = stop.rule, cohort = 1, samedose = TRUE, nbb = 100)
```

\$safe.d indicates whether each dose (column) was declared safe in stage 1 (1 = yes, 0 = no) for each trial (row).  

```{r}
head(simulations$safe.d)
head(simulations$sim.Y)
head(simulations$sim.d)
```

\$sim.Y gives the observed outcomes, where each column corresponds to 1 patient (maximum of N columns), and each row is a simulated trial. Correspondingly, \$sim.d gives the dose allocation for each patient (column) in each trial (row). 

To see the rate at which we've designated each dose as safe in stage 1, we can simply take the column totals of the \$safe.d matrix:
```{r}
colSums(simulations$safe.d) / sims
```

Simulation results are also summarized below.

```{r}
sim.summary(simulations)
```