4  Sample size and power for a group sequential design

We extend the fixed design calculations performed using nSurv() to compute sample size and power for a group sequential design. Classically, the required number of events is increased by a factor depending on the group sequential bounds to obtain an inflated statistical information. Under the null hypothesis of equal failure distributions, this is proportional to statistical information. While this adds another layer of approximation, we will follow this practice here. Thus, we will increase either accrual rates, accrual duration or follow-up duration to achieve the desired expected number of events under the alternate hypothesis. There are two required computations we will apply in deriving a group sequential design: first, inflating the expected number of events and second, deriving the times at which the expected number of events for each interim analysis equal the desired number.

4.1 Inflating the expected number of events

There are three cases to consider corresponding to the three ways we solved the sample size for a fixed design. These involve solving for 1) accrual rate, 2) accrual duration or 3) minimum follow-up duration required to obtain a desired number of events. We explain here as if R is a scalar; for a vector R the calculations are more complicated, but analogous.

  1. Solving for accrual rate (when both T and minfup are provided): The accrual duration R is fixed at T - minfup, and the accrual rate(s) gamma are increased proportionally to achieve the inflated number of events. This also increases the expected sample size by the same proportion.

  2. Solving for accrual duration (when T is NULL but minfup is provided): The accrual rate(s) gamma are kept fixed, and the accrual duration R is extended or shortended to achieve the inflated number of events. The study duration T is then computed as R + minfup.

  3. Solving for minimum follow-up duration (when minfup is NULL): Both accrual rate(s) gamma and accrual duration R are kept fixed, and the minimum follow-up duration minfup is extended or shortended to achieve the inflated number of events. The study duration T is then computed as R + minfup.

The first case is straightforward: to increase the expected number of events by an inflation factor when the study duration and relative enrollments are fixed over time requires only to increase the enrollment rate(s) by the same proportion. The other two computations are simple applications of the function KTZ().

4.2 Source code for inflating the expected number of events

We consider fixed design and a targeted number of events greater than that for a fixed design.

gsnSurv <- function(x, nEvents) {
  if (x$variable == "Accrual rate") {
    Ifct <- nEvents / x$d
    rval <- list(
      lambdaC = x$lambdaC, etaC = x$etaC, etaE = x$etaE,
      gamma = x$gamma * Ifct, ratio = x$ratio, R = x$R, S = x$S, T = x$T,
      minfup = x$minfup, hr = x$hr, hr0 = x$hr0, n = x$n * Ifct, d = nEvents,
      eDC = x$eDC * Ifct, eDE = x$eDE * Ifct, eDC0 = x$eDC0 * Ifct,
      eDE0 = x$eDE0 * Ifct, eNC = x$eNC * Ifct, eNE = x$eNE * Ifct,
      variable = x$variable
    )
  } else if (x$variable == "Accrual duration") {
    y <- KT(
      n1Target = nEvents, minfup = x$minfup, lambdaC = x$lambdaC,
      etaC = x$etaC, etaE = x$etaE,
      R = x$R, S = x$S, hr = x$hr, hr0 = x$hr0, gamma = x$gamma,
      ratio = x$ratio, tol = x$tol
    )
    rval <- list(
      lambdaC = x$lambdaC, etaC = x$etaC, etaE = x$etaE,
      gamma = x$gamma, ratio = x$ratio, R = y$R, S = x$S, T = y$T,
      minfup = y$minfup, hr = x$hr, hr0 = x$hr0,
      n = y$n, d = nEvents,
      eDC = y$eDC, eDE = y$eDE, eDC0 = y$eDC0,
      eDE0 = y$eDE0, eNC = y$eNC, eNE = y$eNE, tol = x$tol,
      variable = x$variable
    )
  } else {
    y <- KT(
      n1Target = nEvents, minfup = NULL,
      lambdaC = x$lambdaC, etaC = x$etaC,
      etaE = x$etaE, R = x$R, S = x$S, hr = x$hr, hr0 = x$hr0,
      gamma = x$gamma, ratio = x$ratio, tol = x$tol
    )
    rval <- list(
      lambdaC = x$lambdaC, etaC = x$etaC, etaE = x$etaE,
      gamma = x$gamma, ratio = x$ratio, R = x$R, S = x$S, T = y$T,
      minfup = y$minfup, hr = x$hr, hr0 = x$hr0, n = y$n,
      d = nEvents, eDC = y$eDC, eDE = y$eDE, eDC0 = y$eDC0,
      eDE0 = y$eDE0, eNC = y$eNC, eNE = y$eNE, tol = x$tol,
      variable = x$variable
    )
  }
  class(rval) <- "gsSize"
  return(rval)
}
x <- nSurv(
  lambdaC = log(2) / matrix(c(3, 4, 5, 6, 8, 10, 9, 12, 15), nrow = 3),
  gamma = matrix(c(2, 4, 8, 3, 6, 10), nrow = 2),
  hr = .6, S = c(3, 6), R = c(3, 15), minfup = NULL, T = 15
)
print(x)
nSurv fixed-design summary (method=LachinFoulkes; target=Follow-up duration)
HR=0.600 vs HR0=1.000 | alpha=0.025 (sided=1) | power=90.0%
N=303.0 subjects | D=160.9 events | T=21.5 study duration | accrual=18.0 Accrual duration | minfup=3.5 minimum follow-up | ratio=1 randomization ratio (experimental/control)

Expected events by stratum (H1):
  control experimental  total
1  27.796       22.718 50.514
2  23.340       17.331 40.671
3  41.294       28.418 69.713

Key inputs (names preserved):
                               desc    item                           value
                    Accrual rate(s)   gamma             2, 4, 8 ... (len=6)
           Accrual rate duration(s)       R                           3, 15
             Control hazard rate(s) lambdaC 0.231, 0.173, 0.139 ... (len=9)
            Control dropout rate(s)     eta             0, 0, 0 ... (len=9)
       Experimental dropout rate(s)    etaE             0, 0, 0 ... (len=9)
 Event and dropout rate duration(s)       S                            3, 6
                                                    input
                   matrix(c(2, 4, 8, 3, 6, 10), nrow = 2)
                                                 c(3, 15)
 log(2)/matrix(c(3, 4, 5, 6, 8, 10, 9, 12, 15), nrow = 3)
                                                      eta
                                                     etaE
                                                  c(3, 6)
gsnSurv(x, 170)
$lambdaC
      Stratum 1  Stratum 2  Stratum 3
0-3   0.2310491 0.11552453 0.07701635
3-9   0.1732868 0.08664340 0.05776227
9-Inf 0.1386294 0.06931472 0.04620981

$etaC
      Stratum 1 Stratum 2 Stratum 3
0-3           0         0         0
3-9           0         0         0
9-Inf         0         0         0

$etaE
      Stratum 1 Stratum 2 Stratum 3
0-3           0         0         0
3-9           0         0         0
9-Inf         0         0         0

$gamma
     Stratum 1 Stratum 2 Stratum 3
0-3          2         8         6
3-18         4         3        10

$ratio
[1] 1

$R
[1]  3 15

$S
[1] 3 6

$T
[1] 22.85288

$minfup
[1] 4.85288

$hr
[1] 0.6

$hr0
[1] 1

$n
[1] 303

$d
[1] 170

$eDC
[1] 28.77412 24.39549 44.04440

