---
title: "Group sequential design with the simulation trio"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Group sequential design with the simulation trio}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  eval     = FALSE
)
```

# Purpose

This vignette demonstrates the three simulation functions working together on
a real phase 3 trial: `simdata_fast()` generates the trial data,
`analysis_fast()` performs the interim and final analyses, and
`simsummary_fast()` aggregates the operating characteristics. We reproduce the
group sequential design of the innovaTV 301 trial and check that the simulated
operating characteristics match the closed-form values from
[gsDesign](https://cran.r-project.org/package=gsDesign) and
[rpact](https://cran.r-project.org/package=rpact).

The point is not that simulation replaces the closed-form calculation under
proportional hazards, where gsDesign and rpact are exact, but that once the
simulation agrees on a case they can handle, the same machinery can be trusted
for cases they cannot, such as non-proportional hazards or data-dependent
analysis timing. The benchmark code is shown but not executed when the vignette
is built, because the simulation uses many replicates; the printed results are
those obtained from an interactive run.

```{r load}
library(FastSurvival)
```

# The innovaTV 301 trial

innovaTV 301 (ENGOT-cx12/GOG-3057) was a phase 3, open-label trial of tisotumab
vedotin versus the investigator's choice of chemotherapy in patients with
recurrent or metastatic cervical cancer (Vergote et al., 2024). The primary end
point was overall survival, with patients randomly assigned in a 1:1 ratio.

The design enrolled approximately 482 patients and was powered at 90% on the
occurrence of 336 total deaths, with one prespecified interim efficacy analysis
at about 75% of information (252 of 336 events). The overall two-sided type I
error was controlled at 5% using the Lan-DeMets spending function with an
O'Brien-Fleming boundary. For planning we take an exponential overall survival
with a median of 12.9 months in the tisotumab vedotin group and 9.0 months in
the chemotherapy group (a hazard ratio of about 0.70), accrual over 23 months,
and a 5% annual dropout rate in each group.

# Simulating the trial

`simdata_fast()` generates the survival, censoring, and entry times for all
replicates in one fused C++ pass. We specify the two group sizes, the accrual
window, the per-group event medians, and a per-group dropout hazard
corresponding to 5% per 12 months.

```{r simulate}
df <- simdata_fast(
  nsim     = 10000,
  n        = c(241, 241),
  a.time   = c(0, 23),
  a.prop   = 1,
  e.median = list(12.9, 9.0),
  d.hazard = list(-log(1 - 0.05) / 12, -log(1 - 0.05) / 12),
  seed     = 1
)
```

# Interim and final analyses

`analysis_fast()` runs the event-driven looks. We analyze at 252 and 336 events,
take the chemotherapy group as the control, and compute both the log-rank and
the Cox statistics with a two-sided test, matching the trial.

```{r analyze}
res <- analysis_fast(
  df, control = 2,
  event.looks = c(252, 336),
  stat = c("logrank", "coxph"), side = 2
)
```

# Spending boundaries

The efficacy boundary is the Lan-DeMets O'Brien-Fleming spending function at the
planned information fractions. We obtain it from gsDesign and convert the upper
Z boundaries to two-sided nominal p-value boundaries, which is the scale
`simsummary_fast()` consumes through its `p.col` argument.

```{r boundaries}
library(gsDesign)

gsd <- gsDesign(
  k      = 2,
  timing = c(252, 336) / 336,
  alpha  = 0.025,
  beta   = 0.1,
  sfu    = sfLDOF,
  test.type = 1
)

spend_alpha <- 2 * pnorm(gsd$upper$bound, lower.tail = FALSE)
spend_alpha
```

# Operating characteristics

`simsummary_fast()` applies the nominal p-value boundaries to the simulated
log-rank p-values and aggregates the crossing probabilities, expected events,
expected sample size, and expected analysis time across the looks.

```{r summary}
simsummary_fast(
  res,
  p.col     = "logrank.p",
  alpha     = spend_alpha,
  direction = "lower"
)
```

The interactive run produces the following output.

```
Group-Sequential Operating Characteristics (simsummary_fast)
  Simulations: 10000
  Boundaries: nominal p-value on 'logrank.p'

