Vaccine efficacy trial design

Exact binomial approach

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 labeled N_E, P_E. Vaccine efficacy is defined as

π = 1 − P_E/P_C.

The parameter π is often labeled 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} $$

As noted, this approximation is dependent on a large sample size and small probability of events. The above 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 π₁ = 0.7 and experimental:control randomization ratio 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.4736842

We use the inversion formula to revert this to π₁ = 0.7

1 - 1 / (ratio * (1 / p1 - 1))

#> [1] 0.7

Letting the null hypothesis vaccine efficacy be π₀ = 0.3, our exact binomial null hypothesis probability that an event is in the experimental group is

pi0 <- .3
p0 <- ratio / (ratio + 1 / (1 - pi0))
p0

#> [1] 0.6774194

We also translate several vaccine efficacy values to proportion of events in the experimental group:

ve <- c(.5, .6, .65, .7, .75, .8)
prob_experimental <- ratio / (ratio + 1 / (1 - ve))
tibble(VE = ve, "P(Experimental)" = prob_experimental) %>%
  gt() %>%
  tab_options(data_row.padding = px(1)) %>%
  fmt_number(columns = 2, decimals = 3)

Chapter 12 of Jennison and Turnbull (2000) walks through how to design and analyze such a study using a fixed or group sequential design. The time-to-event approximation provides an initial approximation to computing bounds; more importantly, it provides sample size and study duration approximations that are not given by the Jennison and Turnbull approach.

The time-to-event design

Now we generate the design. If resulting alpha and beta do not satisfy requirements, adjust parameters above until a satisfactory result is obtained.

# Derive Group Sequential Design
# This determines final sample size
x <- gsSurv(
  k = k, test.type = test.type, 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
)

Now we convert this to a design with integer event counts at analyses. This is achieved by rounding interim analysis event counts from the above design and rounding up the final analysis event count. This will result in a slight change in event fractions at interim analyses as well as a slight change from the targeted 90% power. We now explain the rationale behind the spending function choices. Recall that the hazard ratio (HR) is 1 minus the VE. The ~HR at bound represents the approximate hazard ratio required to cross a bound. Thus, small HR’s at the interim analyses along with small cumulative α-spending suggest crossing an interim efficacy bound would provide a result strong enough to potentially justify the new treatment. The hazard ratio of ~0.69 (VE ~ 0.31) for the interim 1 futility bound mean that the efficacy trend would be essentially no better than the null hypothesis if the futility bound were crossed. The second analysis futility bound with approximate VE of 0.5 would be worth discussion with a data monitoring committee as well as other planners for the trial; a custom spending function could be used to set both the first and second interim bounds to desired levels.

xx <- toInteger(x)
gsBoundSummary(xx,
  tdigits = 1, logdelta = TRUE, deltaname = "HR", Nname = "Events",
  exclude = c("B-value", "CP", "CP H1", "PP")
) %>%
  gt() %>%
  tab_header(
    title = "Initial group sequential approximation",
    subtitle = "Integer event counts at analyses"
  ) %>%
  tab_options(data_row.padding = px(1))

A textual summary for the design is:

cat(summary(xx, timeunit = "months"))

Asymmetric two-sided group sequential design with non-binding futility bound, 3 analyses, time-to-event outcome with sample size 3606 and 68 events required, 90 percent power, 2.5 percent (1-sided) Type I error to detect a hazard ratio of 0.3 with a null hypothesis hazard ratio of 0.7. Enrollment and total study durations are assumed to be 8 and 24.2 months, respectively. Efficacy bounds derived using a Hwang-Shih-DeCani spending function with gamma = -3. Futility bounds derived using a Hwang-Shih-DeCani spending function with gamma = -3.

Combined summary table

We produce a summary table. Code for this is provided in the R markdown for this vignette provided with the package. This combines information from the time-to-event design for calendar timing of analyses (Time) and expected sample size at each analysis (N) along with bounds and operating characteristics for the design.