$eDE
[1] 23.89893 18.31081 30.57626

$eDC0
[1] 26.85077 21.75258 37.84991

$eDE0
[1] 26.85077 21.75258 37.84991

$eNC
[1] 33.0 34.5 84.0

$eNE
[1] 33.0 34.5 84.0

$tol
[1] 0.0001220703

$variable
[1] "Follow-up duration"

attr(,"class")
[1] "gsSize"
nSurv fixed-design summary (method=LachinFoulkes; target=Follow-up duration)
HR=0.600 vs HR0=1.000 | alpha=0.025 (sided=1) | power=90.0%
N=303.0 subjects | D=160.9 events | T=21.5 study duration | accrual=18.0 Accrual duration | minfup=3.5 minimum follow-up | ratio=1 randomization ratio (experimental/control)

Expected events by stratum (H1):
  control experimental  total
1  27.796       22.718 50.514
2  23.340       17.331 40.671
3  41.294       28.418 69.713

Key inputs (names preserved):
                               desc    item                           value
                    Accrual rate(s)   gamma             2, 4, 8 ... (len=6)
           Accrual rate duration(s)       R                           3, 15
             Control hazard rate(s) lambdaC 0.231, 0.173, 0.139 ... (len=9)
            Control dropout rate(s)     eta             0, 0, 0 ... (len=9)
       Experimental dropout rate(s)    etaE             0, 0, 0 ... (len=9)
 Event and dropout rate duration(s)       S                            3, 6
                                                    input
                   matrix(c(2, 4, 8, 3, 6, 10), nrow = 2)
                                                 c(3, 15)
 log(2)/matrix(c(3, 4, 5, 6, 8, 10, 9, 12, 15), nrow = 3)
                                                      eta
                                                     etaE
                                                  c(3, 6)

4.3 Enrollment at interim analysis

We begin with a simple function to compute enrollment given enrollment rates, enrollment periods and duration of enrollment. The input variable y should be an object of class nSurv, gsSize or gsSurv. This specifies enrollment and event rates. nEventsIA() computes the expected number of events at an interim analysis performed at a specific time specified in the variable tIA. tEventsIA() uses root-finding to derive the time at which the expected number of events is equal to a planned interim proportion (specified in the input timing) of the final planned events. The input variable tol will normally be left as the default by the user; this is used in the root-finding routine uniroot(). Note that nEventsIA() was updated in August, 2015 so that the simple = TRUE option works correctly for stratified populations; this impacted tEventsIA() as well.

# tEventsIA roxy [sinew] ----
#' @export
#' @rdname nSurv
#' @seealso
#'  \code{\link[stats]{uniroot}}
#' @importFrom stats uniroot
# tEventsIA function [sinew] ----
tEventsIA <- function(x, timing = .25,
                      tol = .Machine$double.eps^0.25) {
  T <- x$T[length(x$T)]
  z <- uniroot(
    f = nEventsIA, interval = c(0.0001, T - .0001), x = x,
    target = timing, tol = tol, simple = TRUE
  )
  return(nEventsIA(tIA = z$root, x = x, simple = FALSE))
}
# nEventsIA roxy [sinew] ----
#' @export
#' @rdname nSurv
# nEventsIA function [sinew] ----
nEventsIA <- function(tIA = 5, x = NULL, target = 0,
                      simple = TRUE) {
  Qe <- x$ratio / (1 + x$ratio)
  eDC <- eEvents(
    lambda = x$lambdaC, eta = x$etaC,
    gamma = x$gamma * (1 - Qe), R = x$R, S = x$S, T = tIA,
    Tfinal = x$T[length(x$T)], minfup = x$minfup
  )
  eDE <- eEvents(
    lambda = x$lambdaC * x$hr, eta = x$etaC,
    gamma = x$gamma * Qe, R = x$R, S = x$S, T = tIA,
    Tfinal = x$T[length(x$T)], minfup = x$minfup
  )
  if (simple) {
    if (class(x)[1] == "gsSize") {
      d <- x$d
    } else if (!is.matrix(x$eDC)) {
      d <- sum(x$eDC[length(x$eDC)] + x$eDE[length(x$eDE)])
    } else {
      d <- sum(x$eDC[nrow(x$eDC), ] + x$eDE[nrow(x$eDE), ])
    }
    return(sum(eDC$d + eDE$d) - target * d)
  } else {
    return(list(
      T = tIA, eDC = eDC$d, eDE = eDE$d,
      eNC = eDC$n, eNE = eDE$n
    ))
  }
}

4.4 The general group sequential design approach

We provide two variations, one with event-based (information-based) timing of analyses (gsSurv()), the other with calendar-based timing (gsSurvCalendar()). gsSurv() is being updated in August, 2023 as it inadvertently did not return expected event counts under null hypothesis when originally written. gsSurvCalendar() is being added in August, 2023. In addition to adding calendar-based timing, gsSurvCalendar() can be used to compute power given enrollment, failure rate, and dropout rate assumptions. This capability is not currently available for gsSurv(). In January, 2026 we are correcting inverse-variance weighting for strata when computing sample size for stratified populations. Also, in January, 2026 we are adding options to compute events required using the Schoenfeld or Freedman methods.

