Why can’t we just peek at our data with α = 0.05 every time? 🤔 Gotta see it to believe it! Simulating the inflated error, seeing why it’s not Bonferroni, and rebuilding O’Brien-Fleming & Lan-DeMets boundaries to match gsDesign. Hope it’s not an illusion of learning 😅

Motivations

I’ve always heard “interim analysis” thrown around in trial papers but never really understood what’s happening under the hood. Why can’t we just peek at 1.96 every time? And isn’t this just Bonferroni? Turns out it’s not! Let’s dig in and simulate our way to understanding.

Objectives:

• What Is Interim Analysis?
• What’s O’Brien-Fleming Method?
• How To Perform It?
• How To Do It From Scratch?
  • O’Brien-Fleming Method
  • Lan De Mat & Obrien Fleming Method
  • Comparison
• Opportunities For Improvement
• Lessons Learnt

What Is Interim Analysis?

An interim analysis is a planned look at the data half way through a study, before all participants have been enrolled or followed up. They do it for practical and ethical reasons: if a treatment is clearly working, or clearly harmful, it isn’t right to keep randomizing people to the worse option, and if it’s clearly going nowhere, they can save time and resources by stopping early. The catch is that every extra look is another chance to find a “significant” result purely by luck. Peeking twice at the usual 5% threshold pushes the false-positive rate up. Interim analysis methods, such as O’Brien-Fleming boundaries, fix this by setting a stricter bar at the early looks so that the overall error rate stays at 5%.

Now you know me, I’m the see it to believe kind of person. Let’s try to simulate this with multiple trials and see if that’s true

library(tidyverse)
n <- 500
prob <- 0.5
peek <- c(n/2,n)
zscore_vec <- c()
trial <- 100000
count_false <- vector(mode = "numeric", length = trial)

for (trial_i in 1:trial) {
  # for reproducibility
  set.seed(trial_i)
  
  # get sim data
  control <- rbinom(n,1,prob)
  rx <- rbinom(n,1,prob)
  
  # look at when n is 250, and when n is 500
  for (i in 1:length(peek)) {
    m  <- peek[i]
    pA <- mean(rx[1:m])
    pB <- mean(control[1:m])
    p_pool <- (pA + pB) / 2
    se <- sqrt(p_pool * (1 - p_pool) * (2 / m))
    
    # calculate z score
    zscore_vec[i] <- abs(pA - pB) / se
  }
  # record zscore if > 1.96 (~pval 0.05)
  if (any(zscore_vec > 1.96)) { count_false[trial_i] <- 1 }
}

# what percent of family-wise error?
(val <- count_false |> mean())

## [1] 0.08952

Lots of code, I hope some of the comments on the codes are self-explanatory. So we basically simulated a no difference result of 2 groups. We used zscore of 1.96 which is our 5% alpha to see if any of these z scores were above the threshold (it’s not supposed to, remember, there is no difference in this truth!). We expect it to be ~5% but we got 8.952% instead! Wow! That’s 3.952% higher than our threshold, not good! 🤔

Note to self: The confusion of false postive/discovery/family-wise examples. If you test 90 placebo and 10 real rx and get 15 “significants”, 5 of them wrong: the false positive rate is 5/90 (how often placebos get credit) aka (1 – specificity), the false discovery rate is 5/15 (how many of your hits are wrong) aka (1 – ppv [aka precision]), and the family-wise error rate is the chance of getting even one false hit across the whole set, which is what interim analysis controls across looks.

wow, doesn’t that remind us of something here? Multiplicity, yea? Maybe, kind of. Now, what if we’re being conservative like bonferonni and split that 5% twice? and see how much false positive did we get?

count_false <- vector(mode = "numeric", length = trial)

for (trial_i in 1:trial) {
  # for reproducibility
  set.seed(trial_i)
  
  # get sim data
  control <- rbinom(n,1,prob)
  rx <- rbinom(n,1,prob)
  
  # look at when n is 250, and when n is 500
  for (i in 1:length(peek)) {
    m  <- peek[i]
    pA <- mean(rx[1:m])
    pB <- mean(control[1:m])
    p_pool <- (pA + pB) / 2
    se <- sqrt(p_pool * (1 - p_pool) * (2 / m))
    
    # calculate z score
    zscore_vec[i] <- abs(pA - pB) / se
  }
  # record zscore if > 2.24, two tails with two looks (divide 5% to 2)
  if (any(zscore_vec > 2.24)) { count_false[trial_i] <- 1 }
}

