#' ---
#' title: "Designing group-sequential trials with two groups and a survival endpoint with rpact"
#' author: "Marcel Wolbers, Gernot Wassmer, and Friedrich Pahlke"
#' date: "Last change: `r format(Sys.time(), '%d %B, %Y')`"
#' number: 4
#' 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: footer2.html
#' ---
#'
#'
#'
#' # Summary {-}
#'
#' This [R Markdown](https://rmarkdown.rstudio.com) document provides examples for
#' designing trials with survival endpoints with [rpact](https://cran.r-project.org/package=rpact).
#'
#'
#' # Introduction
#'
#' The examples in this document are not intended to replace the official rpact
#' documentation and help pages but rather to supplement them. They also only cover
#' a selection of all rpact features.
#'
#' General convention: In rpact, arguments containing the **index "2"** always refer
#' to the **control group**, **"1"** refer to the **intervention group**, and
#' **treatment effects compare treatment versus control**.
#'
#' **First, load the rpact package**
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
library(rpact)
packageVersion("rpact") # version should be version 2.0.5 or later
#'
#'
#' # Overview of relevant rpact functions
#'
#' Sample size calculation for a group-sequential trial with a survival endpoints
#' is performed in **two steps**:
#'
#' 1. **Define the (abstract) group-sequential design** using the function
#' `getDesignGroupSequential()`. For details regarding this step, see the R
#' markdown file "Defining group-sequential boundaries with rpact". This step
#' can be omitted if the trial has no interim analyses.
#'
#' 2. **Calculate sample size** for the survival endpoint by feeding the abstract
#' design into function `getSampleSizeSurvival()`. This step is described in
#' detail below.
#'
#' Other relevant rpact functions for survival are:
#'
#' - `getPowerSurvival()`: This function is the analogue to `getSampleSizeSurvival()`
#' for the calculation of power rather than sample size.
#' - `getEventProbabilities()`: Calculates the probability of an event depending on
#' follow-up time. This is useful for aligning interim analyses for different
#' time-to-event endpoints.
#' - `getSimulationSurvival()`: This function simulates group-sequential trials.
#' For example, it allows to assess the power of trials with delayed treatment
#' effects or to assess the data-dependent variability of the timing of interim
#' analyses even if the protocol assumptions are perfectly fulfilled. It also
#' allows to simulate hypothetical datasets for trials stopped early.
#'
#' This document describes all functions mentioned above except for trial
#' simulation (`getSimulationSurvival()`) which is described in the document
#' [Simulation-based design of group-sequential trials with a survival endpoint
#' with rpact](https://vignettes.rpact.org/html/rpact_survival_examples.html).
#'
#' However, before describing the functions themselves, the document describes how
#' survival functions, drop-out, and accrual can be specified in rpact which is
#' required for all of these functions.
#'
#' # Specifying survival distributions in rpact
#'
#' rpact allows to specify survival distributions with exponential, piecewise
#' exponential, and Weibull distributions. Exponential and piecewise exponential
#' distributions are described below.
#'
#' Weibull distributions are specified in the same way as exponential distributions
#' except that an additional scale parameter `kappa` needs to be provided which is
#' 1 for the exponential distribution. Note that the parameters `shape` and
#' `scale` in the standard R functions for the Weibull distribution in the
#' stats-library (such as `dweibull`) are equivalent to `kappa` and 1/`lambda`,
#' respectively, in rpact.
#'
#' ## Exponential survival distributions
#'
#' ### Event probability at a specific time point known
#'
#' - The time point is given by argument `eventTime`.
#' - The probability of an event (i.e., 1 minus survival function) in the control
#' group is given by argument `pi2`.
#' - The probability of an event in the intervention arm can either be provided
#' explicitly through argument `pi1` or, alternatively, implicitly by specifying
#' the target hazard ratio `hazardRatio`.
#'
#' **Example:** If the intervention is expected to improve survival at 24 months
#' from 70% to 80%, this would be expressed through arguments
#' `eventTime=24, pi2=0.3, pi1=0.2`.
#'
#' ### Exponential parameter $\lambda$ known
#'
#' - The constant hazard function in the control arm can be provided as argument
#' `lambda2`.
#' - The hazard in the intervention arm can either be provided explicitly through
#' the argument `lambda1` or, alternatively, implicitly by specifying the target
#' hazard ratio `hazardRatio`.
#'
#' ### Median survival known
#'
#' Medians cannot be specified directly. However, one can exploit that the hazard
#' rate $\lambda$ of the exponential distribution is equal to
#' $\lambda=\log(2)/\text{median}$.
#'
#' **Example:** If the intervention is expected to improve the median survival from
#' 60 to 75 months, this would be expressed through the arguments
#' `lambda2=log(2)/60, lambda1=log(2)/75`. Alternatively, it could be specified via
#' `lambda2=log(2)/60, hazardRatio=0.8` (as the hazard ratio is 60/75=0.8).
#'
#' ## Weibull
#'
#' Weibull distributions are specified in the same way as exponential distributions
#' except that an additional scale parameter `kappa` needs to be provided which is
#' 1 for the exponential distribution.
#'
#'
#' ## Piecewise exponential survival
#'
#' ### Piecewise constant hazard rate $\lambda$ in each interval known
#'
#' **Example:** If the survival function in the control arm is assumed to be 0.03
#' (events/time-unit of follow-up) for the first 6 months, 0.06 for months 6-12,
#' and 0.02 from month 12 onwards, this can be specified using the argument
#' `piecewiseSurvivalTime` as follows:
#'
#' ```
#' piecewiseSurvivalTime=list("0 - <6" =0.03,
#' "6 - <12"=0.06,
#' ">=12" =0.02)
#' ```
#' Alternatively, the same distribution can be specified by giving the start times
#' of each interval as argument `piecewiseSurvivalTime` and the actual hazard rate
#' in that interval as argument `lambda2`. I.e., the relevant arguments for this
#' example would be:
#'
#' ```
#' piecewiseSurvivalTime=c(0,6,12), lambda2=c(0.03,0.06,0.02)
#' ```
#'
#' For the intervention arm, one could either explicitly specify the hazard rate in
#' the same time intervals through the argument `lambda1` or, more conveniently,
#' specify the survival function in the intervention arm indirectly via the target
#' hazard ratio (argument `hazardRatio`).
#'
#' Note: The sample size calculation functions discussed in this document **assume
#' that the proportional hazards assumption holds**, i.e., that `lambda1` and
#' `lambda2` are proportional if both are provided. Otherwise, the function call to
#' `getSampleSizeSurvival()` gives an error. For situations with non-proportional
#' hazards, please see the separate document which discusses the simulation tool
#' using the function `getSimulationSurvival()`.
#'
#' ### Survival function at different time points known
#'
#' Piecewise exponential distribution are useful to **approximate survival functions**.
#' In principle, it is possible to approximate any distribution function well with
#' a piecewise exponential distribution if suitably narrow intervals are chosen.
#'
#' Assume that the survival function is $S(t_1)$ at time $t_1$ and $S(t_2)$ at time
#' $t_2$ with $t_2>t_1$ and that the hazard rate in the interval $[t_1,t_2]$ is
#' constant. Then, the hazard rate in the interval can be derived as
#' $(\log(S(t_1))-\log(S(t_2)))/(t_2-t_1)$.
#'
#' **Example:** Assume that it is known that the survival function is 1, 0.9, 0.7,
#' and 0.5 at months 0, 12, 24, 48 in the control arm. Then, an interpolating
#' piecewise exponential distribution can be derived as follows:
## ---- include=TRUE, echo=TRUE--------------------------------------------
t <- c(0,12,24,48) # time points at which survival function is known (must include 0)
S <- c(1,0.9,0.7,0.5) # Survival function at timepoints t
# derive hazard in each intervals per the formula above
lambda <- -diff(log(S))/diff(t)
# Define parameters for piecewise exponential distribution in the control arm
# interpolating the targeted survival values
# (code for lambda1 below assumes that the hazard afer month 48 is
# identical to the hazard in interval [24,48])
piecewiseSurvivalTime <- t
lambda2 <- c(lambda,tail(lambda,n=1))
lambda2 # print hazard rates
#'
#' ### Appendix: Random numbers, cumulative distribution, and quantiles for the piecewise exponential distribution {.unnumbered}
#'
#' rpact also provides general functionality for the piecewise exponential distribution,
#' see `?getPiecewiseExponentialDistribution` for details. Below is some example
#' code which shows that the derivation of the piecewise exponential distribution
#' in the previous example is correct.
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
# plot the piecewise exponential survival distribution from the example above
tp <- seq(0,72,by=0.01)
plot(tp,
1-getPiecewiseExponentialDistribution(time=tp,
piecewiseSurvivalTime = piecewiseSurvivalTime,
piecewiseLambda = lambda2),
type="l",xlab="Time (months)",ylab="Survival function S(t)",lwd=2,ylim=c(0,1),axes=F,
main="Piecewise exponential distribution for the example")
axis(1,at=seq(0,72,by=12)); axis(2,at=seq(0,1,by=0.1))
abline(v=seq(0,72,by=12),h=seq(0,1,by=0.1),col=gray(0.9))
# Calculate median survival for the example (which should give 48 months here)
getPiecewiseExponentialQuantile(0.5,piecewiseSurvivalTime = piecewiseSurvivalTime,
piecewiseLambda = lambda2)
#'
#' # Specifying dropout in rpact
#'
#' rpact models dropout with an independent exponentially distributed variable.
#' Dropout is specified by giving the probability of dropout at a specific time
#' point. For example, an annual (12-monthly) dropout probability of 5% in both
#' treatment arms can be specified through the following arguments:
#' ```
#' dropoutRate1 = 0.05, dropoutRate2 = 0.05, dropoutTime = 12
#' ```
#'
#' # Specifying accrual in rpact
#'
#' rpact allows to specify arbitrarily complex recruitment scenarios. Some examples
#' are provided below.
#'
#' ## Absolute accrual intensity and accrual duration known {#accrualDurationKnown}
#'
#' **Example 1:** A **constant accrual intensity** of 24 subjects/months over 30
#' months (i.e., a maximum sample size of 24*30=720 subjects) can be specified
#' through the arguments below:
#' ```
#' accrualTime = c(0,30), accrualIntensity=24
#' ```
#'
#' **Example 2:** A **flexible accrual intensity** of 20 subjects/months in the
#' first 6 months, 25 subjects/months for months 6-12, and 30 subjects per month
#' for months 12-24 can be specified through either of the two equivalent options
#' below:
#'
#'
#' Option 1: List-based definition:
#' ```
#' accrualTime = list("0 - <6" = 20,
#' "6 - <12" = 25,
#' "12- <24" = 30)
#' ```
#'
#' Option 2: Vector-based definition (note that the length of the `accrualTime`
#' vector is 1 larger than the length of the `accrualIntensity` as the end time of
#' the accrual period is also included):
#' ```
#' accrualTime = c(0,6,12,24), accrualIntensity=c(20,25,30)
#' ```
#'
#' ## Absolute accrual intensity known, accrual duration unknown
#'
#' **Example 1:** A **constant accrual intensity** of 24 subjects/months with
#' unspecified end of recruitment is specified through the arguments below:
#' ```
#' accrualTime = 0, accrualIntensity=24
#' # Note: accrualTime is the start time of accrual here which must be explicitly set to 0
#' ```
#'
#' **Example 2:** A **flexible accrual intensity** of 20 subjects/months in the
#' first 6 months, 25 subjects/months from month 6-12, and 30 subjects thereafter
#' can be specified through either of the two equivalent options below:
#'
#'
#' Option 1: List-based definition:
#' ```
#' accrualTime = list("0 - <6" = 20,
#' "6 - <12" = 25,
#' ">=12" = 30)
#' ```
#'
#' Option 2: Vector-based definition (note that the length of the `accrualTime`
#' vector is equal to the length of the accrualIntensity vector as only the start
#' time of the last accrual intensity is provided):
#' ```
#' accrualTime = c(0,6,12), accrualIntensity=c(20,25,30)
#' ```
#'
#' # Sample size calculation for superiority trials without interim analyses
#'
#' Given the specification of the survival distributions, drop-out and accrual as
#' described above, sample size calculations can be readily performed in rpact
#' using the function `getSampleSizeSurvival()`. Importantly, in survival trials,
#' the **number of events (and not the sample size) determines the power** of the
#' trial. Hence, the **choice of the sample size** to reach this target number of
#' events **is based on study-specific trade-offs** between the costs per recruited
#' patient, study duration, and desired minimal follow-up duration.
#'
#' In practice, **a range of sample sizes need to be explored manually or via
#' suitable graphs to find an optimal trade off**. Given a plausible absolute
#' accrual intensity (usually provided by operations), a simple approach in rpact
#' is to use one or both of the two options below:
#'
#' 1. Perform the calculation by trying multiple maximum sample sizes (argument
#' `maxNumberOfSubjects`).
#' 2. It often makes sense to require a minimal follow-up time (argument
#' `followUpTime`) for all patients in the trial at the time of the analysis. In
#' this case, rpact automatically determines the maximum sample size.
#'
#' ## Example - exponential survival, flexible accrual intensity, no interim analyses {#sampleSizeGroupSeq}
#'
#' - Exponential PFS with a median PFS of 60 months in control
#' (`lambda2=log(2)/60`) and a target hazard ratio of 0.74 (`hazardRatio=0.74`).
#' - Log-rank test at the two-sided 5%-significance level (`sided=2, alpha=0.05`),
#' power 80% (`beta=0.2`).
#' - 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 (`allocationRatioPlanned=1`). This is the default and
#' is thus not explicitly set in the function call below.
#' - Two sample size choices will be initially explored:
#' + A fixed total sample size of 1200 (`maxNumberOfSubjects = 1200`).
#' + Alternatively, the total sample size will be implicitly determined by
#' specifying that every subject must have a minimal follow-up duration of
#' at 12 months at the time of the analysis (`followUpTime=12`).
#'
#' Based on this, the required number of events and timing of interim analyses for
#' the **fixed total sample size** of 1200 can be determined as follows:
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
sampleSize1 <- getSampleSizeSurvival(sided=2,alpha=0.05,beta=0.2,
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)
sampleSize1
#' Thus, the required number of events is
#' `r ceiling(sampleSize1$maxNumberOfEvents)` and the MDD corresponds to an
#' observed HR of
#' `r formatC(sampleSize1$criticalValuesEffectScale,digits=2,format="f")`.The
#' `r ceiling(sampleSize1$maxNumberOfSubjects)` subjects will be recruited over
#' `r formatC(sampleSize1$totalAccrualTime,digits=2,format="f")` months and the
#' total study duration is
#' `r formatC(sampleSize1$studyDuration,digits=2,format="f")` months. The user
#' could now vary `maxNumberSubject` to further optimize the trade-off between
#' sample size and study duration.
#'
#' Alternatively, **specifying a minimum follow-up duration** of 12 months leads to
#' the following result:
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
sampleSize2 <- getSampleSizeSurvival(sided=2,alpha=0.05,beta=0.2,
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),
followUpTime=12)
sampleSize2
#'
#'
#' This specification leads to a higher sample size
#' of `r ceiling(sampleSize2$maxNumberOfSubjects)` subjects which will be recruited
#' over `r formatC(sampleSize2$totalAccrualTime,digits=2,format="f")` months and
#' the total study duration is
#' `r formatC(sampleSize2$studyDuration,digits=2,format="f")` months.
#'
#' With the generic `summary()` function the two calculations are summarized as
#' follows:
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
summary(sampleSize1)
summary(sampleSize2)
#'
#'
#' To further explore the possible trade-offs, one could visualize recruitment and
#' study durations for a range of sample sizes as illustrated in the code below:
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
# set up data frame which contains sample sizes and corresponding durations
sampleSizeDuration <- data.frame(
maxNumberSubjects=seq(600,1800,by=50),
accrualTime=NA,
studyDuration=NA)
# calculate recruitment and study duration for each sample size
for (i in 1:nrow(sampleSizeDuration)){
sampleSizeResult <- getSampleSizeSurvival(
sided = 2,alpha = 0.05,beta = 0.2,
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 = sampleSizeDuration$maxNumberSubjects[i])
sampleSizeDuration$accrualTime[i] <- sampleSizeResult$totalAccrualTime
sampleSizeDuration$studyDuration[i] <- sampleSizeResult$maxStudyDuration
}
# plot result
plot(sampleSizeDuration$maxNumberSubjects,
sampleSizeDuration$studyDuration,type="l",
xlab="Total sample size",
ylab="Duration (months)",
main="Recruitment and study duration vs sample size",
ylim=c(0,max(sampleSizeDuration$studyDuration)),
col="blue",lwd=1.5)
lines(sampleSizeDuration$maxNumberSubjects,
sampleSizeDuration$accrualTime,col="red",lwd=1.5)
legend(x=1000,y=100,
legend=c("Study duration under H1",
"Recruitment duration"),
col=c("blue","red"),lty=1,lwd=1.5)
#'
#' # Sample size calculation for non-inferiority trials without interim analyses
#'
#' Sample size calculation proceeds in the same fashion as for superiority trials
#' except that the role of the null and the alternative hypothesis are reversed.
#' I.e., in this case, the non-inferiority margin $\Delta$ corresponds to the
#' treatment effect under the null hypothesis (`thetaH0`). Testing for
#' non-inferiority trials is always one-sided.
#'
#' In the example code below, we assume that the non-inferiority margin corresponds
#' to a 20% increase of the hazard function and that "in reality", i.e., under the
#' alternative hypothesis, the survival functions in the two groups are equal. In
#' this case, the "null hypothesis" is `thetaH0=1.2` versus the "alternative"
#' `hazardRatio=1`. For illustration, all other assumptions were chosen as per the
#' sample size calculations for superiority trials above.
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
sampleSizeNonInf <- getSampleSizeSurvival(sided = 1,alpha = 0.025,beta = 0.2,
lambda2 = log(2)/60,
thetaH0=1.2, hazardRatio = 1,
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),
followUpTime = 12)
sampleSizeNonInf
#' In this example, the required number of events is
#' `r ceiling(sampleSizeNonInf$maxNumberOfEvents)`. The specification above leads
#' to a sample size of `r ceiling(sampleSizeNonInf$maxNumberOfSubjects)` subjects
#' recruited over
#' `r formatC(sampleSizeNonInf$totalAccrualTime,digits=2,format="f")` months and an
#' expected total study duration of
#' `r formatC(sampleSizeNonInf$studyDuration,digits=2,format="f")` months.
#'
#' # Sample size calculation for trials with interim analyses
#'
#' Calculations can be performed as for trials without interim analyses except that
#' the group-sequential design needs to be specified using
#' `getDesignGroupSequential` before calling `getSampleSizeSurvival`. For details
#' regarding the specification of the group-sequential design, see the R Markdown
#' file "Defining group-sequential boundaries with rpact".
#'
#' In general, rpact supports both one-sided and two-sided group-sequential
#' designs. However, if futility boundaries are specified, only one-sided tests are
#' permitted. **For simplicity, it is often preferred to use one-sided tests for
#' group-sequential designs** (typically with $\alpha=0.025$).
#'
#' One additional consideration for phase III trials with **efficacy interim
#' analyses** is that they **should only be performed at time-points when the data
#' is sufficiently mature** to allow for a potential filing. For example, it is
#' often sensible to perform efficacy interim analyses only after all subjects have
#' been recruited. This aspect can be further explored by looking at simulated data
#' for (hypothetical) trials that would have stopped at the interim analysis. This
#' can easily be implemented using the **simulation tool** in rpact and we refer to
#' the R Markdown document [Simulation-based design of group-sequential trials with
#' a survival endpoint with
#' rpact](https://vignettes.rpact.org/html/rpact_survival_simulation_examples.html)
#' for further details.
#'
#' ## Example - piece-wise exponential survival, constant accrual intensity, group-sequential design
#'
#' - 1:1 randomized.
#' - 80% power at the one-sided 2.5% significance level.
#' - Efficacy interim analyses at 50% and 75% of total information using an
#' O'Brien&Fleming type $\alpha$-spending function, no futility interim analyses.
#' - Target HR for primary endpoint (PFS) is 0.75.
#' - PFS in the control arm follows a piece-wise exponential distribution, with
#' the hazard rate h(t) estimated using historical controls as follows:
#' + h(t) = 0.025 for t between 0 and 6 months,
#' + h(t) = 0.04 for t between 6 and 9 months,
#' + h(t) = 0.015 for t between 9 and 15 months,
#' + h(t) = 0.01 for t between 15 and 21 months,
#' + h(t) = 0.007 for t beyond 21 months.
#' - An annual dropout probability of 5%.
#' - Recruitment of 42 patients/month up to a maximum of 1000 subjects.
#'
#' ### Sample size calculation
#' The **group-sequential design** (including Type I and II error) can be specified
#' as follows:
## ---- include=TRUE, echo=TRUE--------------------------------------------
design <- getDesignGroupSequential(sided = 1, alpha = 0.025, beta = 0.2,
informationRates = c(0.5,0.75,1),
typeOfDesign="asOF")
#'
#' The **piecewise exponential survival** is defined as described above:
## ---- include=TRUE, echo=TRUE--------------------------------------------
piecewiseSurvivalTime <- list(
"0 - <6" = 0.025,
"6 - <9" = 0.04,
"9 - <15" = 0.015,
"15 - <21" = 0.01,
">=21" = 0.007)
#'
#' Now, the **sample size** calculation can be performed:
## ---- include=TRUE, echo=TRUE--------------------------------------------
sampleSizeResult <- getSampleSizeSurvival(
design = design,
piecewiseSurvivalTime = piecewiseSurvivalTime,hazardRatio = 0.75,
dropoutRate1 = 0.05, dropoutRate2 = 0.05, dropoutTime = 12,
accrualTime = 0, accrualIntensity=42, maxNumberOfSubjects = 1000)
#' The sample size result can be summarized by displaying the **default rpact output**:
## ---- include=TRUE, echo=TRUE--------------------------------------------
sampleSizeResult
#' In addition, the **generic `summary()` function** summarizes the results per
#' stage in a different format and also adds results which are stored in the
#' underlying group-sequential `design` such as exit probabilities:
## ---- include=TRUE, echo=TRUE--------------------------------------------
# Summarize design
summary(sampleSizeResult)
#'
#' ### Sample size plots
#'
#' rpact provides a large number of plots, see the R Markdown document [How to
#' create admirable plots with
#' rpact](https://vignettes.rpact.org/html/rpact_plot_examples.html)
#' for examples of all plots. Below, we just show two selected plots.
## ---- include=TRUE, echo=TRUE--------------------------------------------
# Boundary on treatment effect scale
plot(sampleSizeResult, type = 2, showSource = TRUE)
# Survival function and piecewise constant hazard
plot(sampleSizeResult, type = 14)
#'
#' # Power calculation
#'
#' The function `getPowerSurvival()` is the analogue of `getSampleSizeSurvival()`
#' for power calculations and has a very similar syntax. We illustrate its usage
#' for the above example below.
#'
#' ## Example - piece-wise exponential survival, constant accrual intensity, group-sequential design (continued)
#'
#' What is the power of the trial from the example above if the true hazard ratio
#' is 0.7 rather than the targeted hazard ratio of 0.75 above?
#'
#' The corresponding power calculation is given below. The function call uses the
#' same arguments as the call to function `getSampleSizeSurvival()` above except
#' that the maximum sample size of 1000 and the maximum number of events of 387
#' from above are provided as arguments and the hazard ratio is changed to 0.7. In
#' addition, the direction of the test (`directionUpper`) needs to be specified.
## ---- include=TRUE, echo=TRUE--------------------------------------------
powerResult <- getPowerSurvival(
design = design,
piecewiseSurvivalTime = piecewiseSurvivalTime,hazardRatio = 0.7,
dropoutRate1 = 0.05, dropoutRate2 = 0.05, dropoutTime = 12,
accrualTime = 0, accrualIntensity=42,
maxNumberOfSubjects = 1000,maxNumberOfEvents=387,directionUpper = FALSE)
powerResult
#' This yields an overall power of
#' `r formatC(powerResult$overallReject,digits=3,format="f")`. Note also that the
#' timing of the interim analyses is slightly deferred now (as the events now occur
#' more slowly in the intervention arm) but the expected study duration is
#' decreased from `r formatC(sampleSizeResult$studyDurationH1,digits=2,format="f")`
#' to `r formatC(powerResult$studyDuration,digits=2,format="f")` months because
#' chances to stop the trial early have increased.
#'
#' The impact of the true hazard ratio on power and study duration can also be
#' calculated by entering the `hazardRatio` argument as a vector. Below, the hazard
#' ratio is varied from 0.6 to 1 (in increments of 0.02) and the results are
#' displayed graphically.
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
powerResult2 <- getPowerSurvival(
design = design,
piecewiseSurvivalTime = piecewiseSurvivalTime,hazardRatio = seq(0.6,1,by=0.02),
dropoutRate1 = 0.05, dropoutRate2 = 0.05, dropoutTime = 12,
accrualTime = 0, accrualIntensity=42,
maxNumberOfSubjects = 1000,maxNumberOfEvents=387,directionUpper = FALSE)
# Plot true hazard ratio vs overall power, probability of an early rejection,
# and expected number of events
plot(powerResult2, type = 6, legendPosition = 3) # legend in left bottom
# Plot hazard ratio vs analysis time
plot(powerResult2, type = 12, legendPosition = 4) # right top
#'
#' For examples of all available plots in rpact, see the R vignette "How to create
#' admirable plots with rpact" (written by rpact,
#' [link](https://vignettes.rpact.org/html/rpact_plot_examples.html)).
#'
#' # Prediction of number of events over time
#'
#' The function `getEventProbabilities` calculates the probability that a randomly
#' chosen subject has an observed event up to a specied `time` after the start of
#' the study. By multiplying this probability with the maximal sample size, we
#' obtain the **expected number of events** up to that time point.
#'
#' The arguments for the function `getEventProbabilities` are similar to those for
#' the function `getSampleSizeSurvival` discussed above. One difference is that the
#' maximal sample size must always be provided in function `getEventProbabilities`.
#' This can either be done explicitly by setting argument `maxNumberOfSubjects` or
#' implicitly by specifying the accrual duration as described in Section
#' \@ref(accrualDurationKnown).
#'
#' For illustration, we revisit the example from Section \@ref(sampleSizeGroupSeq).
#' In that example, the sample size was `r sampleSizeResult$maxNumberOfSubjects`
#' and the expected timing of the final analysis after observing
#' `r ceiling(sampleSizeResult$maxNumberOfEvents)` events was after
#' `r formatC(sampleSizeResult$maxStudyDuration,digits=2,format="f")` months. The
#' code below re-calculates the expected number of events until that time point.
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
# Sample size per arm and overall for the example
maxNumberOfSubjects1 <- 500
maxNumberOfSubjects2 <- 500
maxNumberOfSubjects <- maxNumberOfSubjects1+maxNumberOfSubjects2
# Probability that a randomly chosen subject has an observed event
# until 60 months after study start
probEvent <- getEventProbabilities(time = 60,
piecewiseSurvivalTime = piecewiseSurvivalTime,hazardRatio = 0.75,
dropoutRate1 = 0.05, dropoutRate2 = 0.05, dropoutTime = 12,
accrualTime = 0, accrualIntensity = 42,maxNumberOfSubjects = maxNumberOfSubjects)
probEvent
# Total expected number of events until 60 months
probEvent$overallEventProbabilities*maxNumberOfSubjects
# Expected number of events per group until 60 months
probEvent$eventProbabilities1*maxNumberOfSubjects1 # intervention
probEvent$eventProbabilities2*maxNumberOfSubjects2 # control
#'
#' The function `getEventProbabilities()` also allows a vector argument for `time`
#' and the code below illustrates how to calculate (and plot) expected event
#' numbers over time (overall and by group). In addition, the number of recruited
#' subjects is added to the plot which can be obtained with function
#' `getNumberOfSubjects()`.
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
time <- 1:60
# Expected number of events overall and per group at 1-60 months after study start
probEvent <- getEventProbabilities(time = time,
piecewiseSurvivalTime = piecewiseSurvivalTime,hazardRatio = 0.75,
dropoutRate1 = 0.05, dropoutRate2 = 0.05, dropoutTime = 12,
accrualTime = 0, accrualIntensity = 42,maxNumberOfSubjects = maxNumberOfSubjects)
expEvent <- probEvent$overallEventProbabilities*maxNumberOfSubjects # overall
expEventInterv <- probEvent$eventProbabilities1*maxNumberOfSubjects1 # intervention
expEventCtrl <- probEvent$eventProbabilities2*maxNumberOfSubjects2 # control
# Cumulative number of recruited patients over time
nSubjects <- getNumberOfSubjects(time=time, accrualTime = 0, accrualIntensity=42,
maxNumberOfSubjects = maxNumberOfSubjects)
# plot result
plot(time,nSubjects$numberOfSubjects,type="l",lwd=2,
main="Number of patients and expected number of events over time",
xlab="Time from study start (months)",
ylab="Number of patients / events")
lines(time,expEvent,col="blue",lwd=2)
lines(time,expEventCtrl,col="green")
lines(time,expEventInterv,col="orange")
legend(x=28,y=900,
legend=c("Recruitment","Expected events - All subjects",
"Expected events - Control","Expected events - Intervention"),
col=c("black","blue","green","orange"),lty=1,lwd=c(2,2,1,1))
#'
#' # Interim analyses for multiple time-to-event endpoints
#'
#' The function `getEventProbabilities()` can also be used to align interim
#' analyses of different time to event endpoints as illustrated in two examples
#' below.
#'
#' ## Example: PFS and OS as co-primary endpoints, OS interim aligned with PFS final analysis
#'
#' - Overall Type I error 5% (two-sided), 1:1 randomized with PFS and OS as
#' co-primary endpoints.
#' - Assumptions for PFS (1% Type I error assigned):
#' + Median PFS control is 6 months, target hazard ratio is 0.65, 95% power at
#' two-sided 1% significance level.
#' + Only one primary analysis.
#' - Assumptions for OS (4% Type I error assigned):
#' + Median OS control is 12 month, target hazard ratio is 0.75, 80% power at
#' two-sided 4% significance level.
#' + Interim analysis at the primary PFS analysis, final analysis when the
#' target number of OS events has been reached.
#' + O'Brien&Fleming type $\alpha$-spending function.
#' - Uniform recruitment of 60 patients/month over 10 months, no dropout.
#'
#' To display results more concisely, all outputs below were created with the
#' generic `summary()` function.
#'
#' The first step is to calculate the **required number of events and the timing of the primary analyses for PFS**:
## ---- include=TRUE, echo=TRUE--------------------------------------------
sampleSizePFS <- getSampleSizeSurvival(beta = 0.05,sided = 2,alpha = 0.01,
lambda2 = log(2)/6,hazardRatio = 0.65,
accrualTime = c(0,10),accrualIntensity = 60)
summary(sampleSizePFS)
#'
#' Based on this, the estimated time of the primary PFS analysis is after
#' `r formatC(sampleSizePFS$analysisTime,2,format="f")` months.
#'
#' Second, we **calculate the expected number of OS events at the final PFS
#' analysis** after `r formatC(sampleSizePFS$analysisTime,2,format="f")` months.
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
# Expected number of OS events at the time of the final PFS analysis
probOS <- getEventProbabilities(time = 16.37,
accrualTime = c(0,10),accrualIntensity = 60,
lambda2 = log(2)/12, hazardRatio = 0.75)
probOS$overallEventProbabilities*600
#'
#' Third, we caculate the **required number of events and the timing of analyses
#' for OS using an initial guess for the information fraction at the OS interim
#' analysis** of 50%:
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
# Preliminary group-sequential design for OS assuming interim at 50% information
designOS1 <- getDesignGroupSequential(sided = 2, alpha = 0.04,beta = 0.2,
informationRates = c(0.5, 1), typeOfDesign = "asOF")
# Preliminary sample size for OS assuming interim at 50% information
sampleSizeOS1 <- getSampleSizeSurvival(designOS1,
lambda2 = log(2)/12,hazardRatio = 0.75,
accrualTime = c(0,10),accrualIntensity = 60)
summary(sampleSizeOS1)
#'
#' Finally, we **recalculate** the required number of events and timing of analyses
#' for OS using **an updated estimate of the information fraction at the interim
#' analysis**. Specifically, we calculate this fraction as the expected number of
#' events at the time of the final PFS analysis (i.e.,
#' `r ceiling(probOS$overallEventProbabilities*600)` events) relative to the
#' maximal number of events based on the calculation for OS above (i.e.,
#' `r ceiling(sampleSizeOS1$maxNumberOfEvents)`):
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
# Updated group-sequential design for OS assuming an updated information fraction
designOS2 <- getDesignGroupSequential(sided = 2,alpha = 0.04,beta = 0.2,
informationRates = c(258/407, 1), typeOfDesign = "asOF")
# Updated sample size for OS
sampleSizeOS2 <- getSampleSizeSurvival(designOS2,
lambda2 = log(2)/12,hazardRatio = 0.75,
accrualTime = c(0,10),accrualIntensity = 60)
summary(sampleSizeOS2)
#'
#' According to the updated calculation, the panned timing of the OS interim
#' analysis is after `r formatC(sampleSizeOS2$analysisTime[1,1],2,format="f")`
#' months which is sufficiently close to the planned timing of the PFS final
#' analysis after `r formatC(sampleSizePFS$analysisTime,2,format="f")` months for planning purposes.
#'
#' ## Example: Power considerations for OS as a secondary endpoint, OS interim aligned with primary PFS analysis
#'
#' - Overall Type I error 5% (two-sided), 1:1 randomized with PFS as the primary
#' endpoint.
#' - Sample size calculation for PFS resulted in the following design:
#' + No interim analyses for PFS.
#' + Uniform recruitment of 80 patients/month over 10 months (800 patients in
#' total), no dropout.
#' + Final PFS analysis with 80% power for the targeted HR is expected after a
#' study duration of 24 months.
#' - Asssumption for OS:
#' + Median OS 30 months, OS HR expected to be around 0.8.
#' + OS is only formally tested if PFS is significant (hierarchical testing).
#' + One interim analysis at the time of the final PFS analysis.
#' + It is anticipated that the OS analysis cannot be fully powered. The target
#' number of OS events at the final OS analysis will be pre-defined after
#' considering all possible trade-offs between study power and duration.
#' + O'Brien&Fleming type $\alpha$-spending function.
#'
#' The code below calulates study duration and power for the final OS analysis
#' (conditional on significance for the primary PFS analysis) depending on the
#' target number of OS events at the final analysis. The resulting plots should
#' facilitate the final choice of the timing of the OS analysis.
#'
## ---- include=TRUE, echo=TRUE--------------------------------------------
# 1) Calculate expected number of OS events at the time of the final PFS analysis
noDeathsAtPFSAnalysis <- getEventProbabilities(time = 24,
lambda2 = log(2)/30, hazardRatio = 0.8,
accrualTime = c(0,20),accrualIntensity = 40)$overallEventProbabilities*800
noDeathsAtPFSAnalysis
# 2) Set up data frame which will contain number of OS events and
# corresponding power and study durations
# (vary target number of OS events from 250-700)
powerOSevents <- data.frame(
maxNumberOSEvents=seq(250,700,by=10),
power=NA,
durationFinalOSAnalysis=NA)
# 3) Calculate power and study durations depending on target number of OS events
for (i in 1:nrow(powerOSevents)){
# Define OF-type design, with first interim after noDeathsAtPFSAnalysis events
# and final analysis after maxNumberOSEvents[i] events
# PS: beta not explicitly set here because boundary and
# subsequent power calculations do not depend on it
informationRates <- c(noDeathsAtPFSAnalysis,powerOSevents$maxNumberOSEvents[i])/
powerOSevents$maxNumberOSEvents[i]
design <- getDesignGroupSequential(sided = 2, alpha = 0.05,typeOfDesign= "asOF",
informationRates = informationRates, twoSidedPower = TRUE)
# Power and study duration after maxNumberOSEvents[i]
powerResult <- getPowerSurvival(
design = design,
lambda2 = log(2)/30, hazardRatio = 0.8,
accrualTime = c(0,20),accrualIntensity = 40,
maxNumberOfEvents = powerOSevents$maxNumberOSEvents[i])
powerOSevents$power[i] <- powerResult$overallReject
powerOSevents$durationFinalOSAnalysis[i] <- powerResult$maxStudyDuration
}
# 4) Plot results
par(mfrow=c(1,2))
# events cs power
plot(powerOSevents$maxNumberOSEvents,
powerOSevents$power,type="l",
xlab="Required # OS events at final analysis",
ylab="Power",ylim=c(0,1),lwd=1.5)
# Duration vs power
plot(powerOSevents$durationFinalOSAnalysis,
powerOSevents$power,type="l",
xlab="Time until final OS analysis (months)",
ylab="Power",ylim=c(0,1),lwd=1.5)
#'
#' # Extracting information from rpact objects
#' `names(rpact_object)` gives variable names of an rpact object which can be
#' accessed via
#' `rpact_object$varname`. For example, `sampleSizeResult$studyDurationH1` contains
#' the expected study duration under H1. For more details on using R generics for
#' rpact objects, please refer to vignette
#' ["Example: using rpact with R
#' generics"](https://vignettes.rpact.org/html/rpact_generic_functions.html)).
#'
#'
#'
#' ***
#'
#' System: rpact `r packageVersion("rpact")`, `r R.version.string`, platform: `r R.version$platform`
#'
## ---- include=TRUE, echo=FALSE, results='asis'---------------------------
printCitation()
#'