R code to compute coverage of 95% CI using PO framework

compute-coverage.R

# load packages
library(tidyverse)

# create data frame of the potential outcomes
po <- tribble(
  ~Name,        ~Y1, ~Y0,
  "Alex Smith",     7,   5,
  "Jamie Doe",      2,   3,
  "Pat Johnson",    7,   7,
  "Jordan Lee",     5,   4,
  "Taylor Green",   3,   6,
  "Morgan White",   4,   4,
  "Casey Brown",    5,   3,
  "Drew Wilson",    3,   3,
  "Chris Bailey",   4,   2,
  "Sam Rivera",     1,   3,
  "Jesse Kim",      7,   1,
  "Robin Parker",   5,   3
) 

captured <- logical(1000)
for (i in 1:1000) {
  
  # random assignment
  zeros_and_ones <- rep(0:1, length.out = nrow(po))
  W <- sample(zeros_and_ones)
  
  # create (fake) observed data set
  observed <- data.frame(
    W_num = W, 
    W_fct = ifelse(W == 1, "Treatment", "Control"),
    Y = po$Y1*W + po$Y0*(1 - W)
  )
  
  # estimate variance
  treat <- filter(observed, W_num == 1)   # separate data frame for treatment group
  control <- filter(observed, W_num == 0) # separate data frame for control group
  var(treat$Y)/6 + var(control$Y)/6  # formula to estimate variance
  
  # 95% CI: estimate +/- 1.96*SE
  estimate <- mean(treat$Y) - mean(control$Y)
  se <- sqrt(var(treat$Y)/6 + var(control$Y)/6) # SE = sqrt(variance)
  lwr95 <- estimate - 1.96*se
  upr95 <- estimate + 1.96*se
  
  # did it work?
  truth <- mean(po$Y1) - mean(po$Y0)
  captured[i] <- (lwr95 < truth) & (upr95 > truth)
  
}

table(captured)