The initial approximation of bounds for the exact binomial design was generated from the time-to-event design as follows. First, we computed nominal p-value 1-sided bounds under the null hypothesis for the efficacy bounds using the normal approximation that the time-to-event design used:

efficacyNominalPValue <- pnorm(-xx$upper$bound)
efficacyNominalPValue

#> [1] 0.003610924 0.007155713 0.021200261

Then we took the inverse binomial distribution for these p-values assuming the targeted total number of cases to obtain:

qbinom(p = efficacyNominalPValue, size = xx$n.I, prob = p0) - 1

#> [1] 12 23 37

This is actually the same as the final bounds computed above:

xb$lower$bound

#> [1] 12 23 37

This and the initial approximation for the futility bound are returned from toBinomialExact():

xb$init_approx$a

#> [1] 12 23 37

xb$init_approx$b

#> [1] 20 29 38

For the futility bound, only a slight adjustment was required for the final bound:

xb$upper$bound

#> [1] 21 30 38

The vaccine efficacy at bounds should be checked to see if the evidence is convincing enough to be accepted as a clinically relevant benefit in addition to statistical benefit (efficacy bounds) or a less than relevant benefit for futility bounds.

Checking design properties

Next, we look at α- and β-spending for the time-to-event design to compare to the exact bound restricted by the discrete possible counts at bounds. For both α- and β-spending, the exact binomial design has no more spending at each analysis than allowed by each spending function. We also confirm that changing any of the exact bounds by 1 exceeds the targeting spending.

# Exact design cumulative alpha-spending at efficacy bounds
# (non-binding)
nb <- gsBinomialExact(k = xb$k, theta = xb$theta, n.I = xb$n.I, b= xb$n.I + 1, a = xb$lower$bound)
cumsum(nb$lower$prob[,1])

#> Analysis  1 Analysis  2 Analysis  3 
#> 0.001619438 0.006447739 0.017397214

# Targeted alpha-spending
xx$upper$sf(alpha, t = xx$timing, xx$upper$param)$spend

#> [1] 0.003610924 0.009107476 0.025000000

# Check that increasing any bound goes above cumulative spend
excess_alpha_spend <- matrix(0, nrow = nb$k, ncol=nb$k)
for(i in 1:xb$k){
  a <- xb$lower$bound
  a[i] <- a[i] + 1
  excess_alpha_spend[i,] <-
    cumsum(gsBinomialExact(k = xb$k, theta = xb$theta, n.I = xb$n.I, b= xb$n.I + 1, a = a)$lower$prob[,1])
}
excess_alpha_spend

#>             [,1]        [,2]       [,3]
#> [1,] 0.004979222 0.008826964 0.01933025
#> [2,] 0.001619438 0.013174419 0.02171228
#> [3,] 0.001619438 0.006447739 0.02856667

# Cumulative beta-spending for exact design
cumsum(xb$upper$prob[,2])

#> Analysis  1 Analysis  2 Analysis  3 
#>  0.01033516  0.02225609  0.09941943

# Targeted beta-spending
xx$lower$sf(beta, t = xx$timing, xx$lower$param)$spend

#> [1] 0.0144437 0.0364299 0.1000000

# Check that increasing any bound goes above cumulative spend
excess_beta_spend <- matrix(0, nrow = nb$k - 1, ncol=nb$k)
for(i in 1:(xb$k - 1)){
  b <- xb$upper$bound
  b[i] <- b[i] - 1
  excess_beta_spend[i,] <-
    cumsum(as.numeric(gsBinomialExact(k = xb$k, theta = xb$theta, n.I = xb$n.I, b = b, a = xb$lower$bound)$upper$prob[,2]))
}
excess_beta_spend

#>            [,1]       [,2]      [,3]
#> [1,] 0.02618462 0.03428681 0.1059802
#> [2,] 0.01033516 0.03746257 0.1034467

ebUpdate <- toBinomialExact(xx, observedEvents = c(20, 78))