• 3:30: Introduction and background theory (30 minutes)

• 4:00: Proportional hazards applications with Shiny app (25 minutes)

• 4:25: Intro to non-proportional hazards (NPH; 5 minutes)

• 4:30: Software and piecewise model (15 min)

• 4:45: Average hazard ratio (AHR; 20 minutes)

• 5:05: Break (10 minutes)

• 5:15: NPH design with logrank test (25 minutes)

• 5:40: Weighted logrank and combination tests (40 minutes)

• 6:20 Summary and questions (10 minutes)

• All opinions expressed are those of the presenters and not Merck Sharp & Dohme Corp., a subsidiary of Merck & Co., Inc., Kenilworth, NJ, USA.

• Some slides need to be scrolled down to see the full content.

For minimum effort and maximum benefit, you will probably at least want to use the Shiny interface for the gsDesign R package at https://rinpharma.shinyapps.io/gsdesign/. This is also available at https://gsdesign.shinyapps.io/prod/, but this site is only licensed for a small number of simultaneous users.

• Book at https://keaven.github.io/gsd-deming/
  • Instructions there to install software, download repository at https://github.com/keaven/gsd-deming
  • Directories there we will use:
    • data/: contains design files for examples; also simulation results
    • vignettes/: reports produced by Shiny app to summarize designs
    • simulation/: R code and simulation data for the last part of course

• Analyze a trial repeatedly at planned intervals
• Group of data added at each analysis
• Group sequential design derives boundaries and sample size to
  • Control Type I error
  • Ensure power
  • Stop early for futility or efficacy finding
• Takes advantage of correlated, group sequential tests

• Asymptotic normal assumption works well for most trials
• Scharfstein et al (1997) demonstrated \(Z = (Z_1, \dots, Z_K)\) is asymptotically normal with independent increments.
• We extend the canonical distribution notation of Jennison and Turnbull (2000):
  • \(Z_k\) is test statistics for treatment effect at analysis \(k=1,\ldots,K\)
  • \((Z_1,\ldots,Z_K)\) is multivariate normal
  • \(E(Z_k) = \theta_k \sqrt{I_k}, k=1,\ldots,K\)
  • \(\hbox{Cov}(Z_{k_1}, Z_{k_2})=\sqrt{I_{k_1}/I_{k_2}},\) \(1\le k_1\le k_2\le K\)
• \(I_k\) is the Fisher information for \(\theta_k, k=1,\ldots, K\).
• For most of this training \(\theta_k=-E(\log(HR_k))\)
• Simulation can be used to examine accuracy of normal approximation

• Lachin and Foulkes (1986) piecewise model used to approximate arbitrary enrollment, survival, dropout patterns
  • Fixed enrollment and follow-up period
    • increase enrollment rates to obtain power
  • Proportional hazards: \(\theta_k=\theta\) (constant)
  • \(I_k\) ~proportional to event count Schoenfeld (1981)
• We generalize to average hazard ratio (AHR; Mukhopadhyay et al (2020)) for non-proportional hazards (NPH) tested with logrank (\(\theta_k\) varies with \(k\))
• Generalized to cover weighted logrank as in Yung and Liu (2020) in final section today
  • Fixed design only for RMST
  • Combination tests with multiple weighted logrank tests have more complex correlation structure (Karrison and others (2016))

• Bounds \(-\infty \le a_k \le b_k \le \infty\) for \(k=1,\dots,K\)
• Null hypothesis \(H_0:\) \(\theta_k=0, k=1,\ldots,K\)
• Alternate hypothesis \(H_1:\) \(\theta_k > 0\) for some \(k=1,\ldots,K\)
• Actions at analysis \(k=1,\ldots,K\):
  • Reject \(H_0\) at analysis \(k\) if \(Z_k\ge b_k\)
  • Do not reject \(H_0\) and consider stopping if \(Z_k<a_k\), \(k<K\)
  • Continue trial if \(a_k\le Z_k\le b_k\), \(k<K\)
• Bounds are generally considered advisory for stopping trial, not binding

• Upper boundary crossing probabilities
  • \(u_k(\theta) = \text{Pr}_\theta(\{Z_k \ge b_k\} \cap_{j=1}^{k-1} \{a_j \le Z_j < b_j\})\)
• Lower boundary crossing probabilities
  • \(l_k(\theta) = \text{Pr}_\theta (\{Z_k < a_k\} \cap_{j=1}^{k-1} \{a_j \le Z_j < b_j\})\)
• Null hypothesis: 1-sided Type I error
  • \(a_k = -\infty\) for all \(k\) generally used for Type I error
    • Non-binding lower bound
  • \(\alpha = \sum_{k=1}^{K} u_k(0) = \sum_{k=1}^{K} \text{Pr}(\{Z_k \ge b_k\} \cap_{j=1}^{k-1} \{a_j \le Z_j \le b_j\} \mid H_0)\)