Stopping Boundaries: Look by Look
 Look Info. Frac. Events (s) Sample (n) Nominal p Cum. Cross. Eff.
    1        0.75      252.0      482.0    0.0193           0.6881
    2        1.00      336.0      482.0    0.0442           0.9000

Overall
  Rejection rate (efficacy):    0.9000
  Expected events at stop:      278.2
  Expected sample size at stop: 482.0
  Expected analysis time at stop:27.31
```

# Comparison with the closed-form design

For an independent closed-form reference we recompute the same design with rpact
and read off the operating characteristics. The full rpact output is long, so we
extract only the quantities needed for the comparison.

```{r rpact}
library(rpact)

design <- getDesignGroupSequential(
  kMax = 2,
  alpha = 0.05,
  beta  = 0.1,
  sided = 2,
  typeOfDesign = "asOF",
  informationRates = c(252, 336) / 336
)

results <- getPowerSurvival(
  design,
  maxNumberOfEvents   = 336,
  median1             = 12.9,
  median2             = 9.0,
  maxNumberOfSubjects = 482,
  accrualTime         = c(0, 23),
  dropoutRate1        = -log(1 - 0.05),
  dropoutRate2        = -log(1 - 0.05),
  allocationRatioPlanned = 1
)
```

The simulated and closed-form operating characteristics agree closely. The
efficacy boundaries are the spending boundaries fed into the simulation, so they
match by construction; the crossing probabilities, expected events, and expected
timing are estimated independently by simulation and line up with the analytic
values.

| Quantity | FastSurvival (10,000 sims) | gsDesign / rpact |
|----------|----------------------------|------------------|
| Interim efficacy boundary (two-sided nominal p) | 0.0193 | 0.0193 |
| Final efficacy boundary (two-sided nominal p) | 0.0442 | 0.0442 |
| Probability of crossing at the interim | 0.688 | 0.698 |
| Overall power | 0.900 | 0.905 |
| Expected number of events at stop | 278 | 277 |
| Expected analysis time at stop (months) | 27.3 | 27.3 |

Remaining differences are Monte Carlo error and shrink as `nsim` increases.

# Beyond proportional hazards

The value of the simulation trio is that nothing in the workflow assumes
proportional hazards. To study a delayed treatment effect, replace the constant
event hazard with a piecewise specification through `e.hazard` and `e.time` and
keep everything else the same. The log-rank statistic loses power under a
delayed effect, and a weighted or max-combo statistic can be substituted at the
`analysis_fast()` step to recover it. Because gsDesign and rpact cannot evaluate
these cases in closed form, the validated simulation machinery becomes the tool
of choice.

```{r nph}
df_delay <- simdata_fast(
  nsim     = 10000,
  n        = c(241, 241),
  a.time   = c(0, 23),
  a.prop   = 1,
  e.hazard = list(c(0.077, 0.045), c(0.077, 0.077)),
  e.time   = c(0, 3, Inf),
  d.hazard = list(-log(1 - 0.05) / 12, -log(1 - 0.05) / 12),
  seed     = 1
)

res_delay <- analysis_fast(
  df_delay, control = 2,
  event.looks = c(252, 336),
  stat = "maxcombo", side = 2
)
```

# References

Vergote, I., González-Martín, A., Fujiwara, K., et al. (2024). Tisotumab
vedotin as second- or third-line therapy for recurrent cervical cancer.
*New England Journal of Medicine*, 391(1), 44-55.

O'Brien, P. C., & Fleming, T. R. (1979). A multiple testing procedure for
clinical trials. *Biometrics*, 35(3), 549-556.

Lan, K. K. G., & DeMets, D. L. (1983). Discrete sequential boundaries for
clinical trials. *Biometrika*, 70(3), 659-663.

Lin, R. S., Lin, J., Roychoudhury, S., et al. (2020). Alternative analysis
methods for time to event endpoints under nonproportional hazards: a
comparative analysis. *Statistics in Biopharmaceutical Research*, 12(2),
187-198.