# gsSurv roxy [sinew] ----
#' @export
# gsSurv function [sinew] ----
gsSurv <- function(k = 3, test.type = 4, alpha = 0.025, sided = 1,
                   beta = 0.1, astar = 0, timing = 1, sfu = gsDesign::sfHSD, sfupar = -4,
                   sfl = gsDesign::sfHSD, sflpar = -2, r = 18,
                   lambdaC = log(2) / 6, hr = .6, hr0 = 1, eta = 0, etaE = NULL,
                   gamma = 1, R = 12, S = NULL, T = NULL, minfup = NULL, ratio = 1,
                   tol = .Machine$double.eps^0.25,
                   usTime = NULL, lsTime = NULL,
                   method = c("LachinFoulkes", "Schoenfeld", "Freedman", "BernsteinLagakos")) 
{
  method <- match.arg(method)
  # Validate ratio is a single positive scalar
  if (!is.numeric(ratio) || length(ratio) != 1 || ratio <= 0) {
    stop("ratio must be a single positive scalar")
  }
  # If both gamma and R are provided (non-NULL) and T is NULL, set T to force solving for accrual rate
  # This matches gsDesign::gsSurv behavior which keeps R fixed and solves for gamma
  if (is.null(T) && !is.null(minfup) && 
      !is.null(R) && length(R) > 0 && 
      !is.null(gamma) && length(gamma) > 0) {
    T <- sum(R) + minfup
  }
  x <- nSurv(
    lambdaC = lambdaC, hr = hr, hr0 = hr0, eta = eta, etaE = etaE,
    gamma = gamma, R = R, S = S, T = T, minfup = minfup, ratio = ratio,
    alpha = alpha, beta = beta, sided = sided, tol = tol, method = method
  )
  y <- gsDesign::gsDesign(
    k = k, test.type = test.type, alpha = alpha / sided,
    beta = beta, astar = astar, n.fix = x$d, timing = timing,
    sfu = sfu, sfupar = sfupar, sfl = sfl, sflpar = sflpar, tol = tol,
    delta1 = log(hr), delta0 = log(hr0),
    usTime = usTime, lsTime = lsTime
  )

  z <- gsnSurv(x, y$n.I[k])
  eDC <- NULL
  eDE <- NULL
  eDC0 <- NULL
  eDE0 <- NULL
  eNC <- NULL
  eNE <- NULL
  T <- NULL
  for (i in 1:(k - 1)) {
    xx <- tEventsIA(z, y$timing[i], tol)
    T <- c(T, xx$T)
    eDC <- rbind(eDC, xx$eDC)
    eDE <- rbind(eDE, xx$eDE)
    eNC <- rbind(eNC, xx$eNC)
    eNE <- rbind(eNE, xx$eNE)
  }
  y$T <- c(T, z$T)
  y$eDC <- rbind(eDC, z$eDC)
  y$eDE <- rbind(eDE, z$eDE)
  y$eNC <- rbind(eNC, z$eNC)
  y$eNE <- rbind(eNE, z$eNE)
  y$hr <- hr
  y$hr0 <- hr0
  y$R <- z$R
  y$S <- z$S
  y$minfup <- z$minfup
  y$gamma <- z$gamma
  y$ratio <- ratio
  y$lambdaC <- z$lambdaC
  y$etaC <- z$etaC
  y$etaE <- z$etaE
  y$variable <- x$variable
  y$tol <- tol
  y$method <- x$method
  y$call <- match.call()
  class(y) <- c("gsSurv", "gsDesign")

  nameR <- nameperiod(cumsum(y$R))
  stratnames <- paste("Stratum", 1:ncol(y$lambdaC))
  if (is.null(y$S)) {
    nameS <- "0-Inf"
  } else {
    nameS <- nameperiod(cumsum(c(y$S, Inf)))
  }
  rownames(y$lambdaC) <- nameS
  colnames(y$lambdaC) <- stratnames
  rownames(y$etaC) <- nameS
  colnames(y$etaC) <- stratnames
  rownames(y$etaE) <- nameS
  colnames(y$etaE) <- stratnames
  rownames(y$gamma) <- nameR
  colnames(y$gamma) <- stratnames
  return(y)
}
# print.gsSurv roxy [sinew] ----
#' @rdname nSurv
#' @export
# print.gsSurv function [sinew] ----
print.gsSurv <- function(x, digits = 3, show_gsDesign = FALSE, show_strata = TRUE, ...) {
  if (!inherits(x, "gsSurv")) {
    stop("Primary argument must have class gsSurv")
  }

  # Helper to preserve input labels
  label_input <- function(nm) {
    if (!is.null(x$call) && !is.null(x$call[[nm]])) {
      paste(deparse(x$call[[nm]]), collapse = "")
    } else {
      nm
    }
  }

  # Helper for compact display (matching print.nSurv style)
  compact <- function(v) {
    if (is.null(v)) return("NULL")
    # Convert matrix to vector for display (row-wise)
    if (is.matrix(v)) {
      v <- as.vector(t(v))
    }
    if (length(v) <= 3) return(paste(round(v, digits), collapse = ", "))
    paste0(
      paste(round(v[1:3], digits), collapse = ", "),
      " ... (len=", length(v), ")"
    )
  }

  # Test type description
  test_type_desc <- switch(as.character(x$test.type),
    "1" = "One-sided (efficacy only)",
    "2" = "Two-sided symmetric",
    "3" = "Two-sided asymmetric (efficacy and futility)",
    "4" = "Two-sided asymmetric with non-binding futility",
    "5" = "Two-sided asymmetric with binding futility",
    paste("Test type", x$test.type)
  )

  # Summary header
  cat(
    "Group sequential design ",
    "(method=", x$method, "; k=", x$k, " analyses; ", test_type_desc, ")\n",
    sep = ""
  )
  cat(
    sprintf(
      "HR=%.3f vs HR0=%.3f | alpha=%.3f (sided=%d) | power=%.1f%%\n",
      x$hr, x$hr0, x$alpha, x$sided, (1 - x$beta) * 100
    )
  )
  
  # Get final analysis values
  final_n <- sum((x$eNC + x$eNE)[x$k, ])
  final_d <- sum((x$eDC + x$eDE)[x$k, ])
  final_T <- x$T[x$k]
  accrual_dur <- if (!is.null(x$R) && length(x$R) > 0) {
    sum(x$R)
  } else {
    final_T - x$minfup
  }
  
  cat(
    sprintf(
      paste(
        "N=%.1f subjects | D=%.1f events |",
        "T=%.1f study duration |",
        "accrual=%.1f Accrual duration |",
        "minfup=%.1f minimum follow-up |",
        "ratio=%s randomization ratio (experimental/control)\n"
      ),
      final_n, final_d, final_T, accrual_dur, x$minfup,
      paste(x$ratio, collapse = "/")
    )
  )

  # Spending function summary (extracted from summary())
  summ_text <- summary(x)
  # Extract spending function sentences
  spending_parts <- unlist(strsplit(summ_text, "\\."))
  spending_parts <- grep("spending function", spending_parts, ignore.case = TRUE, value = TRUE)
  if (length(spending_parts) > 0) {
    cat("\nSpending functions:\n")
    for (part in spending_parts) {
      cat("  ", trimws(part), ".\n", sep = "")
    }
  }
  
  # Analysis summary table using gsBoundSummary
  cat("\nAnalysis summary:\n")
  print(gsDesign::gsBoundSummary(x))

  # Optional: show gsDesign details
  if (show_gsDesign) {
    cat("\n--- gsDesign details ---\n")
    gsDesign:::print.gsDesign(x)
  }

  # Key inputs table (matching print.nSurv style)
  cat("\nKey inputs (names preserved):\n")
  items <- c("gamma", "R", "lambdaC", "eta", "etaE", "S")
  desc_map <- list(
    gamma = "Accrual rate(s)",
    R = "Accrual rate duration(s)",
    lambdaC = "Control hazard rate(s)",
    eta = "Control dropout rate(s)",
    etaE = "Experimental dropout rate(s)",
    S = "Event and dropout rate duration(s)"
  )
  tab <- data.frame(
    desc = unlist(desc_map[items], use.names = FALSE),
    item = items,
    value = vapply(
      items,
      function(nm) {
        switch(nm,
          lambdaC = compact(x$lambdaC),
          gamma = compact(x$gamma),
          eta = compact(x$etaC),
          etaE = compact(x$etaE),
          R = compact(x$R),
          S = compact(x$S),
          compact(NULL)
        )
      },
      character(1)
    ),
    input = vapply(items, function(nm) {
      # For "eta", check both "eta" and "etaC" in the call
      if (nm == "eta") {
        if (!is.null(x$call) && "eta" %in% names(x$call)) {
          return(label_input("eta"))
        } else if (!is.null(x$call) && "etaC" %in% names(x$call)) {
          return(label_input("etaC"))
        } else {
          return("eta")
        }
      }
      label_input(nm)
    }, character(1)),
    stringsAsFactors = FALSE
  )
  print(tab, row.names = FALSE)

  # Optional: show strata details
  if (show_strata && !is.null(x$eDC) && ncol(x$eDC) > 1) {
    cat("\nExpected events by stratum (H1) at final analysis:\n")
    evt <- data.frame(
      stratum = paste("Stratum", seq_len(ncol(x$eDC))),
      control = round(x$eDC[x$k, ], digits),
      experimental = round(x$eDE[x$k, ], digits),
      total = round(x$eDC[x$k, ] + x$eDE[x$k, ], digits)
    )
    print(evt, row.names = FALSE)
  }

  invisible(x)
}