• Alternate hypothesis: Type II error \(\beta= 1 - \hbox{power}\)

  \[-\infty\le a_k<b_k, k=1,\ldots,K-1,\] \[a_K\le b_K\] If \(a_K = b_K\) then total Type II error is \[\beta = \sum_{k=1}^{K} l_k = \sum_{k=1}^{K} \text{Pr}(\{Z_k < a_k\} \cap_{j=1}^{k-1} \{a_j \le Z_j \le b_j\}\mid H_1)\]

Test each treatment for superiority vs the other

Usually not of interest in pharmaceutical industry

Give up if experimental arm not trending in favor of control?

• Ethics of continuing trial if unlikely to show superiority?
• Excess risk of crossing lower bound too soon?

Give up if experimental arm trending worse than control

• Given boundary computation, sample size solves for enrollment rate
• Fixing relative enrollment rates, dropout rates and trial duration enable Lachin and Foulkes (1986) and our AHR methods
• Under proportional hazards, you can also fix max enrollment rate and solve for trial duration (Kim and Tsiatis 1990)
  • You can easily create scenarios here for which there is no solution
  • Error message An error has occurred. Check your logs or contact the app author for clarification.
  • Need to adjust parameters until a solution is found or use Lachin and Foulkes (1986)

• Approaches to calculate decision boundary:

  • The error spending approach: specify boundary crossing probabilities at each analysis. This is most commonly done with the error spending function approach (Lan and DeMets 1983).

  • The boundary family approach: specify how big boundary values should be relative to each other and adjust these relative values by a constant multiple to control overall error rates. The commonly applied boundary family include:

    • Haybittle-Peto boundary (Haybittle (1971), Peto et al (1977))
    • Wang-Tsiatis boundary (Wang and Tsiatis 1987)
      • Pocock boundary (Pocock 1977)
      • O’Brien and Fleming boundary (O’Brien and Fleming 1979)
    • Slud and Wei (1982) and Fleming et al (1984) set fixed IA boundary crossing probabilities
      • Hybrid of boundary family/spending approach

• Main idea:

  • Interim Z-score boundary: 3
  • Final Z-score boundary: 1.96 (slight \(\alpha\) inflation)

• Modified Haybittle-Peto procedure 1:

  Bonferroni adjustment:

  • For the first \(K-1\) analyses, significant p-value is set as 0.001;
  • For the final analysis, significant p-value is set as \(0.05 - 0.01 \times (K-1)\);
    • More generous final bound if you adjust for test correlations

• Advantages:

  • Avoid type I error inflation
  • Does not require equally spaced analyses

• Definition:

  For 2-sided testing, Wang and Tsiatis (1987) defined the boundary function for the \(k\)-th look as \[ \Gamma(\alpha, K, \Delta) k^{\Delta - 0.5}, \] where \(\Gamma(\alpha, K, \Delta)\) is a constant chosen so that the level of significance is equal to \(\alpha\).

• Two special cases:

  • \(\Delta = 0.5\): Pocock bounds
  • \(\Delta = 0\): O’Brien-Fleming bounds.

For 2-sided testing, the Pocock procedure rejects at the \(k\)-th equally-spaced of \(K\) looks if \[|Z_k| > c_P(K),\] where \(c_P(K)\) is fixed given \(K\) such that \(\text{Pr}(\cup_{k=1}^{K} |Z_k| > c_P(K)) = \alpha\).

• Early stage: very conservative \(\Rightarrow\) large boundary at the begining;
• Final stage: nominal value close to the overall value of the design \(\Rightarrow \approx 1.96\).
• Regulators generally like these bounds

• Information fraction, Lan and DeMets (1983)
  • Time-to-event: fraction of planned final events in analysis
  • Normal or binomial: fraction of planned final sample size in analysis
  • Usually expected by regulators
• Calendar fraction, Lan and DeMets (1989)
  • Fraction of trial planned calendar duration at analysis
• Minimum of planned and actual information fraction
  • Probably not advised yet

• Sample size and power derivation
• Time-to-event endpoint
• 2-arm trial
• Logrank (or Cox model coefficient) to test treatment effect
• Constant treatment effect over time (proportional hazards or PH)
• Number of events drives power, regardless of study duration

• Shiny app available at:
  • https://rinpharma.shinyapps.io/gsdesign/ (preferred)
  • https://gsdesign.shinyapps.io/prod/ (not preferred)
  • Video intro
• Allows broad specification of designs
• Much simpler than writing code
• It will write code and reports for you
• Saves and re-loads designs
• Some flexibility of gsDesign not in app to simplify interface

