3  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.

3.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. 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. This also increases the expected sample size by the same proportion. The other two computations are simple applications of the function KTZ().

3.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
)
x
$alpha
[1] 0.025

$sided
[1] 1

$beta
[1] 0.09999999

$power
[1] 0.9

$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] 21.53628

$minfup
[1] 3.536276

$hr
[1] 0.6

$hr0
[1] 1

$n
[1] 303

$d
[1] 160.8979

$tol
[1] 0.0001220703

$eDC
[1] 27.79551 23.34026 41.29445

$eDE
[1] 22.71825 17.33116 28.41828

$eDC0
[1] 25.75003 20.70618 35.33690

$eDE0
[1] 25.75003 20.70618 35.33690

$eNC
[1] 33.0 34.5 84.0

$eNE
[1] 33.0 34.5 84.0

$variable
[1] "Follow-up duration"

attr(,"class")
[1] "nSurv"
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"

3.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 <- stats::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
    ))
  }
}

3.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().

# 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) # KA: last 2 arguments added 10/8/2017
{
  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
  )
  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
  ) # KA: last 2 arguments added 10/8/2017

  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)
    # Following a placeholder
    #    eDC0 <- rbind(eDC0, xx$eDC0) # Added 6/15/2023
    #    eDE0 <- rbind(eDE0, xx$eDE0) # Added 6/15/2023
    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)
  # Following is a placeholder
  #  y$eDC0 <- rbind(eDC0, z$eDC0) # Added 6/15/2023
  #  y$eDE0 <- rbind(eDE0, z$eDE0) # Added 6/15/2023
  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
  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 = 2, ...) {
  cat("Time to event group sequential design with HR=", x$hr, "\n")
  if (x$hr0 != 1) cat("Non-inferiority design with null HR=", x$hr0, "\n")
  if (min(x$ratio == 1) == 1) {
    cat("Equal randomization:          ratio=1\n")
  } else {
    cat(
      "Randomization (Exp/Control):  ratio=",
      x$ratio, "\n"
    )
    if (length(x$ratio) > 1) cat("(randomization ratios shown by strata)\n")
  }
  gsDesign:::print.gsDesign(x)
  if (x$test.type != 1) {
    y <- cbind(
      x$T, (x$eNC + x$eNE) %*% rep(1, ncol(x$eNE)),
      (x$eDC + x$eDE) %*% rep(1, ncol(x$eNE)),
      round(gsDesign::zn2hr(x$lower$bound, x$n.I, x$ratio, hr0 = x$hr0, hr1 = x$hr), 3),
      round(gsDesign::zn2hr(x$upper$bound, x$n.I, x$ratio, hr0 = x$hr0, hr1 = x$hr), 3)
    )
    colnames(y) <- c("T", "n", "Events", "HR futility", "HR efficacy")
  } else {
    y <- cbind(
      x$T, (x$eNC + x$eNE) %*% rep(1, ncol(x$eNE)),
      (x$eDC + x$eDE) %*% rep(1, ncol(x$eNE)),
      round(gsDesign::zn2hr(x$upper$bound, x$n.I, x$ratio, hr0 = x$hr0, hr1 = x$hr), 3)
    )
    colnames(y) <- c("T", "n", "Events", "HR efficacy")
  }
  rnames <- paste("IA", 1:(x$k))
  rnames[length(rnames)] <- "Final"
  rownames(y) <- rnames
  print(y)
  cat("Accrual rates:\n")
  print(round(x$gamma, digits))
  cat("Control event rates (H1):\n")
  print(round(x$lambda, digits))
  if (max(abs(x$etaC - x$etaE)) == 0) {
    cat("Censoring rates:\n")
    print(round(x$etaC, digits))
  } else {
    cat("Control censoring rates:\n")
    print(round(x$etaC, digits))
    cat("Experimental censoring rates:\n")
    print(round(x$etaE, digits))
  }
}

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, ...
  ))
}

3.5 Examples

3.5.1 Simple

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