# what percent of family wise error rate?
(val <- count_false |> mean())

## [1] 0.04093

What !?!?! 4.093% !?! It’s a bit lower than expected and we may be wasting 0.907% of our threshold! This is also not good! If there truly is a signal, we may be missing it because we are being too conservative. So, what if we use the O’Brien-Fleming method to set our threshold? Let’s see how that works out.

What’s O’Brien-Fleming Method?

O’Brien-Fleming is a way of setting stopping boundaries for interim analyses, and it’s the most widely used method in clinical trials. The idea is to be very strict early and lenient late. At the first look, when you have little data, you need overwhelming evidence to stop, a z-score of about 2.80 (pval = 0.005110261) with two looks. By the final analysis, the bar drops back to about 1.98 (pval = 0.04770353), barely above the usual 1.96. Because so little alpha is spent early, the trial almost never stops on a fluke, and the final analysis keeps nearly all of its power. 🙌 Let’s take a look of its family wise error rate like before

count_false <- vector(mode = "numeric", length = trial)

for (trial_i in 1:trial) {
  # for reproducibility
  set.seed(trial_i)
  
  # get sim data
  control <- rbinom(n,1,prob)
  rx <- rbinom(n,1,prob)
  
  # look at when n is 250, and when n is 500
  for (i in 1:length(peek)) {
    m  <- peek[i]
    pA <- mean(rx[1:m])
    pB <- mean(control[1:m])
    p_pool <- (pA + pB) / 2
    se <- sqrt(p_pool * (1 - p_pool) * (2 / m))
    
    # calculate z score
    zscore_vec[i] <- abs(pA - pB) / se
  }
  # record zscore if > O'Brien-Fleming, two looks, notice the zscores are no longer 1.96
  bounds <- c(2.80, 1.98) 
  if (any(zscore_vec > bounds)) {count_false[trial_i] <- 1}
}

# what percent of family wise error?
(val <- count_false |> mean())

## [1] 0.0489

Alright! That’s much better! 👍 Closer to the 5%. As I was learning this, I just realized that this is very different from multiplicity because the second look technically is correlated to the first look because it looks at the same data again! The two looks share half their data, so their z-scores have a correlation of √0.5 ≈ 0.71. Bonferroni treats them as if they were barely related, so it over-adjusts. 🤔 Wow, that totally makes sense now! Ok, now how do we get sample calculation + obrien fleming method with gsDesign? Take note below, the history of obrien-fleming and Lan De Mat

How To Perform It?

Let’s say we want to design a RCT of cefdinir vs cephalexin. Looking at this study, without adjustment, cefdinir had a 23.4% treatment failure compared to cephalexin 12.5%. With 80% power, what sample size, what’s the upper boundary of z score we need if we want 3 looks on our interim analyses?

library(gsDesign)
p1 <- 0.234
p2 <- 0.125

# fixed design: total sample size, both groups together
(n_fix <- nBinomial(p1 = p1, p2 = p2, alpha = 0.025, beta = 0.2))

## [1] 386.8202

Total sample size needed is 386.820172, will have to divide this by 2. Let’s have gsDesign tell us the upper boundary zscore and when to look.

# group sequential: 3 equally spaced looks, O'Brien-Fleming spending
d <- gsDesign(k = 3, test.type = 2, alpha = 0.025, beta = 0.2,
              n.fix = n_fix, sfu = sfLDOF)

sample <- d$n.I / 2           
bound <- d$upper$bound           
pval <- 2 * (1 - pnorm(d$upper$bound))  
(table <- tibble(number_per_group=sample,zscore_needed=bound,p_value=pval))

## # A tibble: 3 × 3
##   number_per_group zscore_needed  p_value
##              <dbl>         <dbl>    <dbl>
## 1             65.3          3.71 0.000207
## 2            131.           2.51 0.0120  
## 3            196.           1.99 0.0463

Per group, we will need the above numbers to take a peek. Hooray! Now we kind of understand the surface of the goal and why we need to do this so we don’t inflate our family-wise error rate. Now, let’s try to do this from scratch and see if we can get the same numbers as gsDesign did. Take note that this time we will be using the first simulate DGP with no diff, not the cefdinir study.

How To Do It From Scratch?