We also provide an xtable function for gsSurv object types

# xtable.gsSurv roxy [sinew] ----
#' @seealso
#'  \code{\link[stats]{Normal}}
#'  \code{\link[xtable]{xtable}}
#' @importFrom stats pnorm
#' @rdname nSurv
#' @export
# xtable.gsSurv function [sinew] ----
xtable.gsSurv <- function(x, caption = NULL, label = NULL, align = NULL, digits = NULL,
                          display = NULL, auto = FALSE, footnote = NULL, fnwid = "9cm", timename = "months", ...) {
  k <- x$k
  stat <- c(
    "Z-value", "HR", "p (1-sided)", paste("P\\{Cross\\} if HR=", x$hr0, sep = ""),
    paste("P\\{Cross\\} if HR=", x$hr, sep = "")
  )
  st <- stat
  for (i in 2:k) stat <- c(stat, st)
  an <- rep(" ", 5 * k)
  tim <- an
  enrol <- an
  fut <- an
  eff <- an
  an[5 * (0:(k - 1)) + 1] <- c(paste("IA ", as.character(1:(k - 1)), ": ",
    as.character(round(100 * x$timing[1:(k - 1)], 1)), "\\%",
    sep = ""
  ), "Final analysis")
  an[5 * (1:(k - 1)) + 1] <- paste("\\hline", an[5 * (1:(k - 1)) + 1])
  an[5 * (0:(k - 1)) + 2] <- paste("N:", ceiling(rowSums(x$eNC)) + ceiling(rowSums(x$eNE)))
  an[5 * (0:(k - 1)) + 3] <- paste("Events:", ceiling(rowSums(x$eDC + x$eDE)))
  an[5 * (0:(k - 1)) + 4] <- paste(round(x$T, 1), timename, sep = " ")
  if (x$test.type != 1) fut[5 * (0:(k - 1)) + 1] <- as.character(round(x$lower$bound, 2))
  eff[5 * (0:(k - 1)) + 1] <- as.character(round(x$upper$bound, 2))
  if (x$test.type != 1) fut[5 * (0:(k - 1)) + 2] <- as.character(round(gsHR(z = x$lower$bound, i = 1:k, x, ratio = x$ratio) * x$hr0, 2))
  eff[5 * (0:(k - 1)) + 2] <- as.character(round(gsHR(z = x$upper$bound, i = 1:k, x, ratio = x$ratio) * x$hr0, 2))
  asp <- as.character(round(stats::pnorm(-x$upper$bound), 4))
  asp[asp == "0"] <- "$< 0.0001$"
  eff[5 * (0:(k - 1)) + 3] <- asp
  asp <- as.character(round(cumsum(x$upper$prob[, 1]), 4))
  asp[asp == "0"] <- "$< 0.0001$"
  eff[5 * (0:(k - 1)) + 4] <- asp
  asp <- as.character(round(cumsum(x$upper$prob[, 2]), 4))
  asp[asp == "0"] <- "$< 0.0001$"
  eff[5 * (0:(k - 1)) + 5] <- asp
  if (x$test.type != 1) {
    bsp <- as.character(round(stats::pnorm(-x$lower$bound), 4))
    bsp[bsp == "0"] <- " $< 0.0001$"
    fut[5 * (0:(k - 1)) + 3] <- bsp
    bsp <- as.character(round(cumsum(x$lower$prob[, 1]), 4))
    bsp[bsp == "0"] <- "$< 0.0001$"
    fut[5 * (0:(k - 1)) + 4] <- bsp
    bsp <- as.character(round(cumsum(x$lower$prob[, 2]), 4))
    bsp[bsp == "0"] <- "$< 0.0001$"
    fut[5 * (0:(k - 1)) + 5] <- bsp
  }
  neff <- length(eff)
  if (!is.null(footnote)) {
    eff[neff] <-
      paste(eff[neff], "\\\\ \\hline \\multicolumn{4}{p{", fnwid, "}}{\\footnotesize", footnote, "}")
  }
  if (x$test.type != 1) {
    xxtab <- data.frame(cbind(an, stat, fut, eff))
    colnames(xxtab) <- c("Analysis", "Value", "Futility", "Efficacy")
  } else {
    xxtab <- data.frame(cbind(an, stat, eff))
    colnames(xxtab) <- c("Analysis", "Value", "Efficacy")
  }
  return(xtable::xtable(xxtab,
    caption = caption, label = label, align = align, digits = digits,
    display = display, auto = auto, ...
  ))
}

Note: gsSurv() and gsSurvCalendar() pass the method argument through to the fixed-design computation (nSurv()), and the returned objects now include $method to indicate which variance/event formulation was used. The same restrictions apply: "Schoenfeld" and "Freedman" methods only support superiority testing (hr0 = 1), and "Freedman" does not support stratified populations.

4.5 Examples

4.5.1 Simple

Fixed accrual duration and follow-up duration, variable accrual rate.

# `method` controls the fixed-design variance formulation:
# "LachinFoulkes" (default), "Schoenfeld", "Freedman", or "BernsteinLagakos".
x <- gsSurv(alpha = 0.025, beta = 0.1,
  lambdaC = log(2) / 6, hr = .6, minfup = 6, T = 36, test.type = 4,
  R = 30, sfu = sfHSD, sfupar = -4, sfl = sfHSD, sflpar = 1,
  method = "LachinFoulkes")
print(x)
Group sequential design (method=LachinFoulkes; k=3 analyses; Two-sided asymmetric with non-binding futility)
N=248.2 subjects | D=196.4 events | T=36.0 study duration | accrual=30.0 Accrual duration | minfup=6.0 minimum follow-up | ratio=1 randomization ratio (experimental/control)

Spending functions:
  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 = 1.

Analysis summary:
    Analysis              Value Efficacy Futility
   IA 1: 33%                  Z   3.0107   0.3779
      N: 138        p (1-sided)   0.0013   0.3528
  Events: 66       ~HR at bound   0.4751   0.9108
   Month: 17   P(Cross) if HR=1   0.0013   0.6472
             P(Cross) if HR=0.6   0.1747   0.0448
   IA 2: 67%                  Z   2.5465   1.2990
      N: 216        p (1-sided)   0.0054   0.0970
 Events: 131       ~HR at bound   0.6408   0.7969
   Month: 26   P(Cross) if HR=1   0.0062   0.9151
             P(Cross) if HR=0.6   0.6531   0.0770
       Final                  Z   1.9992   1.9992
      N: 250        p (1-sided)   0.0228   0.0228
 Events: 197       ~HR at bound   0.7518   0.7518
   Month: 36   P(Cross) if HR=1   0.0199   0.9801
             P(Cross) if HR=0.6   0.9000   0.1000

Key inputs (names preserved):
                               desc    item value    input
                    Accrual rate(s)   gamma 8.274    gamma
           Accrual rate duration(s)       R    30       30
             Control hazard rate(s) lambdaC 0.116 log(2)/6
            Control dropout rate(s)     eta     0      eta
       Experimental dropout rate(s)    etaE     0     etaE
 Event and dropout rate duration(s)       S  NULL        S

4.5.2 How T, R, and minfup interact to determine sample size

