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 compare the simulated operating characteristics with the closed-form values computed with gsDesign.

The point is not that simulation replaces the closed-form calculation under proportional hazards, where gsDesign and similar packages are accurate, 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 whole vignette, with 10,000 simulated trials for the main design, runs in a few seconds.

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. Group 1 is the tisotumab vedotin group and group 2 the chemotherapy group.

nsim <- 10000
eta  <- -log(1 - 0.05) / 12
df <- simdata_fast(
  nsim     = nsim,
  n        = c(241, 241),
  a.time   = c(0, 23),
  a.prop   = 1,
  e.median = list(12.9, 9.0),
  d.hazard = list(eta, eta),
  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.

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.

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

spend_alpha <- 2 * pnorm(gsd$upper$bound, lower.tail = FALSE)
spend_alpha
#> [1] 0.01929865 0.04424346

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.

oc <- simsummary_fast(res, p.col = "logrank.p", alpha = spend_alpha)
oc
#> 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
#> 
#> Events, Sample Size, Dropouts, Pipeline and Analysis Times: Look by Look
#>  Look Info. Frac. Sample (n) Events (s) Dropouts (d) Pipeline Analysis Time
#>     1        0.75      482.0      252.0         16.7    213.3         24.83
#>     2        1.00      482.0      336.0         22.3    123.7         32.86
#>  Cross. Eff.
#>       0.6881
#>       0.2119
#> 
#> 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

The closed-form reference uses the usual large-sample approximation for the log-rank statistic (Schoenfeld, 1983): with 1:1 allocation and d events, the standardized statistic is approximately normal with mean log(HR) * sqrt(d) / 2 in absolute value and unit variance, and the statistics at the two looks have correlation sqrt(252 / 336). gsDesign::gsProbability() gives the crossing probabilities for this drift. The expected calendar time of each look is the time at which the expected number of deaths, under uniform accrual, exponential survival, and exponential dropout, reaches the event target.

looks <- c(252, 336)
theta <- log(12.9 / 9.0) / 2
gp <- gsDesign::gsProbability(k = 2, theta = theta, n.I = looks,
                              a = rep(-20, 2), b = gsd$upper$bound)
p_cross <- gp$upper$prob[, 1]

# Expected number of deaths by calendar time t: uniform accrual of 482
# subjects over 23 months, exponential survival and dropout in each group.
lam <- log(2) / c(12.9, 9.0)
exp_events <- function(t, n_tot = 482, r_acc = 23) {
  a_end <- min(t, r_acc)
  sum(vapply(lam, function(l) {
    cc <- l + eta
    (n_tot / 2) / r_acc * l / cc *
      (a_end - (exp(-cc * (t - a_end)) - exp(-cc * t)) / cc)
  }, numeric(1)))
}
t_look <- vapply(looks, function(d) {
  stats::uniroot(function(t) exp_events(t) - d, c(1, 200))$root
}, numeric(1))

sim_val <- function(lk, col) as.numeric(oc[oc$look == lk, col])
z1 <- -res$logrank.z[res$look == 1]
cmp <- data.frame(
  quantity = c("Mean log-rank Z at the interim",
               "Probability of crossing at the interim",
               "Overall power",
               "Expected number of events at stop",
               "Expected analysis time at stop (months)"),
  simulation = c(mean(z1),
                 sim_val("1", "prob.stop.efficacy"),
                 sim_val("overall", "cum.reject"),
                 sim_val("overall", "n.event.mean"),
                 sim_val("overall", "cutoff.mean")),
  closed_form = c(theta * sqrt(looks[1]),
                  p_cross[1],
                  sum(p_cross),
                  sum(looks * c(p_cross[1], 1 - p_cross[1])),
                  sum(t_look * c(p_cross[1], 1 - p_cross[1])))
)
knitr::kable(cmp, digits = 3,
             col.names = c("Quantity",
                           paste0("FastSurvival (", nsim, " trials)"),
                           "Closed form"))

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. The Monte Carlo standard error of a probability near 0.7 from 10,000 trials is about 0.005, and that of the mean of the interim statistic is about 0.01. The closed form also relies on the large-sample approximation for the log-rank statistic, whose accuracy the first row shows, so the remaining differences reflect both sources of error. The expected analysis time of the closed form is evaluated at the expected number of events, whereas the simulation averages the realized analysis times.

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 (Lin et al., 2020) can be substituted at the analysis_fast() step. Analytic approximations exist for some non-proportional hazards settings and some of these tests, but simulation evaluates any combination of data-generating model, test, and analysis timing in the same way, and it checks such approximations.

The max-combo p-value requires a multivariate normal integral for every simulated trial and look, which dominates the computing time. When only the decisions at the nominal levels are needed, mc.alpha computes the integral only for the p-values that the Bonferroni bounds do not already place on one side of the level; the decisions are the same as with every p-value integrated. The nominal levels below are those of the log-rank design and are used for illustration only: the correlation between the max-combo statistics at the two looks differs from that of the log-rank statistic, so a formal design would derive its boundaries for the max-combo statistic itself.

nsim_nph <- 2000
df_delay <- simdata_fast(
  nsim     = nsim_nph,
  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(eta, eta),
  seed     = 1
)

set.seed(1)
res_delay <- analysis_fast(
  df_delay, control = 2,
  event.looks = looks,
  stat = c("logrank", "maxcombo"), side = 2,
  mc.alpha = spend_alpha
)

oc_lr <- simsummary_fast(res_delay, p.col = "logrank.p",  alpha = spend_alpha)
oc_mc <- simsummary_fast(res_delay, p.col = "maxcombo.p", alpha = spend_alpha)
data.frame(
  test  = c("Log-rank", "Max-combo"),
  power = round(c(oc_lr[oc_lr$look == "overall", "cum.reject"],
                  oc_mc[oc_mc$look == "overall", "cum.reject"]), 3)
)
#>        test power
#> 1  Log-rank 0.924
#> 2 Max-combo 0.959

Of the 4,000 max-combo p-values (two looks per simulated trial), 10% required the multivariate normal integral.

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.

Schoenfeld, D. A. (1983). Sample-size formula for the proportional-hazards regression model. Biometrics, 39(2), 499-503.

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.