#' ---
#' title: "Planning a survival trial with rpact"
#' author: "Gernot Wassmer and Friedrich Pahlke"
#' date: "`r format(Sys.time(), '%d %B, %Y')`"
#' number: 15
#' header-includes:
#' - \usepackage{fancyhdr}
#' - \pagestyle{fancy}
#' - \setlength{\headheight}{23pt}
#' - \fancyfoot[C]{www.rpact.com}
#' - \fancyfoot[L]{\thepage}
#' output:
#' bookdown::html_document2:
#' highlight: pygments
#' number_sections: yes
#' self_contained: yes
#' fig_caption: yes
#' toc: yes
#' toc_depth: 3
#' toc_float: yes
#' css: style.css
#' includes:
#' before_body: header.html
#' after_body: footer.html
#' ---
#'
## ----setup, include=FALSE------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)
options(scipen=999)
#'
#'
#'
#'
#' # Summary {-}
#'
#' This R markdown file provides an example for planning a trial with a survival
#' endpoint using rpact thereby illustrating the different ways of entering
#' recruitment schemes. It also demonstrates the use of the survival simulation
#' function.
#'
#' # Load rpact package {-}
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
# Load rpact
library(rpact)
packageVersion("rpact") # version should be version 2.0.3 or later
#'
#' # The design
#'
#' A survival trial is planned to be performed with one interim stage and using an
#' O'Brien & Fleming type $\alpha$-spending approach at $\alpha = 0.025$. The
#' interim is planned to be performed after half of the necessary events were
#' observed. It is assumed that the median survival time is 18 months in the
#' treatment group, and 12 months in the control. We assume that the drop-out rate
#' is 5% after 1 year and the drop-out time is exponentially distributed.
#'
#' ## Accrual and follow-up time given {#Case1}
#'
#' The patients should be recruited within 12 months assuming uniform accrual.
#' We assume an additional follow-up time of 12 months, i.e., the study should be
#' conducted within 2 years. We also assume the survival time to be exponentially
#' distributed.
#'
#' In this simplest example, accrual and follow-up time needs to be specified and
#' the necessary number of events and patients (total and per month) in order to
#' reach power 90% with the assumed median survival times will be calculated.
#'
#' The effect size is defined in terms of `lambda1` and `lambda2` (you can
#' also specify `lambda2` and `hazardRatio`) and the function `getLambdaByMedian`
#' is used (for `lambda2` the direct calculation is illustrated):
#'
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
dGS <- getDesignGroupSequential(kMax = 2, typeOfDesign = "asOF", beta = 0.1)
x1 <- getSampleSizeSurvival(dGS, lambda1 = getLambdaByMedian(18), lambda2 = log(2)/12,
dropoutRate1 = 0.05, dropoutRate2 = 0.05, dropoutTime = 12,
accrualTime = 12, followUpTime = 12)
x1
#'
#' We conclude that a maximum number of `r ceiling(x1$maxNumberOfEvents)` events
#' needs to be observed to reach 90% power. Under the assumed accrual and
#' follow-up time (and the median survival times) this number of events are
#' expected to be observed if `r ceiling(x1$maxNumberOfSubjects)` subjects are
#' observed at the end of the trial. This corresponds to about
#' `r ceiling(x1$maxNumberOfSubjects)/12`
#' (= `accrualIntensity`) patient per month to be accrued.
#'
#' The interim analysis is performed after `r round(x1$eventsPerStage[1],1)`
#' events (which in practice needs to be rounded up or down). This number is
#' expected to be observed at analysis time `r round(x1$analysisTime[1],2)`.
#'
#' Note that in this case you can also enter the accrual time as
#' `accrualTime = c(0,12)` and so `accrualTime = 12` is just a shortcut (see
#' below).
#'
#' ## Accrual time and maximum number of patients given {#Case2}
#'
#' Assume that accrual stops after 16 months with 25 patients per month, i.e.,
#' after 400 patients were recruited. We now can estimate the necessary follow-up
#' time such that `r ceiling(x1$maxNumberOfEvents)` events will be observed at the
#' end of the trial.
#'
#' There are two ways of defining the situation of *given maximum number of
#' subjects*: You can either specify `accrualTime` and `maxNumberOfSubjects`
#' directly; or you specify `accrualTime` and `accrualIntensity`.
#'
#' In both cases, at given accrual time and number of patients, the follow-up time
#' is calculated.
#'
#' This is the code and the output for the first case:
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
x3 <- getSampleSizeSurvival(dGS, lambda1 = log(2)/18, lambda2 = log(2)/12,
dropoutRate1 = 0.05, dropoutRate2 = 0.05, dropoutTime = 12,
accrualTime = c(0, 16), accrualIntensity = 25)
ceiling(x3$maxNumberOfSubjects)
x3$followUpTime
x3$analysisTime
#'
#' This is the code and the output for the second case:
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
x3 <- getSampleSizeSurvival(dGS, lambda1 = log(2)/18, lambda2 = log(2)/12,
dropoutRate1 = 0.05, dropoutRate2 = 0.05, dropoutTime = 12,
accrualTime = c(0, 16), maxNumberOfSubjects = 400)
x3$accrualIntensity
x3$followUpTime
x3$analysisTime
#'
#' Note that for these two cases entering the accrual time in the form of `c(0,16)`
#' is mandatory, i.e., you cannot enter just the end of accrual which was the
#' shortcut before.
#'
#' Alternatively, you might enter the `accrualTime` as a list in the form
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
atList <- list(
"0 - < 16" = 25)
#'
#' and then enter `accrualTime = at` in `getSampleSizeSurvival()`:
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
x3 <- getSampleSizeSurvival(dGS, lambda1 = log(2)/18, lambda2 = log(2)/12,
dropoutRate1 = 0.05, dropoutRate2 = 0.05, dropoutTime = 12,
accrualTime = atList)
x3$accrualIntensity
x3$followUpTime
x3$analysisTime
#'
#'
#' ## Follow-up time and absolute intensity given {#Case3}
#'
#' Assume now that 25 patients can be recruited each month and that there is
#' uniform accrual. We now want to estimate the necessary accrual time if
#' the *planned follow-up time is given* and equal to 12 as before.
#'
#' Here the end of accrual and the number of patients is calculated at given
#' follow-up time and absolute accrual intensity:
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
x2 <- getSampleSizeSurvival(dGS, hazardRatio = 2/3, lambda2 = log(2)/12,
dropoutRate1 = 0.05, dropoutRate2 = 0.05, dropoutTime = 12,
accrualTime = 0, accrualIntensity = 25, followUpTime = 12)
ceiling(x2$maxNumberOfSubjects)
x2$accrualTime
x2$analysisTime
#'
#'
#' ## Staggered patient entry
#'
#' This case illustrates how the accrual of subjects can be entered if there is
#' piecewise accrual over predefined intervals where it is assumed that accrual
#' within the intervals is uniform.
#'
#' For example, assume that in the first 3 months 15 patients, in the second 3
#' months 20 patients, and after 6 months 25 patients per month can be accrued.
#'
#' If the follow-up time at given number of patients is to be calculated (Section
#' \@ref(Case2)), it is simply to enter
#' `accrualTime = c(0, 3, 6, 16)` and `accrualIntensity = c(15, 20, 25)` as follows
#' (`maxNumberOfSubjects` will be calculated):
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
x5 <- getSampleSizeSurvival(dGS, lambda1 = log(2)/18, lambda2 = log(2)/12,
dropoutRate1 = 0.05, dropoutRate2 = 0.05, dropoutTime = 12,
accrualTime = c(0, 3, 6, 16), accrualIntensity = c(15, 20, 25))
ceiling(x5$maxNumberOfSubjects)
x5$followUpTime
x5$analysisTime
#'
#' Alternatively, you can enter
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
atList <- list(
"0 - < 3" = 15,
"3 - < 6" = 20,
"6 - <= 16" = 25)
#'
#' and run
#'
## ---- include=TRUE, echo=TRUE, eval = FALSE------------------------------
## x5 <- getSampleSizeSurvival(dGS, lambda1 = log(2)/18, lambda2 = log(2)/12,
## dropoutRate1 = 0.05, dropoutRate2 = 0.05, dropoutTime = 12,
## accrualTime = atList)
##
## ceiling(x5$maxNumberOfSubjects)
## x5$followUpTime
## x5$analysisTime
#'
#' If `maxNumberOfSubjects` is specified, the corresponding list definition is as
#' follows:
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
atList <- list(
"0 - < 3" = 15,
"3 - < 6" = 20,
"6 - Inf" = 25)
#'
#' and
#'
## ---- include=TRUE, echo=TRUE, eval = FALSE------------------------------
## x5 <- getSampleSizeSurvival(dGS, lambda1 = log(2)/18, lambda2 = log(2)/12,
## dropoutRate1 = 0.05, dropoutRate2 = 0.05, dropoutTime = 12,
## accrualTime = atList, maxNumberOfSubjects = 355)
##
## ceiling(x5$accrualTime)
## x5$followUpTime
## x5$analysisTime
#'
#' does the job.
#'
#'
#' If the follow-up time is given and the end of accrual it to be calculated
#' (Section \@ref(Case3)), this can be achieved as follows:
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
x4 <- getSampleSizeSurvival(lambda1 = log(2)/18, lambda2 = log(2)/12,
accrualTime = c(0, 3, 6), accrualIntensity = c(15, 20, 25), followUpTime = 12)
ceiling(x4$maxNumberOfSubjects)
x4$accrualTime
x4$analysisTime
#'
#'
#' # Verify results by simulation
#'
#' Assume that the study is planned with 257 events and 400
#' patients under the assumptions that accrual stops after 16 months with 25
#' patients per month and the staggered patient entry from above. We now want
#' to verify by simulation the correctness of the results obtained by the
#' analytical formulae.
#'
#' For this, we first use the function `getPowerSurvival()` in order to obtain
#' the test characteristics if the study is performed with 257 events. The
#' analysis time is obtained by the analytical formulae and we verify that the
#' power is indeed exceeding 90% (note that we have to specify `directionupper = FALSE`
#' in order to indicate that p-values are getting small for *negative* values of
#' the test statistics, i.e., for observed hazard ratios smaller than 1):
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
y3 <- getPowerSurvival(dGS, lambda1 = log(2)/18, lambda2 = log(2)/12,
dropoutRate1 = 0.05, dropoutRate2 = 0.05, dropoutTime = 12,
accrualTime = c(0, 3, 6, 16), accrualIntensity = c(15, 20, 25),
maxNumberOfEvents = 257, directionUpper = FALSE)
y3$analysisTime
y3$overallReject
#'
#' Practically the same result is obtained with the simulation function
#' `getSimulationSurvival()` (note that `plannedEvents = c(129,257)`
#' needs to be specified in teh simulation tool instead of `maxNumberOfEvents =
#' 257` in the power calculation too:
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
maxiter <- 10000
z3 <- getSimulationSurvival(dGS, lambda1 = log(2)/18, lambda2 = log(2)/12,
dropoutRate1 = 0.05, dropoutRate2 = 0.05, dropoutTime = 12, maxNumberOfIterations = maxiter,
accrualTime = c(0, 3, 6, 16), accrualIntensity = c(15, 20, 25),
plannedEvents = c(129, 257), directionUpper = FALSE, seed = 12345)
z3$analysisTime
z3$overallReject
#'
#' # Assess sample size reassessment in adaptive survival design
#'
#' Assume now that a sample size increase up to a ten-fold of the originally
#' planned number of events is foreseen. Conditional power 90% *based on the
#' observed hazard ratios* is used to increase the number of events. We can
#' now assess by simulation the magnitute of power increase when using the
#' appropriate method.
#'
#' Note that we have to make sure that enough subjects are used in the simulation.
#' For this, we set `maxNumberOfSubjects = 3000` and no drop-outs.
#'
#' Since we perform a data driven sample size recalculation (SSR), we define an
#' inverse normal design with the same parameters as the original group sequential
#' design:
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
dIN <- getDesignInverseNormal(kMax = 2, typeOfDesign = "asOF", beta = 0.1)
#'
#' The following code simulates the SSR with conditional power 90% under the
#' observed hazard ratio. `minNumberOfEventsPerStage = c(NA,128)`,
#' and `maxNumberOfEventsPerStage = 10*c(NA,128)` defines:
#'
#' - no sample size decrease
#'
#' - up to a 10 fold increase of number of events for the final stage.
#'
#' This is compared with the situation without SSR:
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
# Without SSR
z0 <- getPowerSurvival(dIN, lambda1 = log(2)/18, lambda2 = log(2)/12,
maxNumberOfIterations = maxiter,
accrualTime = c(0,16), maxNumberOfSubjects = 3000, maxNumberOfEvents = 257,
directionUpper = FALSE)
z0$analysisTime
z0$overallReject
#'
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
# With SSR
z4 <- getSimulationSurvival(dIN, lambda1 = log(2)/18, lambda2 = log(2)/12,
maxNumberOfIterations = maxiter,
accrualTime = c(0,16), maxNumberOfSubjects = 3000, plannedEvents = c(129, 257),
directionUpper = FALSE, conditionalPower = 0.9,
minNumberOfEventsPerStage = c(NA,128),
maxNumberOfEventsPerStage = 10 * c(NA,128), seed = 23456)
z4$analysisTime
z4$overallReject
#'
#' Finally, we simulate the Type I error rate when using
#'
#' - the group sequential method
#'
#' - the inverse normal method
#'
#' for analysing the data. The difference here is that the group sequential
#' design uses the logrank tests for the test decision whereas with the
#' combination test the *increments* of the sequence of logrank statistics are
#' combined with the use of the inverse normal combination function. This assures
#' Type I error control also the data-driven SSR.
#'
#' The following simulation compares the Type I error rate of the inverse normal
#' method with the Type I error rate of the (inappropriate) use of the group
#' sequential method:
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
maxiter <- 10000
dGS <- getDesignGroupSequential(kMax = 2, typeOfDesign = "asOF")
dIN <- getDesignInverseNormal(kMax = 2, typeOfDesign = "asOF")
IN <- getSimulationSurvival(dIN, hazardRatio = 1,
maxNumberOfIterations = maxiter,
accrualTime = c(0, 16), maxNumberOfSubjects = 3000, plannedEvents = c(129, 257),
directionUpper = FALSE, conditionalPower = 0.9,
minNumberOfEventsPerStage = c(NA, 128),
maxNumberOfEventsPerStage = 10 * c(NA,128), seed = 34567)
GS <- getSimulationSurvival(dGS, hazardRatio = 1,
maxNumberOfIterations = maxiter,
accrualTime = c(0,16), maxNumberOfSubjects = 2000, plannedEvents = c(129, 257),
directionUpper = FALSE, conditionalPower = 0.9, minNumberOfEventsPerStage = c(NA,128),
maxNumberOfEventsPerStage = 10 * c(NA,128), seed = 45678)
IN$overallReject
GS$overallReject
#'
#' We see that indeed the Type I error rate is controlled for the inverse normal
#' but nor for the group sequential method (see also [this
#' vignette](https://vignettes.rpact.org/html/rpact_ggplot_examples.html)).
#'
#'
#' # References
#'
#' 1. Gernot Wassmer and Werner Brannath, *Group Sequential and Confirmatory Adaptive
#' Designs in Clinical Trials*, Springer 2016, ISBN 978-3319325606
#'
#' ***
#'
#' System: rpact `r packageVersion("rpact")`, `r R.version.string`, platform: `r R.version$platform`
#'
## ---- include=TRUE, echo=FALSE, results='asis'---------------------------
print(citation("rpact"), bibtex = FALSE)
#'