The parameters T (study duration), R (accrual duration), and minfup (minimum follow-up duration) interact to determine which variable is solved for in the sample size calculation:

  1. When T and minfup are both provided: The accrual rate (gamma) is solved for to achieve the required number of events. The accrual duration R is fixed at T - minfup.

  2. When T is NULL but minfup is provided: The accrual duration R is solved for to achieve the required number of events. The study duration T is then R + minfup.

  3. When both gamma and R are provided and T is NULL: The accrual rate (gamma) is solved for to achieve the required number of events, keeping R fixed. This matches the behavior of gsDesign::gsSurv(). The study duration T is set to R + minfup internally.

  4. When minfup is NULL: The minimum follow-up duration is solved for to achieve the required number of events, keeping R and T fixed.

The relationship between these parameters is: T = R + minfup (when all are specified). For group sequential designs, the required number of events is inflated from the fixed design, and the appropriate variable is adjusted proportionally to achieve the inflated event count.

4.5.3 Alternative methods: Schoenfeld and Freedman

The Schoenfeld method uses a simpler variance formulation based on the hazard ratio contrast. Note that Schoenfeld method only supports superiority testing (hr0 = 1).

x_schoen <- gsSurv(alpha = 0.025, beta = 0.1,
  lambdaC = log(2) / 6, hr = .6, minfup = 6, T = 36, test.type = 4,
  R = 30, sfu = sfHSD, sfupar = -4, sfl = sfHSD, sflpar = 1,
  method = "Schoenfeld"
)
print(x_schoen)
Group sequential design (method=Schoenfeld; k=3 analyses; Two-sided asymmetric with non-binding futility)
N=250.2 subjects | D=198.0 events | T=36.0 study duration | accrual=30.0 Accrual duration | minfup=6.0 minimum follow-up | ratio=1 randomization ratio (experimental/control)

Spending functions:
  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 = 1.

Analysis summary:
    Analysis              Value Efficacy Futility
   IA 1: 33%                  Z   3.0107   0.3779
      N: 138        p (1-sided)   0.0013   0.3528
  Events: 66       ~HR at bound   0.4765   0.9112
   Month: 17   P(Cross) if HR=1   0.0013   0.6472
             P(Cross) if HR=0.6   0.1747   0.0448
   IA 2: 67%                  Z   2.5465   1.2990
      N: 218        p (1-sided)   0.0054   0.0970
 Events: 132       ~HR at bound   0.6419   0.7976
   Month: 26   P(Cross) if HR=1   0.0062   0.9151
             P(Cross) if HR=0.6   0.6531   0.0770
       Final                  Z   1.9992   1.9992
      N: 252        p (1-sided)   0.0228   0.0228
 Events: 198       ~HR at bound   0.7526   0.7526
   Month: 36   P(Cross) if HR=1   0.0199   0.9801
             P(Cross) if HR=0.6   0.9000   0.1000

Key inputs (names preserved):
                               desc    item value    input
                    Accrual rate(s)   gamma  8.34    gamma
           Accrual rate duration(s)       R    30       30
             Control hazard rate(s) lambdaC 0.116 log(2)/6
            Control dropout rate(s)     eta     0      eta
       Experimental dropout rate(s)    etaE     0     etaE
 Event and dropout rate duration(s)       S  NULL        S

Comparison with rpact for Schoenfeld method (non-binding futility):

suppressPackageStartupMessages({
  if (!require("rpact", quietly = TRUE)) {
    install.packages("rpact", repos = "https://cran.rstudio.com/")
  }
  library(rpact)
})

# Match gsSurv parameters: k=3, alpha=0.025, sided=1, beta=0.1
# Non-binding futility (test.type=4)
# Match spending functions: sfHSD with gamma=-4 (upper) and gamma=1 (lower)
design_gs <- getDesignGroupSequential(
  kMax = 3, alpha = 0.025, beta = 0.1, sided = 1,
  informationRates = x_schoen$timing,
  typeOfDesign = "asHSD",
  gammaA = -4,
  typeBetaSpending = "bsHSD",
  gammaB = 1
)

# Get accrual and follow-up from gsSurv result
accrual_time <- sum(x_schoen$R)
followup_time <- x_schoen$minfup
lambda2 <- log(2) / 6

rp_schoen_gs <- getSampleSizeSurvival(
  design = design_gs,
  lambda2 = lambda2, hazardRatio = 0.6,
  accrualTime = accrual_time, followUpTime = followup_time,
  dropoutRate1 = 0, dropoutRate2 = 0
)

# Compare at each analysis: sample size, events, z-value bounds, and calendar timing
library(gt)
n_gsSurv <- as.vector((x_schoen$eNC + x_schoen$eNE) %*% rep(1, ncol(x_schoen$eNE)))
events_gsSurv <- as.vector((x_schoen$eDC + x_schoen$eDE) %*% rep(1, ncol(x_schoen$eDE)))
n_rpact <- as.vector(rp_schoen_gs$numberOfSubjects[, 1])
events_rpact <- as.vector(rp_schoen_gs$cumulativeEventsPerStage[, 1])
z_lower_rpact <- c(design_gs$futilityBounds, design_gs$criticalValues[3])
time_gsSurv <- x_schoen$T
time_rpact <- as.vector(rp_schoen_gs$analysisTime[, 1])
comparison_schoen <- data.frame(
  Analysis = c("IA 1", "IA 2", "Final"),
  `N (gsSurv)` = round(n_gsSurv, 1),
  `N (rpact)` = round(n_rpact, 1),
  `Events (gsSurv)` = round(events_gsSurv, 1),
  `Events (rpact)` = round(events_rpact, 1),
  `Time (gsSurv)` = round(time_gsSurv, 1),
  `Time (rpact)` = round(time_rpact, 1),
  `Z upper (gsSurv)` = round(x_schoen$upper$bound, 3),
  `Z upper (rpact)` = round(design_gs$criticalValues, 3),
  `Z lower (gsSurv)` = round(x_schoen$lower$bound, 3),
  `Z lower (rpact)` = round(z_lower_rpact, 3),
  check.names = FALSE
)
comparison_schoen |>
  gt() |>
  tab_header(
    title = "Schoenfeld method comparison",
    subtitle = "gsSurv vs rpact (non-binding futility)"
  ) |>
  fmt_number(columns = c(2:5), decimals = 1) |>
  fmt_number(columns = c(6:7), decimals = 1) |>
  fmt_number(columns = c(8:11), decimals = 3)
Schoenfeld method comparison
gsSurv vs rpact (non-binding futility)
Analysis N (gsSurv) N (rpact) Events (gsSurv) Events (rpact) Time (gsSurv) Time (rpact) Z upper (gsSurv) Z upper (rpact) Z lower (gsSurv) Z lower (rpact)
IA 1 137.7 137.8 66.0 66.0 16.5 16.5 3.011 3.011 0.378 0.378
IA 2 216.5 216.5 132.0 132.0 26.0 26.0 2.547 2.547 1.299 1.299
Final 250.2 250.2 198.0 198.0 36.0 36.0 1.999 1.999 1.999 1.999

The Freedman method uses a different variance formulation based on event proportions. Note that Freedman method only supports superiority testing (hr0 = 1) and does not support stratified populations.