• KEYNOTE 189 trial (Gandhi et al (2018))
  • Endpoints: progression free survival (PFS) and overall survival (OS) in patients
  • Indication: previously untreated metastatic non-small cell lung cancer (NSCLC)
  • Treatments: chemotherapy +/- pembrolizumab
  • Randomized 2:1 to an add-on of pembrolizumab or placebo
  • Type I error (1-sided): 0.025 familywise error rate (FWER) split between PFS (\(\alpha=0.0095\)) and OS (\(\alpha=0.0155\))
  • Graphical method \(\alpha\)-control for group sequential design of Maurer and Bretz (2013) used

Key aspects of the design as documented in the protocol accompanying Gandhi et al (2018).

• \(\alpha=0.0155\)
• Control group survival: exponential median=13 months
• Exponential dropout rate of 0.133% per month
• 90% power to detect a hazard ratio (HR) of 0.70025
• 2:1 randomization, experimental:control
• Enrollment over 1 year
• While not specified in the protocol, we have further assumed:
  • Trial duration: 35 months.
  • Observed deaths of 240, 332 and 416 at the 3 planned analyses
  • A one-sided bound using the Lan and DeMets (1983) spending approximating an O’Brien-Fleming bound.

• AFCAPS/TEXCAPS trial: use of lovastatin to reduce cardiovascular outcomes
• Design described in Downs et al (1997)
• Results reported in Downs et al (1998)
• Reproduction here not exact mainly due to choice of Lachin and Foulkes (1986)
  • Little difference between methods
  • Efficacy bounds will be exactly as proposed

• 5 years minimum follow-up of all patients enrolled
• Interim analyses after 0.375 and 0.75 of final planned event count has accrued
• 2-sided bound using the Hwang et al (1990) spending function with parameter \(\gamma = -4\) to approximate an O’Brien-Fleming bound
• We arbitrarily set the following parameters to match design:
  • Power of 90% for a hazard ratio of 0.6921846; this is slightly different than the 0.7 hazard ratio suggested in Downs et al (1997)
  • Enrollment duration of 1/2 year with constant enrollment.
  • An exponential failure rate of 0.01131 per year which is nearly identical to the annual failure rate of 0.01125.
  • An exponential dropout rate of 0.004 per year which is nearly identical to the annual dropout rate of 0.00399.

• Indication: treatment of diabetes
• Treatments: DPP4 inhibitor alogliptin compared to placebo
• Primary endpoint: major cardiovascular outcomes (MACE)
• Objective: establish non-inferiority
• Results in White et al (2013)
• Design in White et al (2011)
• We approximate the design and primary analysis evaluation here.
• Software and design assumptions not completely clear; not exact design reproduction

• Primary analysis: stratified Cox model for MACE
  • 1-sided repeated confidence interval for HR at each analysis.
  • Analysis to rule out HR > 1.3, but also tests superiority.
  • Analyses planned after 550, 600, 650 MACE events.
  • O’Brien-Fleming-like spending function Lan and DeMets (1983).
  • 2.5% Type I error.
  • Approximately 91% power.
  • 3.5% annual MACE event rate.
  • Uniform enrollment over 2 years.
  • 4.75 years trial duration.
  • 1% annual loss-to-follow-up rate.
  • Software: EAST 5 (Cytel).

Poisson mixture cure model we consider:

