---
title: "An Example of Simulating a Trial with Adaptive Design"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{An Example of Simulating a Trial with Adaptive Design}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  markdown: 
    wrap: 72
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  cache.path = 'cache/adaptiveDesign/',
  comment = '#>',
  dpi = 300,
  out.width = '100%'
)
```

```{r setup, echo = FALSE, message = FALSE}
library(dplyr)
library(R6)
library(kableExtra)
library(rlang)
library(survival)
library(survminer)
library(ggplot2)
library(TrialSimulator)
```

In this vignette, we illustrate how to use `TrialSimulator` to simulate
a trial with seamless adaptive design of multiple endpoints, where dose
selection, futility analysis, interim analysis are triggered by
user-defined milestones. Family-wise error rate are controlled by closed
test with Dunnett's test, with $\alpha$ split for the endpoints. Please
refer to another vignette where graphical testing procedure is adopted
instead.

## Simulation Settings

-   Trial consists of two active arms of high or low dose, and a placebo
    arm.
-   Patients are randomized into the trial with ratio `1:1:1`.
-   30 patients are randomized per month for the first 10 months, and 50
    patients per month after then until 1000 patients.
-   Dropout rate is 8% at month 12, and 18% at month 18, which is
    modeled by a Weibull distribution. The scale and shape parameters of
    the Weibull distribution can be computed using helper function

```{r daoif, class.source='fold-show'}
weibullDropout(time = c(12, 18), dropout_rate = c(.08, .18))
```

-   Two time-to-event endpoints (`PFS` and `OS`) and a binary
    `surrogate` are under consideration, where a higher rate of
    `surrogate` is better.
-   `surrogate` has a readout time of 5 weeks.
-   `surrogate` is used for dose selection using a hypothetical rule.
    Let $z_h$ (or $z_l$) be the z statistics of `surrogate` between high
    (or low) dose and placebo.
    -   low dose arm is selected if $z_l > 1.28$;
    -   high dose arm is selected if $z_h > 1.28 > z_l$;
    -   keep both arms otherwise.
-   Rates of `surrogate` are
    -   0.05 in placebo
    -   0.12 in low dose arm
    -   0.13 in high dose arm
-   Medians of PFS are
    -   5 months in placebo
    -   6.7 months in low dose arm (hazard ratio 0.75)
    -   7.1 months in high dose arm (hazard ratio 0.70) Exponential
        distribution is assumed.
-   Medians of OS are
    -   14 months in placebo
    -   17.5 months in low dose arm (hazard ratio 0.80)
    -   18.2 months high dose arms (hazard ratio 0.77)
-   Dose selection is done when 300 patients are randomized with readout
    of `surrogate` are ready for analysis.
-   An interim is set to test `PFS` when 300 `PFS` events are observed.
    A non-binding futility analysis is also done with a boundary 0.5.
-   Final analysis is set for `PFS` and `OS` when all planned 1000
    patients are randomized and at least 300 `OS` events are observed.
    Meanwhile, the trial has reached at least 28 months or at least 520
    events are observed for `PFS`.
-   To control the family-wise error rate accounting for `PFS` and `OS`,
    as well as adaptation of dose-selection, we split $\alpha$ between
    `PFS` (0.5%) and `OS` (2%), and adopt closed test with Dunnett's
    test for intersection hypotheses.

## Define Three Arms

In this example, we call `endpoint()` twice for each arm to define `PFS` 
and `OS`. This implicitly assumes independence between the two time-to-event 
endpoints. `TrialSimulator` offers a built-in generator `CorrelatedPfsAndOs3()` 
to generate `PFS` and `OS` with correlation based on the ill-death model, 
which enjoys several nice properties. For selection of its arguments (i.e. 
transition hazards), refer to vignette 
[Simulate Correlated Progression-Free Survival and Overall Survival as Endpoints in Clinical Trials](simulatePfsAndOs.html). 

```{r ljgai, class.source="fold-show"}
#' define three arms
pbo <- arm(name = 'placebo')
low <- arm(name = 'low dose')
high <- arm(name = 'high dose')