x_freed <- gsSurv(alpha = 0.025, beta = 0.1,
  lambdaC = log(2) / 6, hr = .6, minfup = 6, T = 36, test.type = 4,
  R = 30, sfu = sfHSD, sfupar = -4, sfl = sfHSD, sflpar = 1,
  method =  "Freedman")
print(x_freed)
Group sequential design (method=Freedman; k=3 analyses; Two-sided asymmetric with non-binding futility)
N=261.2 subjects | D=206.7 events | T=36.0 study duration | accrual=30.0 Accrual duration | minfup=6.0 minimum follow-up | ratio=1 randomization ratio (experimental/control)

Spending functions:
  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 = 1.

Analysis summary:
    Analysis              Value Efficacy Futility
   IA 1: 33%                  Z   3.0107   0.3779
      N: 144        p (1-sided)   0.0013   0.3528
  Events: 69       ~HR at bound   0.4841   0.9130
   Month: 17   P(Cross) if HR=1   0.0013   0.6472
             P(Cross) if HR=0.6   0.1747   0.0448
   IA 2: 67%                  Z   2.5465   1.2990
      N: 226        p (1-sided)   0.0054   0.0970
 Events: 138       ~HR at bound   0.6480   0.8014
   Month: 26   P(Cross) if HR=1   0.0062   0.9151
             P(Cross) if HR=0.6   0.6531   0.0770
       Final                  Z   1.9992   1.9992
      N: 262        p (1-sided)   0.0228   0.0228
 Events: 207       ~HR at bound   0.7572   0.7572
   Month: 36   P(Cross) if HR=1   0.0199   0.9801
             P(Cross) if HR=0.6   0.9000   0.1000

Key inputs (names preserved):
                               desc    item value    input
                    Accrual rate(s)   gamma 8.705    gamma
           Accrual rate duration(s)       R    30       30
             Control hazard rate(s) lambdaC 0.116 log(2)/6
            Control dropout rate(s)     eta     0      eta
       Experimental dropout rate(s)    etaE     0     etaE
 Event and dropout rate duration(s)       S  NULL        S

Comparison with rpact for Freedman method (non-binding futility):

# Get accrual and follow-up from gsSurv result
accrual_time <- sum(x_freed$R)
followup_time <- x_freed$minfup

# Use the same design as for Schoenfeld comparison
# (design_gs was created in the previous chunk)
rp_freed_gs <- getSampleSizeSurvival(
  design = design_gs,
  lambda2 = lambda2, hazardRatio = 0.6,
  accrualTime = accrual_time, followUpTime = followup_time,
  dropoutRate1 = 0, dropoutRate2 = 0,
  typeOfComputation = "Freedman"
)

# Compare at each analysis: sample size, events, z-value bounds, and calendar timing
n_gsSurv <- as.vector((x_freed$eNC + x_freed$eNE) %*% rep(1, ncol(x_freed$eNE)))
events_gsSurv <- as.vector((x_freed$eDC + x_freed$eDE) %*% rep(1, ncol(x_freed$eDE)))
n_rpact <- as.vector(rp_freed_gs$numberOfSubjects[, 1])
events_rpact <- as.vector(rp_freed_gs$cumulativeEventsPerStage[, 1])
z_lower_rpact <- c(design_gs$futilityBounds, design_gs$criticalValues[3])
time_gsSurv <- x_freed$T
time_rpact <- as.vector(rp_freed_gs$analysisTime[, 1])
comparison_freed <- data.frame(
  Analysis = c("IA 1", "IA 2", "Final"),
  `N (gsSurv)` = round(n_gsSurv, 1),
  `N (rpact)` = round(n_rpact, 1),
  `Events (gsSurv)` = round(events_gsSurv, 1),
  `Events (rpact)` = round(events_rpact, 1),
  `Time (gsSurv)` = round(time_gsSurv, 1),
  `Time (rpact)` = round(time_rpact, 1),
  `Z upper (gsSurv)` = round(x_freed$upper$bound, 3),
  `Z upper (rpact)` = round(design_gs$criticalValues, 3),
  `Z lower (gsSurv)` = round(x_freed$lower$bound, 3),
  `Z lower (rpact)` = round(z_lower_rpact, 3),
  check.names = FALSE
)
comparison_freed |>
  gt() |>
  tab_header(
    title = "Freedman method comparison",
    subtitle = "gsSurv vs rpact (non-binding futility)"
  ) |>
  fmt_number(columns = c(2:5), decimals = 1) |>
  fmt_number(columns = c(6:7), decimals = 1) |>
  fmt_number(columns = c(8:11), decimals = 3)
Freedman method comparison
gsSurv vs rpact (non-binding futility)
Analysis N (gsSurv) N (rpact) Events (gsSurv) Events (rpact) Time (gsSurv) Time (rpact) Z upper (gsSurv) Z upper (rpact) Z lower (gsSurv) Z lower (rpact)
IA 1 143.8 143.8 68.9 68.9 16.5 16.5 3.011 3.011 0.378 0.378
IA 2 226.0 226.0 137.8 137.8 26.0 26.0 2.547 2.547 1.299 1.299
Final 261.2 261.2 206.7 206.7 36.0 36.0 1.999 1.999 1.999 1.999

4.5.4 Piecewise enrollment and failure rates

x <- gsSurv(
  lambdaC = log(2) / c(6, 8, 10), hr = .6, gamma = c(2, 4),
  S = c(3, 6), R = c(3, 3), minfup = 6, T = 20
)
print(x)
Group sequential design (method=LachinFoulkes; k=3 analyses; Two-sided asymmetric with non-binding futility)
N=306.0 subjects | D=172.3 events | T=20.0 study duration | accrual=14.0 Accrual duration | minfup=6.0 minimum follow-up | ratio=1 randomization ratio (experimental/control)

Spending functions:
  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.

Analysis summary:
    Analysis              Value Efficacy Futility
   IA 1: 33%                  Z   3.0107  -0.2388
      N: 204        p (1-sided)   0.0013   0.5944
  Events: 58       ~HR at bound   0.4518   1.0650
   Month: 10   P(Cross) if HR=1   0.0013   0.4056
             P(Cross) if HR=0.6   0.1412   0.0148
   IA 2: 67%                  Z   2.5465   0.9410
      N: 308        p (1-sided)   0.0054   0.1733
 Events: 115       ~HR at bound   0.6217   0.8389
   Month: 14   P(Cross) if HR=1   0.0062   0.8347
             P(Cross) if HR=0.6   0.5815   0.0437
       Final                  Z   1.9992   1.9992
      N: 308        p (1-sided)   0.0228   0.0228
 Events: 173       ~HR at bound   0.7374   0.7374
   Month: 20   P(Cross) if HR=1   0.0233   0.9767
             P(Cross) if HR=0.6   0.9000   0.1000

Key inputs (names preserved):
                               desc    item               value
                    Accrual rate(s)   gamma      12.242, 24.483
           Accrual rate duration(s)       R               3, 11
             Control hazard rate(s) lambdaC 0.116, 0.087, 0.069
            Control dropout rate(s)     eta             0, 0, 0
       Experimental dropout rate(s)    etaE             0, 0, 0
 Event and dropout rate duration(s)       S                3, 6
              input
            c(2, 4)
            c(3, 3)
 log(2)/c(6, 8, 10)
                eta
               etaE
            c(3, 6)

4.5.5 Stratified population