w <- LFPWE(lambdaC = log(2) / 3, hr = .6, minfup = 6, T = 36)
y <- nSurv(lambdaC = log(2) / 3, hr = .6, minfup = 6, T = 36, R = Inf)
x <- gsSurv(lambdaC = log(2) / 3, hr = .6, minfup = 6, T = 36)

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

x <- gsSurv(lambdaC = log(2) / 6, hr = .6, gamma = 8, minfup = 6)
x
Time to event group sequential design with HR= 0.6 
Equal randomization:          ratio=1
Asymmetric two-sided group sequential design with
90 % power and 2.5 % Type I Error.
Upper bound spending computations assume
trial continues if lower bound is crossed.

                 ----Lower bounds----  ----Upper bounds-----
  Analysis  N    Z   Nominal p Spend+  Z   Nominal p Spend++
         1  58 -0.24    0.4056 0.0148 3.01    0.0013  0.0013
         2 115  0.94    0.8267 0.0289 2.55    0.0054  0.0049
         3 172  2.00    0.9772 0.0563 2.00    0.0228  0.0188
     Total                     0.1000                 0.0250 
+ lower bound beta spending (under H1):
 Hwang-Shih-DeCani spending function with gamma = -2.
++ alpha spending:
 Hwang-Shih-DeCani spending function with gamma = -4.

Boundary crossing probabilities and expected sample size
assume any cross stops the trial

Upper boundary (power or Type I Error)
          Analysis
   Theta      1      2      3  Total  E{N}
  0.0000 0.0013 0.0049 0.0171 0.0233  99.9
  0.2564 0.1412 0.4403 0.3185 0.9000 126.5

Lower boundary (futility or Type II Error)
          Analysis
   Theta      1      2      3  Total
  0.0000 0.4056 0.4290 0.1420 0.9767
  0.2564 0.0148 0.0289 0.0563 0.1000
             T        n    Events HR futility HR efficacy
IA 1  15.47476 123.7981  57.00202       1.065       0.450
IA 2  24.17700 193.4160 114.00405       0.838       0.621
Final 33.50127 220.0102 171.00607       0.737       0.737
Accrual rates:
       Stratum 1
0-27.5         8
Control event rates (H1):
      Stratum 1
0-Inf      0.12
Censoring rates:
      Stratum 1
0-Inf         0

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

x <- gsSurv(lambdaC = log(2) / 6, hr = .6, gamma = 8, R = 25, minfup = NULL)

3.5.2 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
)
x
Time to event group sequential design with HR= 0.6 
Equal randomization:          ratio=1
Asymmetric two-sided group sequential design with
90 % power and 2.5 % Type I Error.
Upper bound spending computations assume
trial continues if lower bound is crossed.

                 ----Lower bounds----  ----Upper bounds-----
  Analysis  N    Z   Nominal p Spend+  Z   Nominal p Spend++
         1  58 -0.24    0.4056 0.0148 3.01    0.0013  0.0013
         2 115  0.94    0.8267 0.0289 2.55    0.0054  0.0049
         3 173  2.00    0.9772 0.0563 2.00    0.0228  0.0188
     Total                     0.1000                 0.0250 
+ lower bound beta spending (under H1):
 Hwang-Shih-DeCani spending function with gamma = -2.
++ alpha spending:
 Hwang-Shih-DeCani spending function with gamma = -4.

Boundary crossing probabilities and expected sample size
assume any cross stops the trial

Upper boundary (power or Type I Error)
          Analysis
   Theta      1      2      3  Total  E{N}
  0.0000 0.0013 0.0049 0.0171 0.0233 100.6
  0.2554 0.1412 0.4403 0.3185 0.9000 127.4

Lower boundary (futility or Type II Error)
          Analysis
   Theta      1      2      3  Total
  0.0000 0.4056 0.4290 0.1420 0.9767
  0.2554 0.0148 0.0289 0.0563 0.1000
              T        n    Events HR futility HR efficacy
