This R Markdown document provides examples for designing trials with survival endpoints using rpact.
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
## [1] '3.2.1'
The sample size calculation for a group-sequential trial with a survival endpoints is performed in two steps:
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.
Calculate sample size for the survival endpoint by feeding the abstract
design into the 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 the sample size.getEventProbabilities()
: Calculates the probability of an event depending on
the time and type of accrual, follow-up time, and survival distribution.
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.
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.
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.
eventTime
.pi2
.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
.
lambda2
.lambda1
or, alternatively, implicitly by specifying the target
hazard ratio hazardRatio
.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 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.
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()
.
Piecewise exponential distributions 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:
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
## [1] 0.008780043 0.020942869 0.014019677 0.014019677
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.
# 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)
## [1] 48
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
rpact allows to specify arbitrarily complex recruitment scenarios. Some examples are provided below.
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)
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 of accrual 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)
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:
maxNumberOfSubjects
).followUpTime
) for all patients in the trial at the time of the analysis. In
this case, rpact automatically determines the maximum sample size.lambda2 = log(2)/60
) and a target hazard ratio of 0.74 (hazardRatio = 0.74
).sided = 2, alpha = 0.05
),
power 80% (beta = 0.2
).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)
)allocationRatioPlanned = 1
). This is the default and
is thus not explicitly set in the function call below.maxNumberOfSubjects = 1200
).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:
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
## Design plan parameters and output for survival data:
##
## Design parameters:
## Critical values : 1.96
## Two-sided power : FALSE
## Significance level : 0.0500
## Type II error rate : 0.2000
## Test : two-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
## Number of events : 346.3
## Total accrual time : 31.57
## Follow up time : 21.54
## Number of events fixed : 346.3
## Number of subjects fixed : 1200.0
## Number of subjects fixed (1) : 600.0
## Number of subjects fixed (2) : 600.0
## Analysis times : 53.11
## Study duration : 53.11
## Lower critical values (treatment effect scale) : 1.234
## Upper critical values (treatment effect scale) : 0.810
## Local two-sided significance levels : 0.0500
##
## Legend:
## (i): values of treatment arm i
Thus, the required number of events is
347 and the MDD corresponds to an
observed HR of
.The
1200 subjects will be recruited over
31.57 months and the
total study duration is
53.11 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:
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
## Design plan parameters and output for survival data:
##
## Design parameters:
## Critical values : 1.96
## Two-sided power : FALSE
## Significance level : 0.0500
## Type II error rate : 0.2000
## Test : two-sided
##
## User defined parameters:
## lambda(2) : 0.0116
## Hazard ratio : 0.740
## Accrual intensity : 6, 12, 18, 24, 30, 36, 42
## Follow up time : 12.00
## 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 : 1433.7
## Number of events : 346.3
## Accrual time : 1.00, 2.00, 3.00, 4.00, 5.00, 6.00, 37.13
## Total accrual time : 37.13
## Number of events fixed : 346.3
## Number of subjects fixed : 1433.7
## Number of subjects fixed (1) : 716.8
## Number of subjects fixed (2) : 716.8
## Analysis times : 49.13
## Study duration : 49.13
## Lower critical values (treatment effect scale) : 1.234
## Upper critical values (treatment effect scale) : 0.810
## Local two-sided significance levels : 0.0500
##
## Legend:
## (i): values of treatment arm i
This specification leads to a higher sample size of 1434 subjects which will be recruited over 37.13 months and the total study duration is 49.13 months.
With the generic summary()
function the two calculations are summarized as
follows:
## Sample size calculation for a survival endpoint
##
## Fixed sample analysis, significance level 5% (two-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,
## 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 Fixed
## Efficacy boundary (z-value scale) 1.960
## Number of subjects 1200.0
## Number of events 346.3
## Expected study duration 53.1
## Two-sided local significance level 0.0500
## Efficacy boundary (t) 1.234 - 0.810
##
## Legend:
## (t): treatment effect scale
## Sample size calculation for a survival endpoint
##
## Fixed sample analysis, significance level 5% (two-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,
## accrual time = c(1, 2, 3, 4, 5, 6, 37.135),
## 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 Fixed
## Efficacy boundary (z-value scale) 1.960
## Number of subjects 1433.7
## Number of events 346.3
## Expected study duration 49.1
## Two-sided local significance level 0.0500
## Efficacy boundary (t) 1.234 - 0.810
##
## Legend:
## (t): treatment effect scale
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:
# 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)