#' define endpoints in placebo
pfs <- endpoint(name = 'pfs', type = 'tte',
                generator = rexp, rate = log(2) / 5)

os <- endpoint(name = 'os', type = 'tte',
               generator = rexp, rate = log(2) / 14)

five_weeks <- 5 / 52 * 12 ## convert it in months
surrogate <- endpoint(name = 'surrogate', type = 'non-tte',
                      readout = c(surrogate = five_weeks),
                      generator = rbinom, size = 1, prob = .05)
pbo$add_endpoints(pfs, os, surrogate)

#' define endpoints in low dose arm
pfs <- endpoint(name = 'pfs', type = 'tte',
                generator = rexp, rate = log(2) / 6.7)

os <- endpoint(name = 'os', type = 'tte',
               generator = rexp, rate = log(2) / 17.5)

surrogate <- endpoint(name = 'surrogate', type = 'non-tte',
                      readout = c(surrogate = five_weeks),
                      generator = rbinom, size = 1, prob = .12)
low$add_endpoints(pfs, os, surrogate)

#' define endpoints in high dose arm
pfs <- endpoint(name = 'pfs', type = 'tte',
                generator = rexp, rate = log(2) / 7.1)

os <- endpoint(name = 'os', type = 'tte',
               generator = rexp, rate = log(2) / 18.2)

surrogate <- endpoint(name = 'surrogate', type = 'non-tte',
                      readout = c(surrogate = five_weeks),
                      generator = rbinom, size = 1, prob = .13)
high$add_endpoints(pfs, os, surrogate)
```

## Define a Trial

With three arms, we can define a trial of class `Trial`. Recruitment
curve are specified through `enroller` with a built-in function
`StaggeredRecruiter` of piecewise constant rate. We set `duration` to be
an arbitrary large number (50) but controlling the end of trial through
pre-defined milestones later. Note that if `seed = NULL`, `TrialSimulator`
will pick a seed for the purpose of reproducibility.

```{r lagieg, class.source="fold-show"}
accrual_rate <- data.frame(end_time = c(10, Inf),
                           piecewise_rate = c(30, 50))
trial <- trial(
  name = 'Trial-3415', n_patients = 1000,
  seed = 1727811904, duration = 40,
  enroller = StaggeredRecruiter, accrual_rate = accrual_rate,
  dropout = rweibull, shape = 2.139, scale = 38.343
)