IA 1   9.827039 203.8729  57.42358       1.065       0.452
IA 2  14.277264 306.0405 114.84716       0.839       0.622
Final 20.000000 306.0405 172.27073       0.737       0.737
Accrual rates:
     Stratum 1
0-3      12.24
3-14     24.48
Control event rates (H1):
      Stratum 1
0-3        0.12
3-9        0.09
9-Inf      0.07
Censoring rates:
      Stratum 1
0-3           0
3-9           0
9-Inf         0

3.5.3 Stratified population

We begin with the example of Bernstein and Lagakos (1978) using the LFPWE() routine and compare to the LFPWEnew() routine. Their results with a single (H1) variance formulation at 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 and 151 events.

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.455
Expected sample size (total):       178.7758
Accrual rates:
        [,1]    [,2]    [,3]
[1,] 35.7552 35.7552 17.8776
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

Now compare to new version.

# NOT READY
LFPWEnew(
  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)
)
x <- gsSurv(
  alpha = .05, beta = .2, lambda = 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:3
)
x
Time to event group sequential design with HR= 0.6666667 
Randomization (Exp/Control):  ratio= 1 2 3 
(randomization ratios shown by strata)
Asymmetric two-sided group sequential design with
80 % power and 5 % Type I Error.
Upper bound spending computations assume
trial continues if lower bound is crossed.

                 ----Lower bounds----  ----Upper bounds-----
  Analysis  N    Z   Nominal p Spend+  Z   Nominal p Spend++
         1  56 -0.40    0.3454 0.0297 2.79    0.0026  0.0026
         2 112  0.67    0.7473 0.0578 2.29    0.0110  0.0099
         3 167  1.68    0.9535 0.1125 1.68    0.0465  0.0375
     Total                     0.2000                 0.0500 
+ lower bound beta spending (under H1):
 Hwang-Shih-DeCani spending function with gamma = -2.
++ alpha spending:
 Hwang-Shih-DeCani spending function with gamma = -4.

Boundary crossing probabilities and expected sample size
assume any cross stops the trial

Upper boundary (power or Type I Error)
          Analysis
   Theta      1      2      3  Total  E{N}
  0.0000 0.0026 0.0099 0.0343 0.0468 104.6
  0.1995 0.0958 0.3382 0.3660 0.8000 130.9

Lower boundary (futility or Type II Error)
          Analysis
   Theta      1      2      3  Total
  0.0000 0.3454 0.4138 0.1940 0.9532
  0.1995 0.0297 0.0578 0.1125 0.2000
             T        n    Events HR futility HR efficacy
IA 1  1.523440 154.7076  55.61556       1.113       0.473
IA 2  2.372963 203.1029 111.23112       0.875       0.631
Final 4.000000 203.1029 166.84668       0.741       0.741
Accrual rates:
    Stratum 1 Stratum 2 Stratum 3
0-2     40.62     40.62     20.31
Control event rates (H1):
      Stratum 1 Stratum 2 Stratum 3
0-Inf         1       0.8       0.5
Censoring rates:
      Stratum 1 Stratum 2 Stratum 3
0-Inf         0         0         0

3.5.4 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
)
x
Time to event group sequential design with HR= 0.6 
Non-inferiority design with null HR= 1.1 
Randomization (Exp/Control):  ratio= 2 
Asymmetric two-sided group sequential design with
90 % power and 2.5 % Type I Error.
Upper bound spending computations assume
trial continues if lower bound is crossed.

                 ----Lower bounds----  ----Upper bounds-----
  Analysis  N    Z   Nominal p Spend+  Z   Nominal p Spend++
         1  45 -0.24    0.4056 0.0148 3.01    0.0013  0.0013
         2  90  0.94    0.8267 0.0289 2.55    0.0054  0.0049
         3 134  2.00    0.9772 0.0563 2.00    0.0228  0.0188
     Total                     0.1000                 0.0250 
+ lower bound beta spending (under H1):
 Hwang-Shih-DeCani spending function with gamma = -2.
