Summary

This R Markdown document describes how to simulate design characterics for survival design under complex settings (incl. non-proportional hazards) in rpact.

1 Introduction

rpact provides the function getSimulationSurvival() for simulation of group sequential trials with a time-to-event endpoint. For a given scenario, getSimulationSurvival() simulates many hypothetical group sequential trials and calculates the test results. Based on this Monte Carlo simulation, estimates of key quantities such as overall study power, stopping probabilities at each interim analysis, timing of analyses etc. can be obtained.

getSimulationSurvival() complements the analytical calculations from the function getPowerSurvival() (and getSampleSizeSurvival()) in multiple ways:

Simulations can be used to assess the accuracy of the analytical formulas.
Simulations allow to answer questions such as the following:
- How variable is the timing of interim analysis (even if all assumptions are correct)?
- How could a dataset of a trial that is stopped early for efficacy at an interim analysis look like?
Simulation is also possible for scenarios that are analytically intractable such as scenarios with delayed treatment effects.

The syntax of function getSimulationSurvival() is very similar to the function getPowerSurvival() and getSampleSizeSurvival(). Hence, this document only provides some examples and expects that the reader is familiar with the R Markdown document Designing group sequential trials with two groups and a survival endpoint with rpact which describes standard designs of a trial with a survival endpoint.

getSimulationSurvival() also supports the usage of adaptive sample size recalculation but this is not covered here. For more details, please also consult the help ?getSimulationSurvival.

First, load the rpact package

library(rpact)
packageVersion("rpact") # version should be version 2.0.5 or later

## [1] '3.3.2'

For this R Markdown document, some additional R packages are loaded which will be used to examine simulated datasets.

# Load ggplot2 for plotting
library(ggplot2) # general plotting
library(survminer) # plotting of KM curves with ggplot2
library(survival) # survival analysis routines

2 Standard analytical calculation

For comparison with the simulation-based analysis, a standard example is first calculated under the following assumptions:

Group-sequential design with one interim analysis after 66% of information using an O’Brien & Fleming type \(\alpha\)-spending function, one-sided Type I error 2.5%, power 80%:

design <- getDesignGroupSequential(
    informationRates = c(0.66, 1),
    typeOfDesign = "asOF", sided = 1, alpha = 0.025, beta = 0.2
)

Exponential PFS with a median PFS of 60 months in control (median2 = 60) and a target hazard ratio of 0.74 (hazardRatio = 0.74).
Annual drop-out of 2.5% in both arms (dropoutRate1 = 0.025, dropoutRate2 = 0.025, dropoutTime = 12).
Recruitment is 42 patients/month from month 6 onwards after linear ramp up. (accrualTime = c(0,1,2,3,4,5,6), accrualIntensity = c(6,12,18,24,30,36,42))
Randomization ratio 1:1 (allocation1 = 1 and allocation2 = 1); this is how subjects are randomized in treatment groups 1 and 2 in a subsequent way). 1:1 allocation is the default and is thus not explicitly set in the function call below.
A fixed total sample size of 1200 (maxNumberOfSubjects = 1200).

As described in the R Markdown document which describes standard designs of a trial with a survival endpoint, sample size calculations for this design can be performed as per the code below:

sampleSizeResult <- getSampleSizeSurvival(design,
    median2 = 60, hazardRatio = 0.74,
    dropoutRate1 = 0.025, dropoutRate2 = 0.025, dropoutTime = 12,
    accrualTime = c(0, 1, 2, 3, 4, 5, 6),
    accrualIntensity = c(6, 12, 18, 24, 30, 36, 42),
    maxNumberOfSubjects = 1200
)

# Summary of results
kable(summary(sampleSizeResult))

Sample size calculation for a survival endpoint

Sequential analysis with a maximum of 2 looks (group sequential design), overall significance level 2.5% (one-sided). The sample size was calculated for a two-sample logrank test, H0: hazard ratio = 1, H1: hazard ratio = 0.74, control median(2) = 60, maximum number of subjects = 1200, accrual time = c(1, 2, 3, 4, 5, 6, 31.571), accrual intensity = c(6, 12, 18, 24, 30, 36, 42), dropout rate(1) = 0.025, dropout rate(2) = 0.025, dropout time = 12, power 80%.

Stage	1	2
Information rate	66%	100%
Efficacy boundary (z-value scale)	2.524	1.992
Overall power	0.4074	0.8000
Expected number of subjects	1200.0
Number of subjects	1200.0	1200.0
Cumulative number of events	231.3	350.5
Expected study duration	47.9
Cumulative alpha spent	0.0058	0.0250
One-sided local significance level	0.0058	0.0232
Efficacy boundary (t)	0.718	0.808
Exit probability for efficacy (under H0)	0.0058
Exit probability for efficacy (under H1)	0.4074

Legend:

(t): treatment effect scale

By design, the power of the trial is 80%. The interim analysis is after 232 events which is expected to occur after 39.5 months, and the final analysis is after 351 events which is expected to occur after 53.7 months.

Use getPowerSurvival() to calculate the power achieved for the ceiled number of events. Note that the direction of the alternative needs to be specified. Here, the alternative is towards hazard ratios < 1 which is specified as directionUpper = FALSE.