trial$add_arms(sample_ratio = c(1, 1, 1), low, high, pbo)
```

## Define Trial Milestones and Action Functions

Next, we define action functions for trial milestones, i.e., dose selection,
interim and final analysis. Note that an action function should always
has arguments `trial` of class `Trial` and a character `milestone_name`. We
calculate z statistics of the Farrington-Manning test for binary
endpoint `surrogate` using helper function `fitFarringtonManning`, which
returns a data frame (see `?fitFarringtonManning`).

```{r ieaong, class.source='fold-show'}
action1 <- function(trial, milestone_name){
  
  locked_data <- trial$get_locked_data(milestone_name)
  
  fit <- fitFarringtonManning(endpoint = 'surrogate', placebo = 'placebo',
                              data = locked_data, alternative = 'greater')
  
  # browser() ## if you want to see what does fit look like
  z_l <- fit$z[fit$arm == 'low dose']
  z_h <- fit$z[fit$arm == 'high dose']
  if(z_l > 1.28){
    trial$remove_arms('high dose')
    trial$save(value = 'low', name = 'kept_arm')
  }else if(z_h > 1.28){
    trial$remove_arms('low dose')
    trial$save(value = 'high', name = 'kept_arm')
  }else{
    trial$save(value = 'both', name = 'kept_arm')
  }

  invisible(NULL)

}
```

At the interim, we test `PFS` using the logrank test to carry out a
non-binding futility analysis in the following action function. Depending
on how many arms are carried forward in dose selection, `fit` may
consists of testing results of one or two dose arms. Note that a formal
test for `PFS` will be done at the end of the trial when we have all
p-values. If no futility analysis is planned, we can simply set
`action = doNothing` in `milestone`.

```{r aotel, class.source="fold-show"}
action2 <- function(trial, milestone_name){
  
  locked_data <- trial$get_locked_data(milestone_name)

  fit <- fitLogrank(Surv(pfs, pfs_event) ~ arm, placebo = 'placebo', 
                    data = locked_data, alternative = 'less')

  ## futility analysis
  if(max(fit$z) < .5){
    trial$save(value = 'negative', name = 'futility')
    
    ## extend duration
    ## trial$set_duration(45)
  }else{
    trial$save(value = 'positive', name = 'futility')
  }

  invisible(NULL)

}
```

At final analysis, we conduct a closed test using Dunnett's test for
intersection hypotheses. The action function is as follows.

```{r alkdae, class.source="fold-show"}
action3 <- function(trial, milestone_name){
  
  locked_data <- trial$get_locked_data(milestone_name)

  ## test PFS
  dt_pfs <- trial$dunnettTest(Surv(pfs, pfs_event) ~ arm, placebo = 'placebo',
                              treatments = c('high dose', 'low dose'),
                              milestones = c('dose selection', 'interim', 'final'),
                              alternative = 'less',
                              planned_info = 'default')
  
  ct_pfs <- trial$closedTest(dt_pfs, treatments = c('high dose', 'low dose'),
                             milestones = c('interim', 'final'),
                             alpha = .005, alpha_spending = 'asOF')

  ## test OS
  dt_os <- trial$dunnettTest(Surv(os, os_event) ~ arm, placebo = 'placebo',
                             treatments = c('high dose', 'low dose'),
                             milestones = c('dose selection', 'final'),
                             alternative = 'less', 
                             planned_info = 'default')
  
  ct_os <- trial$closedTest(dt_pfs, treatments = c('high dose', 'low dose'),
                            milestones = c('final'),
                            alpha = .02, alpha_spending = 'asOF')
  
  ## we only save testing decision here
  ## You can save whatever you want for summarizing things later, e.g. reject time
  trial$save(value = ct_pfs$decision[ct_pfs$arm == 'high dose'], 
             name = 'pfs_high_dose_decision')
  
  trial$save(value = ct_pfs$decision[ct_pfs$arm == 'low dose'], 
             name = 'pfs_low_dose_decision')
  
  trial$save(value = ct_os$decision[ct_os$arm == 'high dose'], 
             name = 'os_high_dose_decision')
  
  trial$save(value = ct_os$decision[ct_os$arm == 'low dose'], 
             name = 'os_low_dose_decision')

  invisible(NULL)
}
```

Next, we register three trial milestones to a listener

```{r}
dose_selection <- milestone(name = 'dose selection', action = action1,   
                            when = eventNumber(endpoint = 'surrogate', n = 300)
                            )

interim <- milestone(name = 'interim', action = action2, 
                     when = eventNumber(endpoint = 'pfs', n = 300)
                     )

final <- milestone(name = 'final', action = action3, 
                   when = enrollment(n = 1000, arms = c('placebo', 'low dose', 'high dose')) & 
                     eventNumber(endpoint = 'os', n = 300) & (
                       calendarTime(time = 28) | 
                         eventNumber(endpoint = 'pfs', n = 520)
                       )
                   )

listener <- listener()
#' register milestones with listener
listener$add_milestones(
  dose_selection,
  interim,
  final
)
```

## Execute a Trial

We can run a trial as follows. By default, `TrialSimulator` generates a
plot showing cumulative number of patients or endpoint events over time
for each of the arms. We can set it to `FALSE` if a massive number of
simulation replicates is run to save time.

```{r eiaaf, dpi = 1200}
controller <- controller(trial, listener)
controller$run(plot_event = TRUE)
```

In custom action functions, we can use `Trial$save()` to save
intermediate results for summary purpose, which can be accessed anytime
and anywhere by

```{r alojfqoitl, class.source="fold-show"}
controller$get_output() %>% 
  kable(escape = FALSE) %>% 
  kable_styling(bootstrap_options = "striped", 
                full_width = FALSE,
                position = "left") %>%
  scroll_box(width = "100%")