\[S(t)= \exp(-\theta (1 - \exp(-\lambda t)).\]

Note that:

• \(1-\exp(-\lambda t)\) is the CDF for an exponential distribution
  • can be replaced by arbitrary continuous CDF.
• As \(t\rightarrow \infty\), \(S(t)\rightarrow\exp(-\theta)\) (cure rate).
• Model useful when historical data suggests plateau in survival.
• PH model: experimental survival \(S_E(t)=S_C(t)^{HR}\).

More details in book.

• Event accumulation over time can be very sensitive to many trial design assumptions.
• Generally, we are trying to mimic a slowing of event accumulation over time.
• Assume 18 month enrollment with 6-month ramp-up.

• Quite possible that event rate design assumptions incorrect
• Good to ensure duration of follow-up is adequate evaluate tail behavior/plateau in survival
• Ensures adequate follow-up to estimate relevant parts of survival curve
• Limits trial to relevant clinical and practical duration
• May be underpowered if events accrue slowly
• Probably less useful for high-risk endpoints (e.g., metastatic cancer)
• Regulatory resistance?

See the following link for the Moderna COVID-19 design replication: https://medium.com/@yipeng_39244/reverse-engineering-the-statistical-analyses-in-the-moderna-protocol-2c9fd7544326

Can you reproduce this using the Shiny interface?

Cholesterol lowering and mortality Miettinen et al (1997)

Scandinavian Simvistatin Survival Study

• What is the impact of (potentially) delayed treatment effect on trial design and analysis?
• Test statistics used (logrank for this section; next section suggests alternatives)
• Sample size vs duration of follow-up
• Timing of analyses
• Futility bounds
• Updating bounds
• Multiplicity adjustments

simtrial, gsDesign2 and gsdmvn are under development and hosted on GitHub/Merck.

Below is the GitHub repo link to the source code

• simtrial: https://github.com/Merck/simtrial
• gsDesign2: https://github.com/Merck/gsDesign2
• gsdmvn: https://github.com/Merck/gsdmvn

How to report issues?

• Create issue under the GitHub repo listed above.

• Time-to-event trial simulation
• Piecewise model
• Logrank and weighted logrank analysis
• Simulating fixed and group sequential design
  • Also potential to simulate adaptive design
• Reverse engineered datasets from NPH Working Group
• Validation near completion
  • Thanks to AP colleagues and Amin Shirazi

• Generating simulated datasets
  • simfix(): Simulation of fixed sample size design for time-to-event endpoint
  • simfix2simPWSurv(): Conversion of enrollment and failure rates from simfix() to simPWSurv() format
  • simPWSurv(): Simulate a stratified time-to-event outcome randomized trial
• Cutting data for analysis
  • cutData(): Cut a dataset for analysis at a specified date
  • cutDataAtCount(): Cut a dataset for analysis at a specified event count
  • getCutDateForCount(): Get date at which an event count is reached
• Analysis functions
  • tenFH(): Fleming-Harrington weighted logrank tests
  • tenFHcorr(): Fleming-Harrington weighted logrank tests plus correlations
  • tensurv(): Process survival data into counting process format
  • pMaxCombo(): MaxCombo p-value
  • pwexpfit(): Piecewise exponential survival estimation
  • wMB(): Magirr and Burman modestly weighted logrank tests
• Lower level functions
  • fixedBlockRand(): Permuted fixed block randomization
  • rpwenroll(): Generate piecewise exponential enrollment
  • rpwexp(): The piecewise exponential distribution

From NPH Working Group

• Ex1delayedEffect: Time-to-event data example 1 for non-proportional hazards working group
• Ex2delayedEffect: Time-to-event data example 2 for non-proportional hazards working group
• Ex3curewithph: Time-to-event data example 3 for non-proportional hazards working group
• Ex4belly: Time-to-event data example 4 for non-proportional hazards working group
• Ex5widening: Time-to-event data example 5 for non-proportional hazards working group
• Ex6crossing: Time-to-event data example 6 for non-proportional hazards working group
• MBdelayed: Simulated survival dataset with delayed treatment effect

• AHR(): Average hazard ratio under non-proportional hazards (test version)
• eAccrual(): Piecewise constant expected accrual
• eEvents_df(): Expected events observed under piecewise exponential model
• ppwe(): Estimate piecewise exponential cumulative distribution function
• s2pwe(): Approximate survival distribution with piecewise exponential distribution
• tEvents(): Predict time at which a targeted event count is achieved

We will focus on AHR(), ppwe(), and s2pwe() in this training.

Power and design functions extending the Jennison and Turnbull (2000) computational model to non-constant treatment effects. Partial list of functions:

• Design and power under average hazard ratio model
  • gs_power_ahr(): Power computation
  • gs_design_ahr(): Design computations
• Bound support
  • gs_b(): direct input of bounds
  • gs_spending_bound(): spending function bounds
• Other tests
  • Design and power with weighted logrank and MaxCombo will be discussed in next section.

• Simple model to approximate arbitrary patterns of
  • Enrollment: piecewise constant enrollment rates
  • Failure rates: piecewise exponential
  • Dropout rates: piecewise exponential
• Combined tools for designing and evaluating designs
  • Asymptotic approach using average hazard ratio (AHR)
  • Simulation tools to confirm asymptotic approximations
  • No requirement for proportional hazards
  • Stick with logrank for this section

• Should we assume all enrollment
  • by targeted time?
  • random, expected enrollment matches target at targeted time?
• simtrial package assumes the latter
  • Also used in AHR approach here
• npsurvSS package assumes the former
  • Used for weighted logrank and combination tests in next section

• Definition of average hazard ratio
• Literature review
• Asymptotic theory for group sequential design
• Verification by simulation

• Geometric mean hazard ratio (Mukhopadhyay et al (2020))
• Exponentiate: average \(\log(\hbox{HR})\) weighted by expected events per interval

\[\hbox{AHR} = \exp\left( \frac{d_1 \log(1) + d_2 \log(0.6)}{d_1 + d_2}\right)\]

• Assume \((t_{m-1},t_m]\) that the hazard rate \(\lambda_m\) is constant
• Assume total time (all observations) in interval is \(T_m\) (total time on test)
• Assume event count in interval is \(D_m\)
• Hazard rate estimate in this interval: \[\hat{\lambda}_m= D_m/T_m\]
• Asymptotic approximation (delta method; see Lu Tian slides) \[\log(\hat\lambda)\sim \hbox{Normal}(\log(\lambda_m),\sigma^2=1/D_m).\]

• Now consider hazard ratio for two treatment arms, \(j=0,1\) (0=control, 1=experimental)
• Same exponential model
• Log hazard ratio in interval \[\beta_m=\log(\lambda_{1m}/\lambda_{0m})=\log(\lambda_{1m}) - \log(\lambda_{0m})\]
• Asymptotic approximation \[\hat{\beta}_m\sim \hbox{Normal}(\beta_m,1/D_{0m}+ 1/D_{1m})\]
• Now we replace \(D_{0m}, D_{1m}\) with their expected values \(E(D_{0m}), E(D_{1m})\)

• Schoenfeld (1981) approximation works well plugging in average hazard ratio
• Use gsDesign::nEvents()

Events <- 332
# if beta is NULL and n= # of events,
# power is computed instead of events required
Power <- gsDesign::nEvents(n = Events, beta = NULL, hr = c(.6, .7, .8))

Assume 332 events

• Some argue this is a bad idea
  • e.g., hazards of hazard ratios (Hernán (2010))
• Pro’s
  • Estimated by Cox regression
  • AHR concept makes more clear what this is
  • Logrank is widely-accepted corresponding test
  • Both asymptotic approximations and simulation supported
    • This includes group sequential design
  • Easy to approximate arbitrary enrollment, failure and dropout patterns
• Cautions
  • No single estimand sufficently describes NPH differences
  • Early interim analysis (futility, efficacy) should anticipate possible reduced effect

Used the following code to get AHR and information at specified times:

analysisTimes <- c(12, 20, 28, 36)
sampleSize <- 500
enrollRates <- tibble(Stratum = "All", duration = 12, rate = sampleSize / 12)
failRates <- tibble(
  Stratum = "All",
  duration = c(4, 100),
  failRate = log(2) / 15,
  hr = c(1, .6),
  dropoutRate = 0.001
)
ahr <- gsDesign2::AHR(
  enrollRates = enrollRates,
  failRates = failRates,
  totalDuration = analysisTimes,
  ratio = 1
)

• Many options for cutting data!
  • Time-based
  • Event-based
  • Max of time- and event-based
  • Limit final trial duration
• Note that
  • Code here is detailed, but (hopefully) not complicated
  • AHR changes if event-based
  • Number of events variable if time-based

# Transform to simPWSurv() format
x <- simfix2simPWSurv(failRates)
nsim <- 10000 # of simulations
# Set up matrix for simulation results
results <- matrix(0, nrow = nsim * 4, ncol = 6)
colnames(results) <- c("Sim", "Analysis", "Events", "beta", "var", "logrank")
ii <- 0 # index for results row
for (sim in 1:nsim) {
  # Simulate a trial
  ds <- simPWSurv(
    n = sampleSize,
    enrollRates = enrollRates,
    failRates = x$failRates,
    dropoutRates = x$dropoutRates
  )
  for (j in seq_along(analysisTimes)) {
    # Cut data at specified analysis times
    # Use cutDataAtCount to cut at event count
    dsc <- ds %>% cutData(analysisTimes[j])
    ii <- ii + 1
    results[ii, 1] <- sim
    results[ii, 2] <- j
    results[ii, 3] <- sum(dsc$event)
    cox <- coxph(Surv(tte, event) ~ Treatment, data = dsc)
    results[ii, 4] <- as.numeric(cox$coefficients)
    results[ii, 5] <- as.numeric(cox$var)
    # Logrank test
    Z <- dsc %>%
      tensurv(txval = "Experimental") %>%
      tenFH(rg = tibble::tibble(rho = 0, gamma = 0))
    results[ii, 6] <- as.numeric(Z$Z)
  }
}

Simulation summary based on 10k replications

results <- tibble::as_tibble(results)
simsum <- results %>%
  group_by(Analysis) %>%
  summarize(
    AHR = exp(mean(beta)), Events = mean(Events),
    info = 1 / mean(var(beta)), info0 = Events / 4
  )

• Simulations represent truth since they are from actual distribution
• Note that asymptotic info0 seems a better approximation of simulation than info
• Using both info0 and info in design will make sample size a little more conservative
  • Not as conservative as simulation info

ggplot(results, aes(x = factor(Analysis), y = beta)) +
  geom_violin() +
  ggtitle("Distribution of Cox Coefficient by Analysis") +
  xlab("Analysis") +
  ylab("Cox coefficient")

9 simulations

Question: Do you really want to adapt sample size based on an early interim estimate of treatment effect?

• Use of Tsiatis (1982) (also extends to weighted logrank)
• Statistical information proportional to expected event counts as in Schoenfeld (1981)
• Natural parameter: \(\log(\hbox{AHR})\)
• Statistical information still proportional to number of events
• Correlation still computed based on statistical information
• Extension of Jennison and Turnbull (2000) calculations to non-constant effect size over time

Statistical information at analysis: \(\mathcal{I}_k\), \(1\le k\le K\)

Proportion of final information at analysis \(k\): \(t_k =\mathcal{I}_k / \mathcal{I}_K\)

\[Z_k\sim \hbox{Normal}(\sqrt{\mathcal{I}_k} \theta(t_k),1)\] Multivariate normal with correlations for \(1\le j\le k\le K\):

\[\hbox{Corr}(Z_j,Z_k)=\sqrt{t_j/t_k}\]

• Previously, Cox \(\hat\beta\) was shown to be approximately a weighted \(\log(AHR)\)
• How do we translate to extended canonical form with \(\theta(t_k)\), \(\mathcal{I_k}\), \(1\le k\le K\)?
  • Assume piecewise model
  • Fix some \(k\) in \(1,2,\ldots,K\)
    • Fix calendar time of analysis \(k\) relative to opening of enrollment; can solve for this if event-based spending fraction desired
    • Compute \(\mathcal{I_k}\) as before based on expected events
    • Compute \(t_k=\mathcal{I}_k/\mathcal{I}_K\)
    • Compute expected Cox coefficient \(\beta\) (\(=\log(AHR)\)) as before and set \(\theta(t_k)=\beta\)
• We now have canonical form for the piecewise model

analysisTimes <- c(12, 20, 28, 36)
sampleSize <- 500

enrollRates

failRates

ahr

Information fraction for interim analyses

## [1] 0.3241690 0.6275343 0.8424726

PH1sided <- gsDesign::gsSurv( # Derive Group Sequential Design
  k = 4, # Number of analyses (interim + final)
  test.type = 1, # use this for 1-sided testing
  alpha = 0.025, # 1-sided Type I error
  beta = 0.1, # Type II error (1 - power)
  timing = timing, # Information fraction for interims
  sfu = sfLDOF, # O'Brien-Fleming spending approximation
  lambdaC = failRates$failRate, # Piecewise control failure rates
  hr = ahr$AHR[4], # Used final analysis AHR
  eta = failRates$dropoutRate, # Piecewise exponential dropout rates
  gamma = enrollRates$rate, # Relative enrollment
  R = enrollRates$duration, # Duration of piecewise enrollment rates
  S = failRates$duration[1], # Duration of piecewise failure rates (K-1)
  T = max(analysisTimes), # Study duration
  minfup = max(analysisTimes) - sum(enrollRates$duration), # Minimum follow-up
  ratio = 1 # Experimental:Control randomization ratio
)

cat(summary(PH1sided))

One-sided group sequential design with 4 analyses, time-to-event outcome with sample size 444 and 297 events required, 90 percent power, 2.5 percent (1-sided) Type I error to detect a hazard ratio of 0.68. Enrollment and total study durations are assumed to be 12 and 36 months, respectively. Efficacy bounds derived using a Lan-DeMets O’Brien-Fleming approximation spending function with none = 1.

library(gsdmvn)
# Spending function setup
upar <- list(sf = gsDesign::sfLDOF, total_spend = 0.025)
NPH1sided <- gs_design_ahr(
  enrollRates = enrollRates,
  failRates = failRates,
  ratio = 1, alpha = .025, beta = 0.1,
  # Information fraction not required (but available!)
  analysisTimes = analysisTimes,
  # Function to enable spending bound
  upper = gs_spending_bound,
  # Spending function and parameters used
  upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025),
  # Lower bound fixed at -infinity
  lower = gs_b, # allows input of fixed bound
  # With gs_b, just enter values for bounds
  lpar = rep(-Inf, 4)
)