++ alpha spending:
 Hwang-Shih-DeCani spending function with gamma = -4.

Boundary crossing probabilities and expected sample size
assume any cross stops the trial

Upper boundary (power or Type I Error)
          Analysis
   Theta      1      2      3  Total E{N}
  0.0000 0.0013 0.0049 0.0171 0.0233 78.2
  0.2897 0.1412 0.4403 0.3185 0.9000 99.1

Lower boundary (futility or Type II Error)
          Analysis
   Theta      1      2      3  Total
  0.0000 0.4056 0.4290 0.1420 0.9767
  0.2897 0.0148 0.0289 0.0563 0.1000
             T        n    Events HR futility HR efficacy
IA 1  14.30629 132.9678  44.64569       1.187       0.423
IA 2  21.89133 207.5835  89.29137       0.891       0.621
Final 30.00000 228.3269 133.93706       0.763       0.763
Accrual rates:
     Stratum 1 Stratum 2 Stratum 3
0-3       1.04      2.07      4.14
3-24      1.55      3.11      5.18
Control event rates (H1):
      Stratum 1 Stratum 2 Stratum 3
0-3        0.23      0.12      0.08
3-9        0.17      0.09      0.06
9-Inf      0.14      0.07      0.05
Censoring rates:
      Stratum 1 Stratum 2 Stratum 3
0-3           0         0         0
3-9           0         0         0
9-Inf         0         0         0

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.66792 41.33584  70.96406
[2,] 32.44934 64.89867 110.23544
[3,] 35.72462 71.44923 121.15304
x$eDC + x$eDE
         [,1]     [,2]     [,3]
[1,] 11.39809 14.59480 18.65280
[2,] 21.48863 29.53394 38.26873
[3,] 29.89590 44.55811 59.48305

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
)
x
Time to event group sequential design with HR= 0.6 
Equal randomization:          ratio=1
Asymmetric two-sided group sequential design with
90 % power and 2.5 % Type I Error.
Upper bound spending computations assume
trial continues if lower bound is crossed.

                 ----Lower bounds----  ----Upper bounds-----
  Analysis  N    Z   Nominal p Spend+  Z   Nominal p Spend++
         1  58 -0.24    0.4056 0.0148 3.01    0.0013  0.0013
         2 115  0.94    0.8267 0.0289 2.55    0.0054  0.0049
         3 173  2.00    0.9772 0.0563 2.00    0.0228  0.0188
     Total                     0.1000                 0.0250 
+ lower bound beta spending (under H1):
 Hwang-Shih-DeCani spending function with gamma = -2.
++ alpha spending:
 Hwang-Shih-DeCani spending function with gamma = -4.

Boundary crossing probabilities and expected sample size
assume any cross stops the trial

Upper boundary (power or Type I Error)
          Analysis
   Theta      1      2      3  Total  E{N}
  0.0000 0.0013 0.0049 0.0171 0.0233 100.5
  0.2556 0.1412 0.4403 0.3185 0.9000 127.3

Lower boundary (futility or Type II Error)
          Analysis
   Theta      1      2      3  Total
  0.0000 0.4056 0.4290 0.1420 0.9767
  0.2556 0.0148 0.0289 0.0563 0.1000
             T        n    Events HR futility HR efficacy
IA 1  10.87225 191.5728  57.37729       1.065       0.452
IA 2  16.39266 296.4606 114.75457       0.839       0.622
Final 22.95897 307.2204 172.13186       0.737       0.737
Accrual rates:
        Stratum 1 Stratum 2 Stratum 3
0-3             2         4         8
3-16.96         3         6        10
Control event rates (H1):
      Stratum 1 Stratum 2 Stratum 3
0-3        0.23      0.12      0.08
3-9        0.17      0.09      0.06
9-Inf      0.14      0.07      0.05
Censoring rates:
      Stratum 1 Stratum 2 Stratum 3
0-3           0         0         0
3-9           0         0         0
9-Inf         0         0         0