```

Here we dive into the action function (`action3`) for the final
analysis. We can literally execute the function line by line with locked
data loaded. Note that the member functions `Trial$dunnettTest` and
`Trial$closedTest` can make use of locked data at multiple stages
(through the `milestones` argument) automatically, thus, unlike in
`action3`, an explicit call of `Trial$get_locked_data(milestone_name)` is
unnecessary.

```{r aief, class.source='fold-show'}
## test PFS
dt_pfs <- trial$dunnettTest(Surv(pfs, pfs_event) ~ arm, placebo = 'placebo',
                            treatments = c('high dose', 'low dose'),
                            milestones = c('dose selection', 'interim', 'final'),
                            alternative = 'less', 
                            planned_info = 'default')
ct_pfs <- trial$closedTest(dt_pfs, treatments = c('high dose', 'low dose'),
                           milestones = c('interim', 'final'),
                           alpha = .005, alpha_spending = 'asOF')

## test OS
dt_os <- trial$dunnettTest(Surv(os, os_event) ~ arm, placebo = 'placebo',
                           treatments = c('high dose', 'low dose'),
                           milestones = c('dose selection', 'final'),
                           alternative = 'less', 
                           planned_info = 'default')
ct_os <- trial$closedTest(dt_pfs, treatments = c('high dose', 'low dose'),
                          milestones = c('final'),
                          alpha = .02, alpha_spending = 'asOF')

print(ct_pfs)
print(ct_os)
```

The two null hypotheses of `PFS` and `OS` are both accepted for the high
dose because it is dropped at dose selection. Thus, `milestone_at_reject` is
`NA`, and the `reject_time` is infinite. In contrast, `PFS` and `OS` are
significant in low does. Note that even if `PFS` are tested at interim
and final, one cannot claim its significance until the final analysis
(`milestone_at_reject`). The `reject_time` (`r ct_os$reject_time[2]`) can be
saved to the output using the member function `Trial$save`.

`TrialSimulator` abstracts the data generation and management to allow
user focus on implementing the statistical analysis. It simulates a
trial on patient level, and provide flexibility in adaptive design.

## Execute Trial Simulation

We can run a massive number of replicates in simulation to study
operating characteristics of a trial design by specifying `n` in
`Controller$run()`. We can set `plot_event = FALSE` to turn off plotting
to save running time.

```{r ioeinaf, eval=FALSE, cache=TRUE, message=FALSE, warning=FALSE, results='hide'}
## reset a controller if $run has been executed before
controller$reset()
controller$run(n = 1000, plot_event = FALSE, silent = TRUE)
output <- controller$get_output()
```

```{r eaiofj, echo=FALSE}
output <- TrialSimulator:::getAdaptiveDesignOutput()
```

```{r liefa}
output %>% 
  head(5) %>% 
  kable(escape = FALSE) %>% 
  kable_styling(bootstrap_options = "striped", 
                full_width = FALSE,
                position = "left") %>%
  scroll_box(width = "100%")
```

```{r iojalf, class.source = "fold-hide"}
output %>% 
  summarise(
    time_dose_selection = mean(`milestone_time_<dose selection>`), 
    time_interim = mean(`milestone_time_<interim>`), 
    time_final = mean(`milestone_time_<final>`), 
    n_dose_selection = mean(`n_events_<dose selection>_<patient_id>`), 
    n_interim = mean(`n_events_<interim>_<patient_id>`), 
    n_final = mean(`n_events_<final>_<patient_id>`), 
    low = mean(kept_arm == 'low') * 100, 
    high = mean(kept_arm == 'high') * 100, 
    both = mean(kept_arm == 'both') * 100
  ) %>% 
  kable(col.names = NULL, digits = 1, align = 'r', 
        caption = 'Number of Randomized Patients at Stages') %>% 
  add_header_above(c(rep(c('Dose Selection', 'Interim', 'Final'), 2), 
                     'Low Dose', 'High Dose', 'Both'), align = 'r') %>% 
  add_header_above(c('Time' = 3, 'Number of Patients' = 3, 'Selected Dose (%)' = 3)) %>% 
  kable_styling(full_width = TRUE)