So who came up with what? 🤔 O’Brien and Fleming came first (1979). Their idea: be super strict early when data is thin, and relax at the end. Every boundary is just one number c divided by √t, and you search for the c that keeps the family-wise error at 5%. Simple, but you had to fix when you’d look ahead of time. Then Lan and DeMets came along (1983) and said, “great idea, but what if we can’t stick to the schedule?” So they made the spending function: decide how much of the 5% you’ve used up at any point, and each look only spends its share. Same strict-early shape, but now you can peek whenever you need to. That’s what gsDesign’s sfLDOF does, and why everyone just calls it “O’Brien-Fleming” anyway! 🙌

Too many words? Yea, I know how you feel. Let’s break things down and code it to understand it.

Let’s start with the simplest O’Brien-Fleming method. We will simulate 3 looks at n=500, with 100,000 trials first and calculate the zscores for each look.

n <- 500
prob <- 0.5
trial <- 100000
k <- 3
zscore_mat <- matrix(nrow = trial, ncol = k)
gp <- round(n/k)
peek <- c(gp, gp*2, n)

for (trial_i in 1:trial) {
  # for reproducibility
  set.seed(trial_i)
  
  # get sim data
  control <- rbinom(n,1,prob)
  rx <- rbinom(n,1,prob)
  
  # peek at respective k
  for (i in 1:k) {
    m  <- peek[i]
    pA <- mean(rx[1:m])
    pB <- mean(control[1:m])
    p_pool <- (pA + pB) / 2
    se <- sqrt(p_pool * (1 - p_pool) * (2 / m))
    
    # calculate z score, record it
    zscore_mat[trial_i, i] <- abs(pA - pB) / se
  }
}

O’Brien-Fleming Method

1. Start with final z score, let’s choose 1.96 as a guess and search upwards.
2. Then back calculate second and first look by dividing by square root of information fraction (t). For example, if we have 3 looks, the first look will be 1/3 of the information, second look will be 2/3 of the information and final look will be 3/3 of the information.
3. Then calculate the total family wise error rates for each trial and call it alpha_final.
4. Then look at the simulated data and filter out alpha_final between 0.045 and 0.05, then see what z scores for first and second look.

# obrien fleming
df <- data.frame(zscore_mat) |>
  as_tibble()

final_z_seq <- seq(1.96, 2.5, 0.001)
t2_z_s_seq <- final_z_seq / sqrt(2/3)
t1_z_s_seq <- final_z_seq / sqrt(1/3)

for (i in 1:length(final_z_seq)) {

# record 1 if zscore passes the respective threshold for each look
alpha_final <- df |>
  mutate(crossed = case_when(
    X3 > final_z_seq[i] ~ 1,
    X2 > t2_z_s_seq[i] ~ 1,
    X1 > t1_z_s_seq[i] ~ 1,
    TRUE ~ 0
  )) |>
  pull(crossed) |>
  mean()

if (alpha_final < 0.05 & alpha_final >= 0.045) { print(paste0("total alpha: ", alpha_final ," .found it, it's on ",i, " ", t1_z_s_seq[i], " ", t2_z_s_seq[i], " ",final_z_seq[i])) }
}

## [1] "total alpha: 0.0494499999999998 .found it, it's on 65 3.50567083451941 2.47888361969658 2.024"
## [1] "total alpha: 0.0463299999999999 .found it, it's on 66 3.50740288532698 2.48010836456797 2.025"
## [1] "total alpha: 0.0454099999999999 .found it, it's on 67 3.50913493613455 2.48133310943936 2.026"

Nice! That worked! Alright, now let’s look at the Lan De Mat & Obrien Fleming method. The idea is to calculate the alpha spent at each look and then find the z scores for each look.

Lan De Mat & Obrien Fleming Method

So Lan De Mat introduced spend function, which makes the O’Brien-Fleming method more flexible. The idea is to calculate the alpha spent at each look and then find the z scores for each look. The formula for the spend function is 2 - 2 * pnorm(qnorm(1 - alpha/2) / sqrt(t)). here is the paper that derives it, might be a good blog post next time to learn. For now we’ll take it as is and roll with it. The idea is that it’s the chance a random walk crosses a fixed line by time t (the reflection principle)

1. Calculate spending for look 1 with the formula. (note the alpha here is 0.025 for one tail). This spending will be in percentage (like p val).
2. Sort descending of the first look zscore from our simulation data, then remove the spending percentage of the first look, but record the first zscore value (that would be our first z score on first look). Save into a new dataframe. PS: this is like saying, I’m OK with the false positive of these much on first look.
3. Calculate the spending for look 2 with the formula, the subtract the first look spending.
4. Sort descending of the second look zscore of the new dataframe, then record the second zscore value, and remove the spending percentage of the second look on the new dataframe.
5. Calculate the spending for look 3 (which is 0.025), then subtract of the third look spending to second+one (cumulative) look spending, then do the same as above, record zscore and filter out the rest based on spending 3 percentage.

