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 getSampleSizeSurvival()
in multiple ways:
The syntax of function getSimulationSurvival()
is very similar to the function 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
## [1] '3.2.1'
For this R Markdown document, some additional R packages are loaded which will be used to examine simulated datasets.
For comparison with the simulation-based analysis, a standard example is first calculated under the following assumptions:
design <- getDesignGroupSequential(informationRates = c(0.66,1),
typeOfDesign = "asOF", sided = 1, alpha = 0.025, beta = 0.2)
lambda2 = log(2)/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:
sampleSizeResult <- getSampleSizeSurvival(design,
lambda2 = log(2)/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)
# Standard rpact display of results
sampleSizeResult
## Design plan parameters and output for survival data:
##
## Design parameters:
## Information rates : 0.660, 1.000
## Critical values : 2.524, 1.992
## Futility bounds (non-binding) : -Inf
## Cumulative alpha spending : 0.005798, 0.025000
## Local one-sided significance levels : 0.005798, 0.023210
## Significance level : 0.0250
## Type II error rate : 0.2000
## Test : one-sided
##
## User defined parameters:
## lambda(2) : 0.0116
## Hazard ratio : 0.740
## Maximum number of subjects : 1200
## Accrual time : 1.00, 2.00, 3.00, 4.00, 5.00, 6.00, 31.57
## Accrual intensity : 6, 12, 18, 24, 30, 36, 42
## Drop-out rate (1) : 0.025
## Drop-out rate (2) : 0.025
##
## Default parameters:
## Theta H0 : 1
## Type of computation : Schoenfeld
## Planned allocation ratio : 1
## kappa : 1
## Piecewise survival times : 0.00
## Drop-out time : 12.00
##
## Sample size and output:
## Direction upper : FALSE
## median(1) : 81.1
## median(2) : 60.0
## lambda(1) : 0.00855
## Maximum number of subjects (1) : 600.0
## Maximum number of subjects (2) : 600.0
## Maximum number of events : 350.5
## Total accrual time : 31.57
## Follow up time : 22.08
## Reject per stage [1] : 0.4074
## Reject per stage [2] : 0.3926
## Early stop : 0.4074
## Analysis times [1] : 39.54
## Analysis times [2] : 53.65
## Expected study duration : 47.91
## Maximal study duration : 53.65
## Cumulative number of events [1] : 231.3
## Cumulative number of events [2] : 350.5
## Expected number of events under H0 : 349.8
## Expected number of events under H0/H1 : 340.5
## Expected number of events under H1 : 302.0
## Number of subjects [1] : 1200.0
## Number of subjects [2] : 1200.0
## Expected number of subjects under H1 : 1200.0
## Critical values (treatment effect scale) [1] : 0.718
## Critical values (treatment effect scale) [2] : 0.808
##
## Legend:
## (i): values of treatment arm i
## [k]: values at stage k
## 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 lambda(2) = 0.012,
## 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.54 months, and the final analysis is after 351 events which is expected to occur after 53.65 months. 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 = 2
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, allocation1 = 1
refers to 2 (intervention):1 (control) randomization and allocation1 = 1, allocation1 = 2
to 1 (intervention):2 (control) randomization.simulationResult <- getSimulationSurvival(design,
lambda2 = log(2)/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 = 12345)
# print the simulation result
simulationResult
## Simulation of survival data (group sequential design):
##
## Design parameters:
## Information rates : 0.660, 1.000
## Critical values : 2.524, 1.992
## Futility bounds (non-binding) : -Inf
## Cumulative alpha spending : 0.005798, 0.025000
## Local one-sided significance levels : 0.005798, 0.023210
## Significance level : 0.0250
## Test : one-sided
##
## User defined parameters:
## Maximum number of iterations : 10000
## Seed : 12345
## Direction upper : FALSE
## Planned cumulative events : 232, 351
## Maximum number of subjects : 1200
## Accrual time : 1.00, 2.00, 3.00, 4.00, 5.00, 6.00, 31.57
## Accrual intensity : 6, 12, 18, 24, 30, 36, 42
## Drop-out rate (1) : 0.025
## Drop-out rate (2) : 0.025
## Piecewise survival times : 0
## lambda(2) : 0.0116
## Hazard ratio : 0.740
##
## Default parameters:
## Planned allocation ratio : 1
## Conditional power : NA
## Drop-out time : 12.00
## Theta H0 : 1
## Allocation 1 : 1
## Allocation 2 : 1
## kappa : 1
##
## Results:
## median(1) : 81.1
## median(2) : 60.0
## lambda(1) : 0.00855
## Analysis times [1] : 39.61
## Analysis times [2] : 53.69
## Expected study duration : 48.04
## Events not achieved [1] : 0.000
## Events not achieved [2] : 0.000
## Number of subjects [1] : 1200.0
## Number of subjects [2] : 1200.0
## Iterations [1] : 10000
## Iterations [2] : 5983
## Overall reject : 0.7949
## Reject per stage [1] : 0.4017
## Reject per stage [2] : 0.3932
## Futility stop per stage : 0.0000
## Early stop : 0.4017
## Expected number of subjects : 1200.0
## Expected number of events : 303.2
## Cumulative number of events [1] : 232.0
## Cumulative number of events [2] : 351.0
## Conditional power (achieved) [1] : NA
## Conditional power (achieved) [2] : 0.5402
##
## Legend:
## (i): values of treatment arm i
## [k]: values at stage k
## 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 lambda(2) = 0.012,
## 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 = 12345.
##
## Stage 1 2
## Fixed weight 0.66 1
## Efficacy boundary (z-value scale) 2.524 1.992
## Overall power 0.4017 0.7949
## Expected number of subjects 1200.0
## Number of subjects 1200.0 1200.0
## Expected number of events 303.2
## Cumulative number of events 232.0 351.0
## Expected study duration 48.0
## Exit probability for efficacy 0.4017 0.3932
According to the output, the simulated overall power is 79.5% and the probability to cross the efficacy boundary at the interim analysis is 40.2%. 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.69 for the final analysis. Both timings differ by <0.1 months from the analytical calculation (Difference analysis times = 0.07, 0.04).
The summary of the simulations above also provide median [range] and mean+/-sd of the trial results across the simulated datasets. Summary results for each trial can be obtained from the simulation object using the function getData()
. Similarly, raw data from individual trials that were stopped at each stage can be obtained using function getRawData()
(if maxNumberOfRawDatasetsPerStage
was set > 0). The format of these datasets is described in the help ?getSimulationSurvival
and illustrated below.
# get aggregate datasets from all simulation runs
aggregateSimulationData <- getData(simulationResult)
# show the first 6 records of the dataset
head(aggregateSimulationData)
## iterationNumber stageNumber pi1 pi2 hazardRatio analysisTime numberOfSubjects
## 1 1 1 NA NA 0.74 40.64430 1200
## 2 2 1 NA NA 0.74 39.98191 1200
## 3 2 2 NA NA 0.74 52.27530 1200
## 4 3 1 NA NA 0.74 41.02950 1200
## 5 3 2 NA NA 0.74 55.05094 1200
## 6 4 1 NA NA 0.74 40.19759 1200
## overallEvents1 overallEvents2 eventsPerStage rejectPerStage eventsNotAchieved
## 1 99 133 232 1 0
## 2 104 128 232 0 0
## 3 160 191 119 1 0
## 4 100 132 232 0 0
## 5 154 197 119 1 0
## 6 107 125 232 0 0
## futilityPerStage testStatistic logRankStatistic conditionalPowerAchieved
## 1 0 2.580807 2.580807 NA
## 2 0 1.873088 1.873088 NA
## 3 0 2.148403 2.148403 0.7039150
## 4 0 2.253289 2.253289 NA
## 5 0 2.644996 2.644996 0.9095044
## 6 0 1.416149 1.416149 NA
## pValuesSeparate trialStop hazardRatioEstimateLR
## 1 NA TRUE 0.7125704
## 2 NA FALSE 0.7819625
## 3 NA TRUE 0.6744299
## 4 NA FALSE 0.7438831
## 5 NA TRUE 0.6157386
## 6 NA FALSE 0.8303156
# The aggregated dataset contains one record for each of the 10'000 simulated
# interim datasets (stageNumber = 1), and one additional record for the final
# analysis (stageNumber = 2) for all simulated trials which were not stopped at
# the interim analysis
table(aggregateSimulationData$stageNumber)
##
## 1 2
## 10000 5983
# One possible analysis which uses the aggregate dataset:
# display the distribution of the timing of the analyses (by stage)
# across simulation runs
p <- ggplot(aggregateSimulationData, aes(factor(stageNumber),analysisTime))
p + geom_boxplot() + labs(x = "Analysis stage",
y = "Timing of analysis (months \nfrom first patient randomized)")
# get all stored raw datasets
rawSimulationDataData <- getRawData(simulationResult)
# chose one dataset which corresponds to a simulated (hypothetical) trial
# which was stopped at the efficacy interim analysis
# (only 1 such dataset exists here, because maxNumberOfRawDatasetsPerStage
# was set to 1 above)
stage1rawData <- subset(rawSimulationDataData,stopStage = 1)
# Perform a Cox regression analysis of the simulated trial
# (treatmentGroup == 1 is the intervention arm)
# Note that the $p$-value from the Cox regression is two-sided whereas
# the group-sequential is one-sided.
summary(coxph(Surv(timeUnderObservation,event)~I(treatmentGroup == 1),data = stage1rawData))
## Call:
## coxph(formula = Surv(timeUnderObservation, event) ~ I(treatmentGroup ==
## 1), data = stage1rawData)
##
## n= 2400, number of events= 583
##
## coef exp(coef) se(coef) z Pr(>|z|)
## I(treatmentGroup == 1)TRUE -0.27458 0.75989 0.08336 -3.294 0.000988 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## I(treatmentGroup == 1)TRUE 0.7599 1.316 0.6453 0.8948
##
## Concordance= 0.542 (se = 0.011 )
## Likelihood ratio test= 10.93 on 1 df, p=0.0009
## Wald test = 10.85 on 1 df, p=0.001
## Score (logrank) test = 10.92 on 1 df, p=0.001
# Show Kaplan-Meier curves of the simulated trial
ggsurvplot(survfit(Surv(timeUnderObservation, event) ~ treatmentGroup,
data = stage1rawData), risk.table = TRUE, break.time.by = 12,
censor.size = 1, risk.table.fontsize = 4, risk.table.pos = "in",
title = paste0("Kaplan-Meier curve of one simulated trial\n which ",
"was stopped at the interim analysis"),
xlab = "Time (months)", data = stage1rawData)