Simulating efficacy trials

library(conflicted)
library(tidyverse)

── Attaching core tidyverse packages ───────────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     

library(tictoc)

Everyone loves a good randomised controlled trial because the mean outcome of people who were exposed to treatment minus the mean outcome of people who were exposed to control – a between-participant difference – is an unbiased estimator of the the mean of within-participant individual treatment effects.

Let’s simulate some data where we know the average treatment effect (ATE). This is an efficacy trial, so we are only interested in the ATE for the study participants. There is still an inferential problem since not everyone gets treatment and not everyone gets control.

set.seed(3137942)
theN <- 500

sim_dat <- tibble(latent_magic = 50 + rnorm(theN, 0, sqrt(.6)),
                  y0 = round(latent_magic + rnorm(theN, 0, sqrt(.4)), 1),
                  y1 = round(latent_magic + 0.2 + rnorm(theN, 0, sqrt(.4)), 1),
                  ITS = y1 - y0) |>
  select(-latent_magic)

Take a look:

sim_dat

y0 is the potential outcome under control, y1 the potential outcome under treatment, and ITS is the individual treatment effect, y1 – y0.

We can get the true average treatment effect of an efficacy trial by averaging all those ITSes:

trueATE <- mean(sim_dat$ITS)
trueATE

[1] 0.2636

That’s the true answer we’re trying to recover with an RCT.

Now let’s simulate an RCT – one randomisation of participants to treat or control. Here I’m using a coin toss.

rct_dat <- sim_dat |>
  mutate(
    treat = rbinom(n(), 1, .5),
    y = treat * y1 + (1 - treat) * y0,
    y1 = ifelse(treat == 1, y1, NA),
    y0 = ifelse(treat == 0, y0, NA)
  ) |>
  select(-ITS)

rct_dat

We only see one of the potential outcomes, depending on which was realised, hence the missing data on y0 and y1 (NAs).

Estimate ATE and its 95% CI using linear regression:

lm_mod <- lm(y ~ treat, data = rct_dat)
cbind(ATE = coef(lm_mod), confint(lm_mod))["treat",]

      ATE     2.5 %    97.5 % 
0.2935164 0.1172047 0.4698281 

Repeat that 2,000 times to see the average of the ATEs across simulated trials:

sim_reses <- data.frame()

for (i in 1:2000) {
  rct_dat <- sim_dat |>
    mutate(treat = rbinom(n(), 1, .5),
           y = treat * y1 + (1 - treat) * y0)
  
  lm_mod <- lm(y ~ treat, data = rct_dat)
  this_res <- cbind(ATE = coef(lm_mod), confint(lm_mod))["treat", ]
  
  sim_reses <- bind_rows(sim_reses, this_res)
}

sim_reses  |>
  mutate(across(everything(), \(x) round(x, 2)))

The mean estimated ATE across simulated trials is really close to the true ATE…

c(
  `mean of ATE estimates` = mean(sim_reses$ATE),
  `true ATE` = trueATE,
  `true ATE - mean of estimates` = trueATE - mean(sim_reses$ATE) 
) |> round(4)

       mean of ATE estimates                     true ATE 
                      0.2605                       0.2636 
true ATE - mean of estimates 
                      0.0031 

… but there’s quite a bit of variation around the mean answer.

sim_reses |>
  ggplot(aes(x = ATE)) +
  geom_histogram(bins = 40) +
  xlab("estimated ATE")

Let’s look at the 95% confidence interval coverage:

sim_reses <- sim_reses |>
  mutate(trueATE_in_interval = trueATE >= `2.5 %` &
                               trueATE <= `97.5 %`)
cover <- mean(sim_reses$trueATE_in_interval)
cover

[1] 0.9745

97% of the simulated studies’ 95% confidence intervals included the true result.

5.79 sec elapsed