```

```{r adalf, class.source = "fold-hide"}
output %>% 
  summarise(
    n_pfs_dose_selection = mean(`n_events_<dose selection>_<pfs>`), 
    n_pfs_interim = mean(`n_events_<interim>_<pfs>`),
    n_pfs_final = mean(`n_events_<final>_<pfs>`),
    n_os_dose_selection = mean(`n_events_<dose selection>_<os>`), 
    n_os_interim = mean(`n_events_<interim>_<os>`),
    n_os_final = mean(`n_events_<final>_<os>`)
  ) %>% 
  kable(col.names = NULL, digits = 1, align = 'r', 
        caption = 'Number of Events of PFS and OS at Stages') %>% 
  add_header_above(rep(c('Dose Selection', 'Interim', 'Final'), 2), align = 'r') %>% 
  add_header_above(c('PFS' = 3, 'OS' = 3)) %>% 
  kable_styling(full_width = TRUE)
```

```{r dgslja, class.source = "fold-hide"}
output %>% 
  summarise(
    power_pfs_low = mean(pfs_low_dose_decision == 'reject') * 100, 
    power_pfs_high = mean(pfs_high_dose_decision == 'reject') * 100, 
    power_pfs_or = mean(pfs_low_dose_decision == 'reject' | pfs_high_dose_decision == 'reject') * 100, 
    power_pfs_and = mean(pfs_low_dose_decision == 'reject' & pfs_high_dose_decision == 'reject') * 100, 
    power_os_low = mean(os_low_dose_decision == 'reject') * 100, 
    power_os_high = mean(os_high_dose_decision == 'reject') * 100, 
    power_os_or = mean(os_low_dose_decision == 'reject' | os_high_dose_decision == 'reject') * 100,
    power_os_and = mean(os_low_dose_decision == 'reject' & os_high_dose_decision == 'reject') * 100
  ) %>% 
  kable(col.names = NULL, digits = 1, align = 'r', 
        caption = 'Power of Testing PFS and OS') %>% 
  add_header_above(rep(c('Low Dose', 'High Dose', 'Low or High', 'Low and High'), 2), align = 'r') %>% 
  add_header_above(c('PFS (%)' = 4, 'OS (%)' = 4)) %>% 
  kable_styling(full_width = TRUE)
```

```{r eioajf, class.source = "fold-hide"}
output %>% 
  summarise(
    power_pfs_not_os = mean((pfs_low_dose_decision == 'reject' | pfs_high_dose_decision == 'reject') & 
                              os_low_dose_decision == 'accept' & os_high_dose_decision == 'accept') * 100, 
    power_os_not_pfs = mean((os_low_dose_decision == 'reject' | os_high_dose_decision == 'reject') & 
                              pfs_low_dose_decision == 'accept' & pfs_high_dose_decision == 'accept') * 100, 
    power_pfs_and_os = mean((pfs_low_dose_decision == 'reject' | pfs_high_dose_decision == 'reject') & 
                              (os_low_dose_decision == 'reject' | os_high_dose_decision == 'reject')) * 100,
    power_pfs_or_os = mean((pfs_low_dose_decision == 'reject' | pfs_high_dose_decision == 'reject') | 
                              (os_low_dose_decision == 'reject' | os_high_dose_decision == 'reject')) * 100
  ) %>% 
  kable(col.names = NULL, digits = 1, align = 'r', 
        caption = 'Power of Testing PFS and OS (Cont.)') %>% 
  add_header_above(c('Reject PFS and Accept OS', 'Accept PFS and Reject OS', 
                     'Reject PFS and OS', 'Reject PFS or OS'), align = 'r') %>% 
  kable_styling(full_width = TRUE)
```