library(gsDesign)
library(gt)
library(dplyr)This article explores a method of approximating a design using the exact binomial method of Chan and Bohidar (1998) by a time-to-event design using the method of Lachin and Foulkes (1986). This allows use of spending functions to derive boundaries for the exact method. The time-to-event design can not only be used to set approximate boundaries for the Chan and Bohidar (1998) method, but to allow specification of enrollment duration and study duration to determine enrollment rates and sample size required. This also allows us to illustrate the concept of super-superiority often used in prevention studies.
We begin with the assumption that we will a require a large sample size due to a endpoint with a small incidence rate. This could apply to a vaccine study or other prevention study with a relatively small number of events expected.
Paralleling the notation of Chan and Bohidar (1998), we assume \(N_C, P_C\) to be binomial sample size and probability of an event for each participant assigned to control; for the experimental treatment group, these are labelled \(N_E, P_E\). Vaccine efficacy is defined as
\[\pi = 1 - P_E/P_C.\]
The parameter \(\pi\) is also often labelled as \(VE\) for vaccine efficacy. Taking into account the randomization ratio \(r\) (experimental / control) the approximate probability that any given event is in the experimental group is
\[ \begin{aligned} p &= rP_E/(rP_E+ P_C)\\ &= r/(r + P_C/P_E)\\ &= r/(r + (1-\pi)^{-1}). \end{aligned} \]
This can be inverted to obtain
\[\pi = 1- \frac{1}{r(1/p-1)}. \]
For our example of interest, we begin with an alternate hypothesis vaccine efficacy \(\pi_1 = 0.7\) and \(r=3\). This converts to an alternate hypothesis approximate probability that any event is in the experimental group of
pi1 <- .7
ratio <- 3
p1 <- ratio / (ratio + 1 / (1 - pi1))
p1#> [1] 0.4736842We use the inversion formula to revert this to \(\pi_1 = 0.7\)
1 - 1 / (ratio * (1 / p1 - 1))#> [1] 0.7Letting the null hypothesis vaccine efficacy be \(\pi_0 = 0.3\), our exact binomial null hypothesis probability is
pi0 <- .3
p0 <- ratio / (ratio + 1 / (1 - pi0))
p0#> [1] 0.6774194Chapter 12 of Jennison and Turnbull (2000) walks through how to design and analyze such a study using a fixed or group sequential design. This may have advantages and disadvantages compared to what is proposed here. However, the time-to-event approximation gives not only proposed bounds, but also sample size and study duration approximations.
For a time-to-event formulation with exponential failure rates \(\lambda_C\) for control and \(\lambda_E\) for experimental group assigned participants, we would define
\[\pi = 1 - \lambda_E / \lambda_E\]
which is 1 minus the hazard ratio often used in time-to-event studies. In the following we examine how closely the time-to-event method using asymptotic distributional assumptions can approximate an appropriate exact binomial design.
We will also define a planned number of events at each of \(K\) planned analyses by \(D_k, 1\le k\le K\).
We begin by specifying parameters. The alpha and beta parameters will not be met exactly due to the discrete group sequential probability calculations performed. Thus, you may need to adjust these parameters slightly to ensure your final design operating characteristics are within the targeted range. The current version includes only designs that use non-binding futility bounds. The design is generated by first using asymptotic theory for a time-to-event design with specified spending functions. This design is then adapted to a design using the exact binomial method of Chan and Bohidar (1998). The randomization ratio (experiemental/control) was assumed to be 3:1 as in the Logunov et al. (2021) trial.
alpha <- 0.023 # Type I error; this was adjusted from .025 to ensure Type I error control
beta <- 0.09 # Type II error (1 - power); this was reduced from .1 to .09 to ensure power
k <- 3 # number of analyses in group sequential design
timing <- c(.5, .8) # Relative timing of interim analyses compared to final
sfu <- sfLDOF # Efficacy bound spending function (O'Brien-Fleming-like here)
sfupar <- 0 # Parameter for efficacy spending function
sfl <- sfHSD # Futility bound spending function (Hwang-Shih-DeCani here)
sflpar <- -12 # Futility bound spending function parameter (conservative)
timename <- "Month" # Time unit
failRate <- .002 # Exponential failure rate
dropoutRate <- .0001 # Exponential dropout rate
enrollDuration <- 4 # Enrollment duration
trialDuration <- 8 # Planned trial duration
VE1 <- .7 # Alternate hypothesis vaccine efficacy
VE0 <- .3 # Null hypothesis vaccine efficacy
ratio <- 3 # Experimental/Control enrollment ratioNow we generate the design. If resulting alpha and beta do not satisfy requirements, adjust parameters above until a satisfactory result is obtained. Press the code button below to review detailed code, which should not require user alteration.
# Derive Group Sequential Design
# This determines final sample size
x <- gsSurv(
k = k, test.type = 4, alpha = alpha, beta = beta, timing = timing,
sfu = sfu, sfupar = sfupar, sfl = sfl, sflpar = sflpar,
lambdaC = failRate, eta = dropoutRate,
# Translate vaccine efficacy to HR
hr = 1 - VE1, hr0 = 1 - VE0,
R = enrollDuration, T = trialDuration,
minfup = trialDuration - enrollDuration, ratio = ratio
)
gsBoundSummary(x, tdigits = 1) %>%
gt() %>%
tab_header(title = "Initial group sequential approximation")| Initial group sequential approximation | |||
|---|---|---|---|
| Analysis | Value | Efficacy | Futility |
| IA 1: 50% | Z | 3.0105 | -1.1284 |
| N: 12111 | p (1-sided) | 0.0013 | 0.8704 |
| Events: 35 | ~HR at bound | 0.2138 | 1.0919 |
| Month: 5 | P(Cross) if HR=0.7 | 0.0013 | 0.1296 |
| P(Cross) if HR=0.3 | 0.2653 | 0.0002 | |
| IA 2: 80% | Z | 2.3042 | 0.6115 |
| N: 12111 | p (1-sided) | 0.0106 | 0.2704 |
| Events: 55 | ~HR at bound | 0.3415 | 0.5786 |
| Month: 6.8 | P(Cross) if HR=0.7 | 0.0110 | 0.7298 |
| P(Cross) if HR=0.3 | 0.7632 | 0.0082 | |
| Final | Z | 2.0610 | 2.0610 |
| N: 12111 | p (1-sided) | 0.0196 | 0.0196 |
| Events: 69 | ~HR at bound | 0.3942 | 0.3942 |
| Month: 8 | P(Cross) if HR=0.7 | 0.0230 | 0.9770 |
| P(Cross) if HR=0.3 | 0.9100 | 0.0900 | |
Now we change the event counts at each analysis by rounding down to move away from the continuous, unrounded event counts above. There are small changes from the above design.
# Round up event counts and update spending based on slightly updated timing and then re-derive bounds.
# This will then be used to set event counts and bounds for the exact binomial bounds
counts <- round(x$n.I) # Round for interim counts
counts[k] <- ceiling(x$n.I[k]) # Round up for final count
timing <- counts / max(counts)
xx <- gsDesign(
k = k, test.type = 4, n.I = counts, maxn.IPlan = x$n.I[k],
alpha = alpha, beta = beta,
delta = x$delta, delta1 = x$delta1, delta0 = x$delta0,
sfu = sfu, sfupar = sfupar, sfl = sfl, sflpar = sflpar,
usTime = timing,
lsTime = timing
)
zupper <- xx$upper$bound # Updated upper bounds
zlower <- xx$lower$bound # Updated lower bounds
# For non-inferiority and super-superiority trials
xx$hr0 <- x$hr0xxsum <- gsBoundSummary(xx, deltaname = "HR", logdelta = TRUE, Nname = "Events")
xxsum %>%
gt() %>%
tab_header(title = "Updated design with spending based on integer event counts")| Updated design with spending based on integer event counts | |||
|---|---|---|---|
| Analysis | Value | Efficacy | Futility |
| IA 1: 49% | Z | 3.0355 | -1.1640 |
| Events: 34 | p (1-sided) | 0.0012 | 0.8778 |
| ~HR at bound | 0.2471 | 1.0435 | |
| P(Cross) if HR=0.7 | 0.0012 | 0.1222 | |
| P(Cross) if HR=0.3 | 0.2532 | 0.0002 | |
| IA 2: 80% | Z | 2.3082 | 0.5995 |
| Events: 55 | p (1-sided) | 0.0105 | 0.2744 |
| ~HR at bound | 0.3756 | 0.5955 | |
| P(Cross) if HR=0.7 | 0.0109 | 0.7258 | |
| P(Cross) if HR=0.3 | 0.7621 | 0.0079 | |
| Final | Z | 2.0601 | 2.0601 |
| Events: 69 | p (1-sided) | 0.0197 | 0.0197 |
| ~HR at bound | 0.4263 | 0.4263 | |
| P(Cross) if HR=0.7 | 0.0230 | 0.9770 | |
| P(Cross) if HR=0.3 | 0.9112 | 0.0888 | |
Now we invert the upper and lower bounds based on the inverse binomial distribution to get \(a_k<b_k, 1\le k\le K\) where \(N_{E,k}\le a_k\) results in a positive finding in favor of the experimental arm and \(N_{E,k}\ge b_k\) results in futility, \(1\le k\le K\):
# Translate vaccine efficacy to exact binomial probabilities
p0 <- (1 - VE0) * ratio / (1 + (1 - VE0) * ratio)
p1 <- (1 - VE1) * ratio / (1 + (1 - VE1) * ratio)
# Lower bound probabilities are for efficacy and Type I error should be controlled under p0
a <- qbinom(p = pnorm(-zupper), size = counts, prob = p0)
a[k] <- a[k] - 1
# Upper bound probabilities are for futility and spending should be under p1
b <- qbinom(p = pnorm(zlower), size = counts, prob = p0, lower.tail = FALSE)
# a < b required for each analysis; subtracting 1 makes one bound or the other
# cross at the end
xxx <- gsBinomialExact(k = k, theta = c(p0, p1), n.I = counts, a = a, b = b)
xxx#> Bounds
#> Analysis N a b
#> 1 34 14 26
#> 2 55 29 35
#> 3 69 38 39
#>
#> Boundary crossing probabilities and expected sample size assume
#> any cross stops the trial
#>
#> Upper boundary
#> Analysis
#> Theta 1 2 3 Total E{N}
#> 0.6774 0.1838 0.6053 0.1869 0.9760 53.9
#> 0.4737 0.0005 0.0107 0.0624 0.0737 51.2
#>
#> Lower boundary
#> Analysis
#> Theta 1 2 3 Total
#> 0.6774 0.0013 0.0134 0.0094 0.0240
#> 0.4737 0.2918 0.5332 0.1013 0.9263Now we make comparisons between the asymptotic, time-to-event design in xx to the exact binomial design in xxx. We begin with comparison of nominal one-sided p-values for each of the bounds.
# Nominal p-values at efficacy bounds
pnorm(-xx$upper$bound) # Asymptotic#> [1] 0.001200857 0.010493965 0.019695361pbinom(xxx$lower$bound, prob = p0, size = xxx$n.I) # Exact binomial#> [1] 0.001274188 0.014367258 0.018708455# Nominal p-values for futility bounds
pnorm(-xx$lower$bound) # Asymptotic#> [1] 0.87779020 0.27443463 0.01969536pbinom(xxx$upper$bound, prob = p0, size = xxx$n.I) # Exact binomial#> [1] 0.90142900 0.30177551 0.03325157Next we examine cumulative non-binding \(\alpha\)-spending for efficacy bounds. While the cumulative interim spending is slightly higher than planned, the cumulative final spending is controlled at the targeted 0.025. If it is desired to lower interim spending to the target, slight adjustments could be made; for the purpose of this example we do not dive into this further.
# Cumulative non-binding alpha-spending at lower bounds
cumsum(xx$upper$spend) # Asymptotic design#> [1] 0.001200857 0.010884213 0.023000000# Exact binomial design
# Compute by removing futility bounds
balt <- xxx$n.I + 1
xxxalt <- gsBinomialExact(k = k, theta = c(p0, p1), n.I = counts, a = a, b = balt)
cumsum(xxxalt$lower$prob[, 1])#> Analysis 1 Analysis 2 Analysis 3
#> 0.001274188 0.014663600 0.024088556Next we look at \(\beta\)-spending for the two bounds. The exact bounds control \(\beta\)-spending close to the asymptotic target at interims and fully controls the complete \(\beta=0.10\) at the final analysis, ensuring power is achieved with the exact design.
# Cumulative non-binding beta-spending at lower bounds
cumsum(xx$lower$spend) # Asymptotic design#> [1] 0.0002039565 0.0078850068 0.0900000000cumsum(xxx$upper$prob[, 2]) # Exact binomial design#> Analysis 1 Analysis 2 Analysis 3
#> 0.000523268 0.011250067 0.073694635Now we examine the approximate vaccine efficacy at efficacy and futility bounds for the asymptotic design vs. the exact vaccine efficacy for the exact binomial design. We begin using the inversion formula from above to convert the proportion of events in the experimental group at the bounds to vaccine efficacy at the bounds.
p_lower <- xxx$lower$bound / xxx$n.I
p_upper <- xxx$upper$bound / xxx$n.I
p_lower#> [1] 0.4117647 0.5272727 0.5507246p_upper#> [1] 0.7647059 0.6363636 0.5652174Now we convert to observed vaccine efficacy at each bound
1 - 1 / (ratio * (1 / p_lower - 1)) # Efficacy bound#> [1] 0.7666667 0.6282051 0.59139781 - 1 / (ratio * (1 / p_upper - 1)) # Futility bound#> [1] -0.08333333 0.41666667 0.56666667This is approximated reasonably by the approximations using the time-to-event design:
# Translate ~HR at efficacy bound to VE for asymptotic design
1 - xxsum[c(3, 8, 13), 3] # Efficacy bound#> [1] 0.7529 0.6244 0.57371 - xxsum[c(3, 8, 13), 4] # Futility bound#> [1] -0.0435 0.4045 0.5737We have provided an extended example to show that a Chan and Bohidar (1998) exact binomial using spending function bounds can be derived in a two-step process that delivers sample size and bounds by 1) deriving a related time-to-event design using asymptotic methods and then 2) converting to an exact binomial design. Small adjustments were made to target Type I and Type II error probabilities in the asymptotic approximation to ensure the exact binomial Type I and Type II error rates were achieve. The method seems a reasonable and straightforward approach to develop a complete design.