library(gsDesign)
library(gt)
library(dplyr)

# Introduction

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.

# Parameterization

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.

## 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 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
#>  0.4736842

We use the inversion formula to revert this to $$\pi_1 = 0.7$$

1 - 1 / (ratio * (1 / p1 - 1))
#>  0.7

Letting 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
#>  0.6774194

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

# The time-to-event approach

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

# Generating a design

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 ratio

## 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. 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$hr0
xxsum <- 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

# Converting to an exact binomial design

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

Now 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
#>  0.001200857 0.010493965 0.019695361
pbinom(xxx$lower$bound, prob = p0, size = xxx$n.I) # Exact binomial #>  0.001274188 0.014367258 0.018708455 # Nominal p-values for futility bounds pnorm(-xx$lower$bound) # Asymptotic #>  0.87779020 0.27443463 0.01969536 pbinom(xxx$upper$bound, prob = p0, size = xxx$n.I) # Exact binomial
#>  0.90142900 0.30177551 0.03325157

Next 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
#>  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.024088556 Next 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 #>  0.0002039565 0.0078850068 0.0900000000 cumsum(xxx$upper$prob[, 2]) # Exact binomial design #> Analysis 1 Analysis 2 Analysis 3 #> 0.000523268 0.011250067 0.073694635 Now 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
#>  0.4117647 0.5272727 0.5507246
p_upper
#>  0.7647059 0.6363636 0.5652174

Now we convert to observed vaccine efficacy at each bound

1 - 1 / (ratio * (1 / p_lower - 1)) # Efficacy bound
#>  0.7666667 0.6282051 0.5913978
1 - 1 / (ratio * (1 / p_upper - 1)) # Futility bound
#>  -0.08333333  0.41666667  0.56666667

This 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
#>  0.7529 0.6244 0.5737
1 - xxsum[c(3, 8, 13), 4] # Futility bound
#>  -0.0435  0.4045  0.5737

# Summary

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