Wow, ok too many words above. Use code to help understand.

# calc spend1 
t <- 1/3                                   
alpha <- 0.025                                 
spend1 <- 2 - 2 * pnorm(qnorm(1 - alpha/2) / sqrt(t))

# filter for those X1 with less than zscore1
df1 <- data.frame(zscore_mat) |>
  as_tibble() |>
  arrange(desc(X1)) |>
  slice(round(spend1*trial*2+1):trial)

zscore1 <- data.frame(zscore_mat) |>
  as_tibble() |>
  arrange(desc(X1)) |>
  slice_head(n = round(spend1*trial*2)) |>
  slice_tail(n = 1) |>
  pull(X1)

# calc spend2
t <- 2/3
spend2_cum <- 2 - 2 * pnorm(qnorm(1 - alpha/2) / sqrt(t))   
spend2 <- spend2_cum - spend1                                

# find how many to look at 
head_num <- round(spend2 * 2 * trial)                                             

# head of however many false positive allowed
zscore2 <- df1 |>
  arrange(desc(X2)) |> 
  slice_head(n = head_num) |>
  slice_tail(n = 1) |>
  pull(X2)

df2 <- df1 |>
  arrange(desc(X2)) |>
  slice((head_num+1):nrow(df1))

# calc spend3
spend3 <- alpha - spend2_cum      
head_num <- round(spend3 * 2 * trial)                

zscore3 <- df2 |>
  arrange(desc(X3)) |>
  slice(head_num) |>
  pull(X3)

zscore1 ; zscore2 ; zscore3

## [1] 3.612396

## [1] 2.480686

## [1] 2.023874

Alright, not too shabby. Looks OK. Let’s use new data to calculate our total family wise error rate.

### check family wise error
count_false <- vector(mode = "numeric", length = trial)

for (trial_i in 1:trial) {
  # for reproducibility
  set.seed(trial_i + trial)
  
  # get sim data
  control <- rbinom(n,1,prob)
  rx <- rbinom(n,1,prob)
  
  # peek at respective k
  for (i in 1:length(peek)) {
    m  <- peek[i]
    pA <- mean(rx[1:m])
    pB <- mean(control[1:m])
    p_pool <- (pA + pB) / 2
    se <- sqrt(p_pool * (1 - p_pool) * (2 / m))
    
    # calculate z score
    zscore_vec[i] <- abs(pA - pB) / se
  }
  # new bounds we found 
  bounds <- c(zscore1, zscore2, zscore3) 
  if (any(zscore_vec > bounds)) {count_false[trial_i] <- 1}
}

# what percent of family wise error?
(val <- count_false |> mean())

## [1] 0.05027

Sweet! 🙌 Alright, let’s see how the professionals do it?

## gsDesign
d0 <- gsDesign(k = 3, test.type = 2, alpha = 0.025, sfu = sfLDOF)
d0$upper$bound      

## [1] 3.710303 2.511427 1.993048

Comparison

Let’s put them all on a table

Not bad at all! 🙌 Our sims land within about 0.1 of the professionals on both methods. The small gaps are from our binary outcome ?maybe. The z-scores come in little jumps, so we can’t hit the boundaries exactly. Also, take note that we got here by brute-force simulation, not the real math. gsDesign uses the multivariate normal distribution with the correlation matrix between looks (remember sqrt(t1/t2)?) to calculate these boundaries exactly. That’s a project for next time!

Opportunities For Improvement

• Calculate the boundaries exactly with mvtnorm::pmvnorm (no simulation noise)
• Simulate power under a real effect, not just the null
• Add a futility boundary (test.type = 4) to stop when it’s clearly going nowhere
• Look at how stopping early overestimates the effect size
• Use a continuous outcome to get rid of the binary lumpiness
• Learn derivation of reflection principle behind the Lan-DeMets formula

Lessons learnt

• learnt the simple interim analysses
• learnt basic gsDesign
• learnt that alpha borrowing in interim analyses is not the same as multiplicity, because the second look is correlated to the first look. Bonferroni treats them as if they were barely related, so it over-adjusts. 🤔
• learnt interim analyses true zscore estimation requires multivariate normal distribution

If you like this article:

• please feel free to send me a comment or visit my other blogs
• please feel free to follow me on BlueSky, twitter, GitHub or Mastodon
• if you would like collaborate please feel free to contact me