powerResult <- getPowerSurvival(design,
    maxNumberOfEvents = ceiling(sampleSizeResult$maxNumberOfEvents),
    directionUpper = FALSE,
    median2 = 60, hazardRatio = 0.74,
    dropoutRate1 = 0.025, dropoutRate2 = 0.025, dropoutTime = 12,
    accrualTime = c(0, 1, 2, 3, 4, 5, 6),
    accrualIntensity = c(6, 12, 18, 24, 30, 36, 42),
    maxNumberOfSubjects = 1200
)

# Summary of results
kable(summary(powerResult))

Power calculation for a survival endpoint

Sequential analysis with a maximum of 2 looks (group sequential design), overall significance level 2.5% (one-sided). The results were calculated for a two-sample logrank test, H0: hazard ratio = 1, power directed towards smaller values, H1: hazard ratio = 0.74, control median(2) = 60, maximum number of subjects = 1200, maximum number of events = 351, accrual time = c(1, 2, 3, 4, 5, 6, 31.571), accrual intensity = c(6, 12, 18, 24, 30, 36, 42), dropout rate(1) = 0.025, dropout rate(2) = 0.025, dropout time = 12.

Stage	1	2
Information rate	66%	100%
Efficacy boundary (z-value scale)	2.524	1.992
Overall power	0.4080	0.8005
Expected number of subjects	1200.0
Number of subjects	1200.0	1200.0
Expected number of events	302.3
Cumulative number of events	231.7	351.0
Expected study duration	47.9
Cumulative alpha spent	0.0058	0.0250
One-sided local significance level	0.0058	0.0232
Efficacy boundary (t)	0.718	0.808
Exit probability for efficacy (under H0)	0.0058
Exit probability for efficacy (under H1)	0.4080

Legend:

(t): treatment effect scale

These numbers will now be compared to simulations.

3 Simulation under proportional hazards

The call getSimulationSurvival() uses the same arguments as getSampleSizeSurvival() with the following changes:

The maximum number of patients (maxNumberOfSubjects = 1200) is always provided to allow the simulation.
The number of events at each analysis is specified as per the analytical calculation above (plannedEvents = c(232,351)).
The direction of the alternative is specified as directionUpper = FALSE.
The number of simulated trials is specified (maxNumberOfIterations = 10000 in the example below).
By default, raw datasets from simulation runs are not extracted. However, in this example, it is specifies that one raw dataset that led to stopping after each stage, respectively, will be stored: maxNumberOfRawDatasetsPerStage = 1.
For reproducibility, it is useful to set the random seed which is set to seed = 123 in the example.
For simulation of trials with unequal randomization, integer arguments allocation1 and allocation2 must be provided to the function getSimulationSurvival (instead of argument allocationRatioPlanned in the function getSampleSizeSurvival). allocation1 and allocation2 specify the number of consecutively enrolled subjects in the intervention and control groups, respectively, before another subject from the opposite group is recruited. For example, allocation1 = 2, allocation2 = 1 refers to 2 (intervention):1 (control) randomization and allocation1 = 1, allocation2 = 2 to 1 (intervention):2 (control) randomization.

simulationResult <- getSimulationSurvival(design,
    median2 = 60, hazardRatio = 0.74,
    dropoutRate1 = 0.025, dropoutRate2 = 0.025,
    dropoutTime = 12,
    accrualTime = c(0, 1, 2, 3, 4, 5, 6),
    accrualIntensity = c(6, 12, 18, 24, 30, 36, 42),
    plannedEvents = c(232, 351),
    directionUpper = FALSE,
    maxNumberOfSubjects = 1200,
    maxNumberOfIterations = 10000,
    maxNumberOfRawDatasetsPerStage = 1,
    seed = 234
)

# Summary of simulation results
kable(summary(simulationResult))

Simulation of a survival endpoint

Sequential analysis with a maximum of 2 looks (group sequential design), overall significance level 2.5% (one-sided). The results were simulated for a two-sample logrank test, H0: hazard ratio = 1, power directed towards smaller values, H1: hazard ratio = 0.74, control median(2) = 60, planned cumulative events = c(232, 351), maximum number of subjects = 1200, accrual time = c(1, 2, 3, 4, 5, 6, 31.571), accrual intensity = c(6, 12, 18, 24, 30, 36, 42), dropout rate(1) = 0.025, dropout rate(2) = 0.025, dropout time = 12, simulation runs = 10000, seed = 234.

Stage	1	2
Fixed weight	0.66	1
Efficacy boundary (z-value scale)	2.524	1.992
Stage Levels	0.0058	0.0232
Overall power	0.3974	0.7959
Expected number of subjects	1200.0
Number of subjects	1200.0	1200.0
Expected number of events	303.7
Cumulative number of events	232.0	119.0
Expected study duration	48.1
Exit probability for efficacy	0.3974	0.3985

According to the output, the simulated overall power is 79.6% and the probability to cross the efficacy boundary at the interim analysis is 39.7%. These are both within 1% of the analytical power.

The mean simulated analysis times are after 39.61 months for the interim analysis and after 53.66 for the final analysis. Both timings differ by <0.1 months from the analytical calculation (Difference analysis times = 0.06, 0.01).

You can show median [range] and mean+/-sd of the trial results across the simulate

Simulation-Based Design of Group Sequential Trials with a Survival Endpoint with rpact

Marcel Wolbers, Gernot Wassmer, and Friedrich Pahlke

Last change: 15 November, 2022

Summary

1 Introduction

2 Standard analytical calculation

3 Simulation under proportional hazards