We begin with the example of Bernstein and Lagakos (1978) using the LFPWE() routine. As with the Lachin and Foulkes (1986) variances are computed under both the null and alternate hypotheses. For the Bernstein and Lagakos (1978) method, the null variance is computed using alternate hypothesis control group event for both the control and experimental groups; this is for superiority trials as a slightly different approach is used for non-inferiority and super-superiority designs. For the Lachin and Foulkes (1986) method the null hypothesis chosen for superiority designs has an average of the failure rates in the control and experimental arms. The Bernstein and Lagakos (1978) method proposes a sample size of 172 and 150 expected events. With gsDesign2::gs_designAHR() with a 2-variance formulation comparable to here, the results are N=182 (rounded up to nearest even integer) and 151 events (rounded up to nearest integer).

LFPWE(
  alpha = .05, beta = .2, hr = 2 / 3, ratio = 1,
  R = 2, S = NULL, T = 4,
  lambdaC = matrix(c(1, .8, .5), nrow = 1),
  etaC = matrix(0, nrow = 1, ncol = 3),
  etaE = matrix(0, nrow = 1, ncol = 3),
  gamma = matrix(c(.4, .4, .2), nrow = 1)
)
Fixed design, two-arm trial with time-to-event
outcome (Lachin and Foulkes, 1986).
Solving for:  Accrual rate 
Hazard ratio                  H1/H0=0.6667/1
Study duration:                   T=4
Accrual duration:                   2
Min. end-of-study follow-up: minfup=2
Expected events (total, H1):        149.4726
Expected sample size (total):       178.797
Accrual rates:
        [,1]    [,2]    [,3]
[1,] 35.7594 35.7594 17.8797
Control event rates (H1):
     [,1] [,2] [,3]
[1,]    1  0.8  0.5
Censoring rates:
     [,1] [,2] [,3]
[1,]    0    0    0
Power:                 100*(1-beta)=80%
Type I error (1-sided):   100*alpha=5%
Equal randomization:          ratio=1
x <- gsSurv(
  alpha = .05, beta = .2, lambdaC = matrix(c(1, .8, .5), nrow = 1),
  hr = 2 / 3, gamma = matrix(c(.4, .4, .2), nrow = 1), minfup = 2, S = NULL, T = 4, ratio = 1,
  method = "BernsteinLagakos"
)
print(x)
Group sequential design (method=BernsteinLagakos; k=3 analyses; Two-sided asymmetric with non-binding futility)
N=185.8 subjects | D=155.3 events | T=4.0 study duration | accrual=2.0 Accrual duration | minfup=2.0 minimum follow-up | ratio=1 randomization ratio (experimental/control)

Spending functions:
  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.

Analysis summary:
    Analysis               Value Efficacy Futility
   IA 1: 33%                   Z   2.7936  -0.3978
      N: 142         p (1-sided)   0.0026   0.6546
  Events: 52        ~HR at bound   0.4600   1.1169
    Month: 2    P(Cross) if HR=1   0.0026   0.3454
             P(Cross) if HR=0.67   0.0958   0.0297
   IA 2: 67%                   Z   2.2890   0.6660
      N: 186         p (1-sided)   0.0110   0.2527
 Events: 104        ~HR at bound   0.6377   0.8773
    Month: 2    P(Cross) if HR=1   0.0125   0.7592
             P(Cross) if HR=0.67   0.4340   0.0875
       Final                   Z   1.6798   1.6798
      N: 186         p (1-sided)   0.0465   0.0465
 Events: 156        ~HR at bound   0.7637   0.7637
    Month: 4    P(Cross) if HR=1   0.0468   0.9532
             P(Cross) if HR=0.67   0.8000   0.2000

Key inputs (names preserved):
                               desc    item                  value
                    Accrual rate(s)   gamma 37.162, 37.162, 18.581
           Accrual rate duration(s)       R                      2
             Control hazard rate(s) lambdaC            1, 0.8, 0.5
            Control dropout rate(s)     eta                0, 0, 0
       Experimental dropout rate(s)    etaE                0, 0, 0
 Event and dropout rate duration(s)       S                   NULL
                              input
 matrix(c(0.4, 0.4, 0.2), nrow = 1)
                                  R
   matrix(c(1, 0.8, 0.5), nrow = 1)
                                eta
                               etaE
                                  S

Expected events by stratum (H1) at final analysis:
   stratum control experimental  total
 Stratum 1  34.987       31.751 66.739
 Stratum 2  33.419       29.298 62.717
 Stratum 3  14.260       11.618 25.878

4.5.6 Stratified population, piecewise enrollment and failure

Here we have three strata (columns of lambdaC and gamma), two enrollment rate periods (\(\leq 3\) months and \(>3\) months in R; rates in rows of gamma) and three failure rate periods (0–3, 3–6, > 6 months); rates in rows of lambdaC). We also have the default of 3 analyses (k = 3) and default spending functions. Here are the failure rate assumptions. Rows represent time periods and columns represent strata.

lambdaC_stratified <- log(2) / matrix(c(3, 4, 5, 6, 8, 10, 9, 12, 15), nrow = 3)
lambdaC_stratified
          [,1]       [,2]       [,3]
[1,] 0.2310491 0.11552453 0.07701635
[2,] 0.1732868 0.08664340 0.05776227
[3,] 0.1386294 0.06931472 0.04620981

Associated with these event rates are durations of time periods for each rate. Again, rows represent time periods and columns represent strata. Note that this is one less row here as the final event rate above issumed to continue indefinitely. Notice also that the durations of rates can, as here, be different for different strata.

gamma_stratified <- matrix(c(2, 4, 8, 3, 6, 10), nrow = 2, byrow = TRUE)
gamma_stratified
     [,1] [,2] [,3]
[1,]    2    4    8
[2,]    3    6   10

We begin with the case where accrual duration and minimum follow-up duration are fixed and we solve for accrual rates.

x <- gsSurv(
  lambdaC = log(2) / matrix(c(3, 4, 5, 6, 8, 10, 9, 12, 15), nrow = 3),
  gamma = matrix(c(2, 4, 8, 3, 6, 10), nrow = 2, byrow = TRUE),
  hr = .6, S = c(3, 6), R = c(3, 3), minfup = 6, T = 30, ratio = 2, hr0 = 1.1
)
print(x)
Group sequential design (method=LachinFoulkes; k=3 analyses; Two-sided asymmetric with non-binding futility)
N=228.4 subjects | D=134.0 events | T=30.0 study duration | accrual=24.0 Accrual duration | minfup=6.0 minimum follow-up | ratio=2 randomization ratio (experimental/control)

Spending functions:
  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.

Analysis summary:
    Analysis              Value Efficacy Futility
   IA 1: 33%                  Z   3.0107  -0.2388
      N: 134        p (1-sided)   0.0013   0.5944
  Events: 45       ~HR at bound   0.4230   1.1866
   Month: 14 P(Cross) if HR=1.1   0.0013   0.4056
             P(Cross) if HR=0.6   0.1412   0.0148
   IA 2: 67%                  Z   2.5465   0.9410
      N: 209        p (1-sided)   0.0054   0.1733
  Events: 90       ~HR at bound   0.6211   0.8906
   Month: 22 P(Cross) if HR=1.1   0.0062   0.8347
             P(Cross) if HR=0.6   0.5815   0.0437
       Final                  Z   1.9992   1.9992
      N: 230        p (1-sided)   0.0228   0.0228
 Events: 135       ~HR at bound   0.7626   0.7626
   Month: 30 P(Cross) if HR=1.1   0.0233   0.9767
             P(Cross) if HR=0.6   0.9000   0.1000