• You will want to round up events and sample size!
• Interim boundary crossing probability much lower than with PH bounds
  • IA 2 crossing probability 50% under PH
• Sample size larger than for PH (N=444, 297 events)

• Moving on to lower bounds
• Common practice is to use binding upper and lower bounds
• Here we use two one-sided tests for \(\alpha-\)spending
• Upper bound \(b_k(\alpha)\) and lower bound \(a_k(\alpha)= - b_k(\alpha)\) defined so that \[ \begin{align} f(s_k,\alpha)-f(s_{k-1},\alpha) =& P_0(\{Z_{k}\geq b_{k}(\alpha)\}\cap_{j=1}^{k-1}\{-b_{j}(\alpha)< Z_{j}< b_{j}(\alpha)\}\\ =& P_0(\{Z_{k}\le -b_{k}(\alpha)\}\cap_{j=1}^{k-1}\{-b_{j}(\alpha)< Z_{j}< b_{j}(\alpha)\} \end{align} \]

library(gsdmvn)

# Spending function and parameters for both bounds
par <- list(sf = gsDesign::sfLDOF, total_spend = 0.025)
NPHsymmetric <- gs_design_ahr(
  enrollRates = enrollRates,
  failRates = failRates,
  ratio = 1, alpha = .025, beta = 0.1,
  # Information fraction not required (but available!)
  analysisTimes = analysisTimes,
  # Function to enable spending bound
  upper = gs_spending_bound,
  lower = gs_spending_bound,
  # Spending function and parameters used
  upar = par,
  lpar = par,
  binding = TRUE, # set lower bound to binding
  h1_spending = FALSE
)

• Use non-binding upper bound with spending function \(f_1(s,\alpha)\)
• Set lower boundary crossing boundary under null hypothesis; may assume binding for lower bound with spending function \(f_2(s,\gamma)\) for some chosen \(0 < \gamma \le 1-\alpha\)
• Set bounds to satisfy \[ \begin{align} f_1(s_k,\alpha)-f_1(s_{k-1},\alpha) & = P_0(\{Z_{k}\geq b_{k}(\alpha)\}\cap_{j=1}^{k-1}\{Z_{j}< b_{j}(\alpha)\}\\ f_2(s_k,\gamma)-f_2(s_{k-1},\gamma) & = P_\theta(\{Z_{k}< a_{k}(\gamma)\}\cap_{j=1}^{k-1}\{a_{j}(\gamma)\le Z_{j}< b_{j}(\alpha)\} \end{align} \]

• Use non-binding upper bound with spending function \(f_1(s,\alpha)\)
• Set lower boundary crossing boundary under alternate hypothesis \(\theta(t_k)\) (denoted \(\theta\) below)
• Assume spending function \(f_2(s,\gamma)\) for \(\gamma\) set to Type II error; power \(=1-\gamma\)
• Set bounds to satisfy

\[ \begin{align} f_1(s_k,\alpha)-f_1(s_{k-1},\alpha) &= P_0(\{Z_{k}\geq b_{k}(\alpha)\}\cap_{j=1}^{k-1}\{Z_{j}< b_{j}(\alpha)\}\\ f_2(s_k,\gamma)-f_2(s_{k-1},\gamma) &= P_\theta(\{Z_{k}< a_{k}(\gamma)\}\cap_{j=1}^{k-1}\{a_{j}(\gamma)\le Z_{j}< b_{j}(\alpha)\} \end{align} \] - Generally, sample size set so that \(a_K=b_K\)

library(gsdmvn)
# Spending function setup
upar <- list(sf = gsDesign::sfLDOF, total_spend = 0.025)
lpar <- list(sf = gsDesign::sfHSD, total_spend = .1, param = -2)
NPHasymmetric <- gs_design_ahr(
  enrollRates = enrollRates,
  failRates = failRates,
  ratio = 1, alpha = .025, beta = 0.1,
  # Information fraction not required (but available!)
  analysisTimes = analysisTimes,
  # Function to enable spending bound
  upper = gs_spending_bound,
  lower = gs_spending_bound,
  # Spending function and parameters used
  upar = upar, lpar = lpar
)

• Not easily done with gsDesign::gsSurv()
• Futility only at interim 1
  • Look for p=0.05 in the wrong direction
• Efficacy only AFTER interim 1
• This is a variation on asymmetric design
• Cannot be done (at least not easily) with gsDesign package
• We will use information fraction instead of calendar times of analysis

# Spending function setup
upar <- list(sf = gsDesign::sfLDOF, total_spend = 0.025)
lpar <- c(qnorm(.05), rep(-Inf, 3))
NPHskip <- gs_design_ahr(
  enrollRates = enrollRates,
  failRates = failRates,
  ratio = 1, alpha = .025, beta = 0.1,
  # Information fraction not required (but available!)
  analysisTimes = analysisTimes,
  # Upper spending bound
  upper = gs_spending_bound, upar = upar,
  # Skip first efficacy analysis
  test_upper = c(FALSE, TRUE, TRUE, TRUE),
  # Spending function and parameters used
  lower = gs_b, lpar = lpar
)

All probabilities are under null hypothesis

• Nominal p-value: probability of rejecting for a specific test conditioning on nothing else
• Repeated p-value: \(\alpha\)-level at which a group sequential test would be rejected at a given analysis for a specific hypothesis
• Sequential p-value: \(\alpha\)-level at which a group sequential test would be rejected for all analyses performed
  • Both sequential and repeated p-values can be performed on interim data
• FWER - Given a set of hypotheses with fixed and/or group sequential designs, the FWER at which an individual hypothesis would be rejected with one of the analyses performed

We consider other alternative tests for group sequential design.

• Review fixed design
  • Weighted logrank test
  • MaxCombo test
  • Illustration using npsurvSS
• Group sequential design with weighted logrank test
  • Under a given boundary
  • Boundary calculation
  • Illustration using gsdmvn
• Group sequential design with MaxCombo test
  • Under a given boundary
  • Boundary calculation
  • Illustration using gsdmvn

For simplicity, we made a few key assumptions.

• Balanced design (1:1 randomization ratio).
• 1-sided test.
• Local alternative: variance under null and alternative are approximately equal.
• Accrual distribution: Piecewise uniform.
• Survival distribution: piecewise exponential.
• Loss to follow-up: exponential.
• No stratification.
• No cure fraction.

The fixed design part largely follows the concept described in Yung and Liu (2019).

• \(\alpha\): Type I error
• \(\beta\): Type II error or power (1 - \(\beta\))
• \(z_\alpha\): upper \(\alpha\) percentile of standard normal distribution
• \(z_\beta\): upper \(\beta\) percentile of standard normal distribution

We considered a 1-sided test with type I error at \(\alpha=0.025\) and \(1-\beta=80\%\) power.

z_alpha <- abs(qnorm(0.025))
z_alpha

## [1] 1.959964

z_beta <- abs(qnorm(0.2))
z_beta

## [1] 0.8416212

• \(\theta\): effect size
• \(n\): total sample size
• \(Z\): test statistics is asymptotic normal
  • Under null hypothesis: \(Z \sim \mathcal{N}(0, \sigma_0^2)\)
  • Under alternative hypothesis: \(Z \sim \mathcal{N}(\sqrt{n}\theta, \sigma_1^2)\)

By assuming local alternative, we have

\[\sigma_0^2 \approx \sigma_1^2 = \sigma^2\] In this simplified case, the sample size can be calculated as

\[ n = \frac{4 (z_{\alpha}+z_{\beta})^{2}}{\theta^2} \]

• Two-sample t-test
• Logrank test

• https://keaven.github.io/gsd-deming/fixed-design.html#non-proportional-hazards

• Examples
• Technical details are skipped

• Examples

Similar to the fixed design, we can define the test statistics for weighted logrank test using counting process formula

\[ Z_k=\sqrt{\frac{n_{0}+n_{1}}{n_{0}n_{1}}}\int_{0}^{t_k}w(t)\frac{\overline{Y}_{0}(t)\overline{Y}_{1}(t)}{\overline{Y}_{0}(t)+\overline{Y}_{0}(t)}\left\{ \frac{d\overline{N}_{1}(t)}{\overline{Y}_{1}(t)}-\frac{d\overline{N}_{0}(t)}{\overline{Y}_{0}(t)}\right\} \]

Note, the only difference is that the test statistics fixed analysis up to \(t_k\) at \(k\)-th interim analysis

• Examples with pre-defined boundary
• Examples with calculated boundary
• Technical details are skipped

• Examples with pre-defined boundary
• Examples with calculated boundary
• Technical details are skipped

• Understand your control group
  • Get simple assumptions to approximate published results?
• Understand AHR
  • Consider treatment effect that is clinically meaningful, but conservative
  • Consider a treatment effect delay that is intermediate
  • Understand likely enrollment pattern and dropout rate
  • Plot AHR, underlying survival differences and expected accrual of patients and events

• Follow-up (trial duration)
  • Get to good part of AHR could (it will plateau under above assumptions)
  • Ensure follow-up long enough to ensure tail characterization
  • Require both minimum follow-up and sufficient events for final analysis

• Start with AHR from above evaluation
• Design under proportional hazards assuming constant AHR
• Ensure early futility bounds are not too aggressive
  • No futility (risk that data monitoring committee will overrule)
  • \(\beta\)-spending: be conservative for early bounds
  • Asymmetric 2-sided test; e.g., Pocock-like bound with moderate total boundary crossing probability
• Efficacy bounds: O’Brien-Fleming always acceptable to regulators
• This is easy to do with gsDesign Shiny app!

• For primary analysis, may wish to stick with logrank test
  • Well-accepted by regulatory agencies
• Futility bounds and power take into account changing effect over time
  • NPH will also impact expected timing of analyses
• More easily allows interim timing based on calendar time
• More options for setting boundaries
  • Fixed futility bounds
  • Haybittle-Peto bounds
  • Futility or efficacy bounds can be eliminated from selected analyses

• RMST has been heavily promoted
  • For delayed effect, power is generally no better than logrank
  • However, some advocate these for primary analysis
• Sensitivity analyses with better power assuming delayed effect
  • Weighted logrank with FH(0, 0.5) or Magirr-Burman test
  • MaxCombo test at final analysis only
• There can be substantial power gains

• Feedback and questions are welcome!
• May wish to submit by issues at GitHub repositories, but e-mail also OK.