This R Markdown document describes how to simulate design characterics for survival design under complex settings (incl. non-proportional hazards) in rpact.
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:
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
For comparison with the simulation-based analysis, a standard example is first calculated under the following assumptions:
<- getDesignGroupSequential(
design informationRates = c(0.66, 1),
typeOfDesign = "asOF", sided = 1, alpha = 0.025, beta = 0.2
)
median2 = 60
) and a target hazard ratio of 0.74
(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)
)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.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:
<- getSampleSizeSurvival(design,
sampleSizeResult 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:
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
.
<- getPowerSurvival(design,
powerResult 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:
These numbers will now be compared to simulations.
The call getSimulationSurvival()
uses the same arguments
as getSampleSizeSurvival()
with the following changes:
maxNumberOfSubjects = 1200
) is always provided to allow
the simulation.plannedEvents = c(232,351)
).directionUpper = FALSE
.maxNumberOfIterations = 10000
in the example below).maxNumberOfRawDatasetsPerStage = 1
.seed = 123
in the example.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.<- getSimulationSurvival(design,
simulationResult 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