Key inputs (names preserved):
                               desc    item                           value
                    Accrual rate(s)   gamma 1.036, 2.072, 4.144 ... (len=6)
           Accrual rate duration(s)       R                           3, 21
             Control hazard rate(s) lambdaC 0.231, 0.116, 0.077 ... (len=9)
            Control dropout rate(s)     eta             0, 0, 0 ... (len=9)
       Experimental dropout rate(s)    etaE             0, 0, 0 ... (len=9)
 Event and dropout rate duration(s)       S                            3, 6
                                                    input
     matrix(c(2, 4, 8, 3, 6, 10), nrow = 2, byrow = TRUE)
                                                  c(3, 3)
 log(2)/matrix(c(3, 4, 5, 6, 8, 10, 9, 12, 15), nrow = 3)
                                                      eta
                                                     etaE
                                                  c(3, 6)

Expected events by stratum (H1) at final analysis:
   stratum control experimental  total
 Stratum 1  10.936       18.976 29.911
 Stratum 2  17.631       26.950 44.581
 Stratum 3  24.442       35.071 59.514

The total number of enrollment and events at each analysis is printed above. Following we see the enrollment and expected number of events by stratum in rows of the following matrices.

x$eNC + x$eNE
         [,1]     [,2]      [,3]
[1,] 20.67854 41.35707  71.00051
[2,] 32.46600 64.93201 110.29207
[3,] 35.74297 71.48593 121.21528
x$eDC + x$eDE
         [,1]     [,2]     [,3]
[1,] 11.40395 14.60229 18.66238
[2,] 21.49967 29.54911 38.28839
[3,] 29.91126 44.58099 59.51361

We proceed to the case where accrual rate and minimum follow-up duration are fixed and solve for accrual duration.

x <- gsSurv(
  lambdaC = log(2) / matrix(c(3, 4, 5, 6, 8, 10, 9, 12, 15), nrow = 3),
  gamma = matrix(c(2, 4, 8, 3, 6, 10), nrow = 2, byrow = TRUE),
  hr = .6, S = c(3, 6), R = c(3, 3), minfup = 6, T = NULL
)
print(x)
Group sequential design (method=LachinFoulkes; k=3 analyses; Two-sided asymmetric with non-binding futility)
N=387.7 subjects | D=173.0 events | T=12.0 study duration | accrual=6.0 Accrual duration | minfup=6.0 minimum follow-up | ratio=1 randomization ratio (experimental/control)

Spending functions:
  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.

Analysis summary:
    Analysis              Value Efficacy Futility
   IA 1: 33%                  Z   3.0107  -0.2388
      N: 332        p (1-sided)   0.0013   0.5944
  Events: 58       ~HR at bound   0.4525   1.0649
    Month: 5   P(Cross) if HR=1   0.0013   0.4056
             P(Cross) if HR=0.6   0.1412   0.0148
   IA 2: 67%                  Z   2.5465   0.9410
      N: 388        p (1-sided)   0.0054   0.1733
 Events: 116       ~HR at bound   0.6224   0.8392
    Month: 8   P(Cross) if HR=1   0.0062   0.8347
             P(Cross) if HR=0.6   0.5815   0.0437
       Final                  Z   1.9992   1.9992
      N: 388        p (1-sided)   0.0228   0.0228
 Events: 173       ~HR at bound   0.7379   0.7379
   Month: 12   P(Cross) if HR=1   0.0233   0.9767
             P(Cross) if HR=0.6   0.9000   0.1000

Key inputs (names preserved):
                               desc    item                             value
                    Accrual rate(s)   gamma 7.833, 15.666, 31.332 ... (len=6)
           Accrual rate duration(s)       R                              3, 3
             Control hazard rate(s) lambdaC   0.231, 0.116, 0.077 ... (len=9)
            Control dropout rate(s)     eta               0, 0, 0 ... (len=9)
       Experimental dropout rate(s)    etaE               0, 0, 0 ... (len=9)
 Event and dropout rate duration(s)       S                              3, 6
                                                    input
     matrix(c(2, 4, 8, 3, 6, 10), nrow = 2, byrow = TRUE)
                                                  c(3, 3)
 log(2)/matrix(c(3, 4, 5, 6, 8, 10, 9, 12, 15), nrow = 3)
                                                      eta
                                                     etaE
                                                  c(3, 6)

Expected events by stratum (H1) at final analysis:
   stratum control experimental  total
 Stratum 1  23.588       18.385 41.973
 Stratum 2  32.902       22.929 55.831
 Stratum 3  45.115       30.080 75.195

4.5.7 Non-inferiority design

Non-inferiority designs test whether the experimental treatment is not substantially worse than control. Here we demonstrate a non-inferiority design with a null hypothesis hazard ratio of 1.3 (non-inferiority margin) and an expected hazard ratio of 1.0 (no difference between treatments). We use a single interim analysis at 50% of information (timing = 0.5).

x <- gsSurv(
  lambdaC = log(2) / 8, hr = 1.0, hr0 = 1.3,
  gamma = 8, R = 18, T = 30, minfup = 6,
  k = 2, timing = 0.5
)
print(x)
Group sequential design (method=LachinFoulkes; k=2 analyses; Two-sided asymmetric with non-binding futility)
N=853.4 subjects | D=639.9 events | T=30.0 study duration | accrual=24.0 Accrual duration | minfup=6.0 minimum follow-up | ratio=1 randomization ratio (experimental/control)

Spending functions:
  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.

Analysis summary:
    Analysis              Value Efficacy Futility
   IA 1: 50%                  Z   2.7500   0.4122
      N: 646        p (1-sided)   0.0030   0.3401
 Events: 320       ~HR at bound   0.9559   1.2414
   Month: 18 P(Cross) if HR=1.3   0.0030   0.6599
               P(Cross) if HR=1   0.3412   0.0269
       Final                  Z   1.9811   1.9811
      N: 854        p (1-sided)   0.0238   0.0238
 Events: 640       ~HR at bound   1.1115   1.1115
   Month: 30 P(Cross) if HR=1.3   0.0239   0.9761
               P(Cross) if HR=1   0.9000   0.1000

Key inputs (names preserved):
                               desc    item  value    input
                    Accrual rate(s)   gamma 35.558        8
           Accrual rate duration(s)       R     24       18
             Control hazard rate(s) lambdaC  0.087 log(2)/8
            Control dropout rate(s)     eta      0      eta
       Experimental dropout rate(s)    etaE      0     etaE
 Event and dropout rate duration(s)       S   NULL        S

4.6 Output options

4.6.2 Plot options

The gsDesign package provides plotting capabilities for group sequential designs. Since gsSurv objects inherit from gsDesign, you can use standard gsDesign plotting functions:

  • plot(): Creates boundary plots showing Z-value boundaries over information fraction
  • plot.gsDesign(): Additional plotting options for gsDesign objects
# Boundary plot
x <- gsSurv(lambdaC = log(2) / 6, hr = .6, gamma = 8, minfup = 6)
plot(x)

# Additional plotting options are available through gsDesign functions
# See ?plot.gsDesign for more details