library(gsDesign)
fixed_design <- nSurv()
fixed_design
#> nSurv fixed-design summary (method=LachinFoulkes; target=Accrual rate)
#> HR=0.600 vs HR0=1.000 | alpha=0.025 (sided=1) | power=90.0%
#> N=250.4 subjects | D=160.5 events | T=18.0 study duration | accrual=12.0 Accrual duration | minfup=6.0 minimum follow-up | ratio=1 randomization ratio (experimental/control)
#>
#> Key inputs (names preserved):
#> desc item value input
#> Accrual rate(s) gamma 20.871 1
#> Accrual rate duration(s) R 12 12
#> Control hazard rate(s) lambdaC 0.116 0.116
#> Control dropout rate(s) eta 0 0
#> Experimental dropout rate(s) etaE 0 etaE
#> Event and dropout rate duration(s) S NULL S5 Time-to-event design and power
We extend the Lachin and Foulkes (Lachin and Foulkes 1986) method to cases where the null hypothesis does not reflect equality. This includes non-inferiority scenarios. For vaccines or other prevention studies this also includes super-superiority. Denote the null hypothesis failure rates for control and experimental treatment groups as \(\lambda_{00}\) and \(\lambda_{01}\), respectively. Denote the alternate hypothesis rates as \(\lambda_{10}\) and \(\lambda_{11}\). Further, denote the alternate hypothesis hazard ratio \(h_1=\lambda_{11}/\lambda_{10}\), and the null hypothesis hazard ratio \(h_0=\lambda_{01}/\lambda_{00}\). We let censoring rates be specific to the control (\(\eta_0\)) and experimental (\(\eta_1\)) groups; these values are only implicit in the equations below. Further, we let \(\xi\) denote the proportion of subjects randomized to the experimental treatment group. Finally, we let \(\eta\) represent an exponential dropout rate independent and the time to dropout is independent of the time until an event. Lachin and Foulkes assumed a null hypothesis with no difference between failure rates in the control and experimental rates and test for superiority. That is, \(\lambda_{00}=\lambda_{01}\) (\(h_0=1\)) and \(\lambda_{10}<\lambda_{01}\) (\(h_1<1\)). They set event rates under the null hypothesis so that the the weighted average event rate is the same under the null and alternate hypotheses:
\[ \lambda_{00}=\lambda_{01}=\bar\lambda=(1-\xi)\lambda_{10}+\xi\lambda_{11}. \tag{5.1}\]
The apparent intent of this is to equalize the variance for the log hazard ratio under null and alternative hypotheses; this will not be exactly the case. We let \(\delta\) represent an indicator that an uncensored event is observed for a patient in a specified treatment group given enrollment, event rate, and dropout rate assumptions. The Lachin and Foulkes power equation for proportional hazards translates in our notation to:
\[ \begin{split} \sqrt{N}\ln(h_1) &= Z_\alpha \sqrt{E\left\{\delta|\bar{\lambda},\eta\right\}^{-1}(\xi^{-1}+(1-\xi)^{-1})} \\ &+Z_\beta\sqrt{ E\left\{\delta|\lambda_1,\eta\right\}^{-1}\xi^{-1}+ E\left\{\delta|\lambda_0,\eta\right\}^{-1}(1-\xi)^{-1}} \end{split} \tag{5.2}\]
Lachin and Foulkes did not cover any cases other than equality under the null hypothesis; i.e., the assumed \(h_0 \neq 1\) (i.e., \(\lambda_{00} \neq \lambda_{01}\)). Equation 5.2 generalizes in this case to
\[ \begin{split} \sqrt{N}\ln\left(\frac{h_1}{h_0}\right) &= Z_\alpha \sqrt{E\left\{\delta|\lambda_{01},\eta_1\right\}^{-1}\xi^{-1}+ E\left\{\delta|\lambda_{00},\eta_0\right\}^{-1}(1-\xi)^{-1}} \\ &+Z_\beta\sqrt{ E\left\{\delta|\lambda_{11},\eta_1\right\}^{-1}\xi^{-1}+ E\left\{\delta|\lambda_{10},\eta_0\right\}^{-1}(1-\xi)^{-1}} \end{split} \tag{5.3}\]
While we have defined null hypothesis assumptions \(\lambda_{00}\) and \(\lambda_{01}\) for an exponential distribution, the gsDesign functions nSurv() and gsSurv() extend the approach above in an analogous fashion to piecewise exponential failure and dropout rates with a common proportional hazard ratio across piecewise intervals.
For a fixed sample size with default arguments, we have:
This intentionally does not round up, so the user needs to round the number of events and sample size up. For gsSurv(), this rounding can be done automatically:
gs_design <- gsSurv() |> toInteger()
gs_design |> gsBoundSummary()
#> Analysis Value Efficacy Futility
#> IA 1: 33% Z 3.0139 -0.2458
#> N: 190 p (1-sided) 0.0013 0.5971
#> Events: 57 ~HR at bound 0.4500 1.0673
#> Month: 8 P(Cross) if HR=1 0.0013 0.4029
#> P(Cross) if HR=0.6 0.1396 0.0147
#> IA 2: 67% Z 2.5428 0.9509
#> N: 270 p (1-sided) 0.0055 0.1708
#> Events: 115 ~HR at bound 0.6224 0.8375
#> Month: 13 P(Cross) if HR=1 0.0063 0.8370
#> P(Cross) if HR=0.6 0.5853 0.0440
#> Final Z 1.9994 1.9994
#> N: 270 p (1-sided) 0.0228 0.0228
#> Events: 172 ~HR at bound 0.7372 0.7372
#> Month: 18 P(Cross) if HR=1 0.0233 0.9767
#> P(Cross) if HR=0.6 0.9003 0.0997A textual summary is also available:
summary(gs_design)#> Asymmetric two-sided group sequential design with
#> non-binding futility bound, 3 analyses, time-to-event
#> outcome with sample size 270 and 172 events required, 90
#> percent power, 2.5 percent (1-sided) Type I error to detect
#> a hazard ratio of 0.6. Enrollment and total study durations
#> are assumed to be 12 and 18 months, respectively. Efficacy
#> bounds derived using a Hwang-Shih-DeCani spending function
#> with gamma = -4. Futility bounds derived using a
#> Hwang-Shih-DeCani spending function with gamma = -2.All the assumptions laid out in this text can be changed as documented in the help file.
5.1 Schoenfeld and Freedman methods
In addition to the Lachin and Foulkes method, the nSurv() and gsSurv() functions support the Schoenfeld and Freedman methods for sample size and power for time-to-event studies under proportional hazards. The method is selected via the method parameter:
-
method = "LF"(default): Lachin and Foulkes method, which accounts for the enrollment pattern and follow-up duration. -
method = "Schoenfeld": Schoenfeld approximation based on number of events. -
method = "Freedman": Freedman method.
These alternative methods can be useful for matching sample size calculations with other software packages:
# Compare methods for the cancer trial example
lf <- nSurv(
lambdaC = log(2) / 6, hr = 0.7, eta = -log(0.95) / 12,
T = 36, minfup = 6
)
sch <- nSurv(
lambdaC = log(2) / 6, hr = 0.7, eta = -log(0.95) / 12,
T = 36, minfup = 6, method = "Schoenfeld"
)
fr <- nSurv(
lambdaC = log(2) / 6, hr = 0.7, eta = -log(0.95) / 12,
T = 36, minfup = 6, method = "Freedman"
)
cat(paste(
"Lachin-Foulkes events:", ceiling(lf$d),
"\nSchoenfeld events: ", ceiling(sch$d),
"\nFreedman events: ", ceiling(fr$d), "\n"
))
#> Lachin-Foulkes events: 330
#> Schoenfeld events: 331
#> Freedman events: 3385.2 Calendar-based design timing
The gsSurvCalendar() function provides an alternative to gsSurv() for designing group sequential trials with calendar-time-based interim analysis timing. While gsSurv() specifies interim analyses by information fractions (proportions of final planned events), gsSurvCalendar() specifies analyses at fixed calendar times. This is often more natural for trial planning, where interim analyses are scheduled at specific time points (e.g., annually).
# Calendar-based design for the cancer trial
xcal <- gsSurvCalendar(
calendarTime = c(18, 24, 36),
alpha = 0.025,
beta = 0.1,
lambdaC = log(2) / 6,
hr = 0.7,
eta = -log(0.95) / 12,
R = 30,
minfup = 6
) |> toInteger()
gsBoundSummary(xcal)
#> Analysis Value Efficacy Futility
#> IA 1: 39% Z 2.9172 0.0077
#> N: 268 p (1-sided) 0.0018 0.4969
#> Events: 137 ~HR at bound 0.6075 0.9987
#> Month: 18 P(Cross) if HR=1 0.0018 0.5031
#> P(Cross) if HR=0.7 0.2045 0.0186
#> IA 2: 60% Z 2.6745 0.6826
#> N: 356 p (1-sided) 0.0037 0.2474
#> Events: 211 ~HR at bound 0.6920 0.9103
#> Month: 24 P(Cross) if HR=1 0.0047 0.7697
#> P(Cross) if HR=0.7 0.4815 0.0366
#> Final Z 1.9928 1.9928
#> N: 444 p (1-sided) 0.0231 0.0231
#> Events: 350 ~HR at bound 0.8081 0.8081
#> Month: 36 P(Cross) if HR=1 0.0235 0.9765
#> P(Cross) if HR=0.7 0.9003 0.0997The calendarTime parameter specifies calendar times (in months from start of enrollment) at which analyses are planned. The function computes the expected number of events at each time point and derives the design accordingly.
5.3 Bidirectional hazard ratio support
The survival design functions now support hr > hr0 for designs where a larger hazard ratio is the alternative hypothesis. This enables direct specification of designs for:
- Time-to-response endpoints where faster response (higher hazard) is beneficial.
- Safety endpoints where the alternative hypothesis is that harm exceeds a threshold.
- Reversed HR conventions where treatment effect is measured with HR > 1.
All sample size methods (Lachin-Foulkes, Schoenfeld, Freedman) handle both directions symmetrically:
# Standard: HR < 1 (treatment reduces hazard)
x_standard <- nSurv(
lambdaC = log(2) / 6, hr = 0.7, hr0 = 1
)
# Reversed: HR > 1 (treatment increases hazard, e.g., time-to-response)
x_reversed <- nSurv(
lambdaC = log(2) / 6, hr = 1 / 0.7, hr0 = 1
)
cat(paste(
"Standard (HR=0.7) events:", ceiling(x_standard$d),
"\nReversed (HR=1.43) events:", ceiling(x_reversed$d), "\n"
))
#> Standard (HR=0.7) events: 330
#> Reversed (HR=1.43) events: 330
5.4 Power computation with gsSurvPower()
The gsSurvPower() function computes power for a group sequential survival design with specified enrollment, dropout, treatment effect, and analysis timing. Unlike gsSurv() and gsSurvCalendar() which solve for sample size to achieve a target power, gsSurvPower() takes fixed assumptions and computes the resulting power. See the gsSurvPower() vignette for complete argument details and additional examples.
This is particularly useful for “what-if” sensitivity analyses:
- What is the power if the true hazard ratio is different from the design assumption?
- What happens if enrollment is slower than planned?
- What is the power at an alternate \(\alpha\) level (e.g., after \(\alpha\) reallocation)?
5.4.1 Role in the design workflow
It is useful to distinguish three questions:
- Use
gsSurv()orgsSurvCalendar()to derive a design that attains a target power. - Use
gsSurvPower()to evaluate expected operating characteristics under explicit alternative assumptions and operational timing rules. - Use simulation, for example with the simtrial package, when trial-to-trial variation in enrollment, failures, dropout, or analysis dates matters.
Thus, gsSurvPower() provides expected analysis times and event counts under one set of assumptions. It does not provide the distribution of realized analysis dates or quantify competition among operational rules in repeated trials. A disciplined scenario analysis starts from a reference design in x, overrides only assumptions that define the scenario, and uses simulation when expected values are not sufficient for the decision. See Appendix A for the class inheritance and methods available for the returned object.
5.4.2 Basic usage
An existing gsSurv design can be passed via the x parameter, with selective overrides:
# Start with a standard design
design <- gsSurv(
k = 3, test.type = 4,
alpha = 0.025, beta = 0.1,
lambdaC = log(2) / 12,
hr = 0.7, hr0 = 1,
eta = 0.01,
gamma = 10, R = 12, T = 36, minfup = 24
)
# Power under the design alternative (should be ~90%)
p1 <- gsSurvPower(x = design)
cat("Power at design HR:", p1$power, "\n")
#> Power at design HR: 0.9
# What-if: power under HR = 0.8
p2 <- gsSurvPower(x = design, hr = 0.8)
cat("Power at HR = 0.8:", p2$power, "\n")
#> Power at HR = 0.8: 0.5380404
# What-if: power under HR = 0.6
p3 <- gsSurvPower(x = design, hr = 0.6)
cat("Power at HR = 0.6:", p3$power, "\n")
#> Power at HR = 0.6: 0.9945628The scenario hazard ratio hr determines the effect under which power is computed. For designs with beta-spending futility bounds, hr1 retains the design alternative used to calibrate those bounds and defaults to the original design hazard ratio when x is supplied. Similarly, design beta is the budget used to construct futility bounds, whereas achieved beta is \(1-\) power under the scenario assumptions. These quantities generally differ after assumptions or timing rules change.
5.4.3 Analysis timing criteria
gsSurvPower() supports flexible analysis timing via several parameters that can be combined:
-
plannedCalendarTime: An exact analysis time when used alone and a calendar floor when combined with other criteria. -
targetEventsandtargetEventsPerStratum: Overall and stratum-specific event requirements. -
minTimeFromPreviousAnalysis: A minimum gap between analyses. -
minNorminNPerStratumwithminFollowUp: Enrollment and subsequent follow-up requirements. -
minfup: An explicit final-analysis floor after enrollment ends. -
maxExtension: A relative cap on delay beyond a timing anchor. -
maxCalendarTime: An absolute calendar-time cap.
These arguments form three groups. Floor criteria define the earliest eligible analysis time. Overall and active stratum-specific event targets are gates that must all be met. Caps stop the wait even if another criterion has not been met; if both relative and absolute caps are supplied, the earlier deadline applies. This provides a direct way to translate protocol language such as “not before,” “after at least,” and “but no later than” into an operating plan.
The choice between plannedCalendarTime and targetEvents has an important consequence for sensitivity analyses:
-
plannedCalendarTimefixes calendar times; expected events are recomputed under the assumed HR. A worse HR produces more events at the same calendar time. -
targetEventsfixes event counts; calendar times adjust. Results match thegsDesignpower plot to numerical precision.
5.4.4 Spending time
Analysis timing and spending time answer different questions. The spending argument selects the clock used to derive bounds:
-
"information"spends according to scenario information fractions; -
"calendar"spends according to calendar-time fractions; -
"min_planned_actual"prevents spending from running ahead of either the reference design or scenario event accumulation.
The selected effective spending-time vector is used for both upper and lower bounds. When incomplete final spending is not intended, set fullSpendingAtFinal = TRUE. Explicit informationRates can instead cap spending by planned information fractions.
5.4.5 Alternate alpha
When only alpha is changed and the effective spending schedule is unchanged, the original futility bounds are preserved while efficacy bounds are recomputed. A lower bound is clipped if it exceeds the new efficacy bound.
# Power at reduced alpha (e.g., after multiplicity adjustment)
p4 <- gsSurvPower(x = design, alpha = 0.01)
cat("Power at alpha = 0.01:", p4$power, "\n")
#> Power at alpha = 0.01: 0.83348845.5 Poisson mixture cure model
The Poisson mixture model (Rodrigues et al. 2009) is a cure model that can be useful when the failure rate in a population is expected to decline substantially over time. It also has the property that if the control group time-to-event follows a Poisson mixture distribution, then a proportional hazards assumption for treatment effect yields another Poisson mixture distribution for the experimental group. The model is flexible and easy to use: the control distribution is specified with two parameters — the cure rate and one other survival rate at an arbitrarily specified time point.
The survival function for a control group is:
\[S_c(t)=\exp(-\theta(1-\exp(-\lambda t)))\]
where \(\theta = -\log(p)\), \(p\) is the cure rate, \(\lambda > 0\) is a constant hazard rate, and \(t \ge 0\). Given the cure rate \(p = S_c(\infty) = e^{-\theta}\) and a survival rate \(S_c(t_1)\) at time \(t_1\), we can solve:
\[\theta = -\log(S_c(\infty)), \quad \lambda = -\log(1 + \log(S_c(t_1)) / \theta) / t_1\]
Under a proportional hazards assumption with hazard ratio \(\gamma > 0\), the experimental group survival is:
\[S_e(t)=\exp(-\theta\gamma(1-\exp(-\lambda t)))\]
5.5.1 Supporting functions
We define pPM() for the Poisson mixture survival function and hPM() for hazard rates:
pPM <- function(x = 0:20, cure_rate = .5, t1 = 10, s1 = .6) {
theta <- -log(cure_rate)
lambda <- -log(1 + log(s1) / theta) / t1
return(exp(-theta * (1 - exp(-lambda * x))))
}
hPM <- function(x = 0:20, cure_rate = .5, t1 = 10, s1 = .6) {
theta <- -log(cure_rate)
lambda <- -log(1 + log(s1) / theta) / t1
return(theta * lambda * exp(-lambda * x))
}5.5.2 Piecewise exponential approximation
To use the Poisson mixture model with gsSurvCalendar(), we approximate the declining hazard rates with a piecewise exponential model. We consider a scenario with a 50% cure rate, 65% survival at 24 months, hazard ratio of 0.7, and enrollment over 12 months with a total study duration of 48 months:
library(dplyr)
library(tibble)
library(ggplot2)
cure_rate <- .5
t1_val <- 24
s1_val <- .65
hr_val <- .7
study_dur <- 48
enroll_dur <- 12
bins <- 5
# Compute piecewise exponential rates from the Poisson mixture model
t_breaks <- seq(0, study_dur + 12, (study_dur + 12) / bins)
surv_control <- pPM(
x = t_breaks, cure_rate = cure_rate, t1 = t1_val, s1 = s1_val
)
cum_haz <- -log(surv_control)
lambdaC <- diff(cum_haz) / diff(t_breaks)
S <- diff(t_breaks)[1:(bins - 1)]The hazard rate declines over time, reflecting the cure model behavior:
hazard_df <- tibble(
Time = t_breaks[-1],
`Piecewise rate` = lambdaC
)
continuous_df <- tibble(
Time = seq(0, study_dur, 0.5),
Rate = hPM(seq(0, study_dur, 0.5),
cure_rate = cure_rate, t1 = t1_val, s1 = s1_val
)
)
ggplot() +
geom_step(
data = hazard_df,
aes(x = Time - diff(t_breaks)[1], y = `Piecewise rate`)
) +
geom_line(data = continuous_df, aes(x = Time, y = Rate), lty = 2) +
ylab("Hazard rate") +
xlab("Time (months)") +
ggtitle(
"Piecewise Approximation of Poisson Mixture Hazard Rate",
subtitle = "Dashed line shows actual hazard rate"
)
5.5.3 Calendar-based design with calendar spending
Due to substantial variability in event accrual caused by possible differences in true event rates, a calendar-based design with calendar-based spending (Lan and DeMets 1989) can make the timing of analyses more predictable. The gsSurvCalendar() function with spending = "calendar" replaces information fractions with calendar time fractions for spending function calculations. This tends to produce more conservative interim bounds since calendar fractions are typically lower than information fractions at early analyses.
design_pm <- gsSurvCalendar(
calendarTime = c(14, 24, 36, 48),
spending = "calendar",
alpha = 0.025,
beta = 0.1,
test.type = 4,
hr = hr_val,
R = enroll_dur,
gamma = 1,
minfup = study_dur - enroll_dur,
sfu = sfHSD,
sfupar = -3,
sfl = sfLDPocock,
lambdaC = lambdaC,
S = S
)
design_pm |>
gsBoundSummary(exclude = c("B-value", "CP", "CP H1", "PP")) |>
gt() |>
tab_header(
title = "Calendar-Based Design with Poisson Mixture Model",
subtitle = "Calendar spending, non-binding futility"
)| Calendar-Based Design with Poisson Mixture Model | |||
| Calendar spending, non-binding futility | |||
| Analysis | Value | Efficacy | Futility |
|---|---|---|---|
| IA 1: 36% | Z | 2.9057 | 0.3778 |
| N: 1038 | p (1-sided) | 0.0018 | 0.3528 |
| Events: 143 | ~HR at bound | 0.6141 | 0.9386 |
| Month: 14 | Spending | 0.0018 | 0.0406 |
| P(Cross) if HR=1 | 0.0018 | 0.6472 | |
| P(Cross) if HR=0.7 | 0.2164 | 0.0406 | |
| IA 2: 67% | Z | 2.7193 | 1.1040 |
| N: 1038 | p (1-sided) | 0.0033 | 0.1348 |
| Events: 265 | ~HR at bound | 0.7156 | 0.8730 |
| Month: 24 | Spending | 0.0027 | 0.0214 |
| P(Cross) if HR=1 | 0.0045 | 0.8841 | |
| P(Cross) if HR=0.7 | 0.5805 | 0.0620 | |
| IA 3: 88% | Z | 2.3641 | 1.6567 |
| N: 1038 | p (1-sided) | 0.0090 | 0.0488 |
| Events: 346 | ~HR at bound | 0.7755 | 0.8368 |
| Month: 36 | Spending | 0.0066 | 0.0208 |
| P(Cross) if HR=1 | 0.0107 | 0.9596 | |
| P(Cross) if HR=0.7 | 0.8215 | 0.0828 | |
| Final | Z | 2.0027 | 2.0027 |
| N: 1038 | p (1-sided) | 0.0226 | 0.0226 |
| Events: 392 | ~HR at bound | 0.8168 | 0.8168 |
| Month: 48 | Spending | 0.0139 | 0.0172 |
| P(Cross) if HR=1 | 0.0203 | 0.9797 | |
| P(Cross) if HR=0.7 | 0.9000 | 0.1000 | |
Considerations for this type of design:
- The futility bounds are advisory only. Late futility bounds may be ignored since extended follow-up may merit continuing the trial.
- The efficacy spending function should ensure evidence required to cross any bound is likely to justify early stopping.
- Substantial deviations in event accumulation do not change timing of analyses from their calendar times.
- Inference at subsequent analyses after crossing an efficacy bound can use repeated \(p\)-values (Jennison and Turnbull 2000) or sequential \(p\)-values.