An example to compute Fisher’s exact p-value and CI

Suppose we conduct experiment with complete randomization, obtain the assignment \(W = (1,0,0,0,0,0,1,1)'\), and the observed outcome for all units are \(Y^{\text{obs}} = (14, 6, 4, 5, 6, 7, 10, 9)'\).

• Impute all the missing potential outcomes under the sharp null hypothesis of no individual causal effects.

W <- c(1,0,0,0,0,0,1,1)

Yobs <- c(14, 6, 4, 5, 6, 7, 10, 9) 

# Under the sharp null of no causal effect
Y0 <- Y1 <- Yobs

• Assume the sharp null hypothesis of no individual causal effects. What is the distribution of the difference-in-means estimator \(\hat{\tau}\).

## observed test statistics value, here you can also replace the difference-in-means estimator with other types of test statistics
tau_obs <- mean(Yobs[W==1]) - mean(Yobs[W==0])

### get all possible assignments
treated_all = combn(8, 3)
W_all = matrix(0, nrow = 56, ncol = 8)
for(m in 1:56){
  W_all[m, treated_all[,m]] = 1
}
dim(W_all)

[1] 56  8

tau_rand <- apply(W_all, 1, function(w) mean(Y1[w==1]) - mean(Y0[w==0]))
hist(tau_rand, breaks = 10, main="Distribution under sharp Null of no individual causal effect", xlab="Obs Diff in Means")

• How extreme is the observed difference-in-means estimator \(\hat{\tau}^{\text{obs}}\)? What is the probability of having at least as extreme as the observed one? More rigorously, specify your test statistic and report your randomization test \(p\)-value.

hist(tau_rand, breaks = 10, main="Distribution under sharp Null of no individual causal effect", xlab="Obs Diff in Means")
abline(v = tau_obs, col = "red")

## calculate p-value 
p <- mean(abs(tau_rand) >= abs(tau_obs))
print(p)

[1] 0.01785714

• Obtain Fisher’s exact p-value under the sharp null hypothesis that the individual treatment effects are all 1.

# Under the sharp null of causal effect is 1
Y0 <- Y1 <- rep(NA, length(Yobs))
Y0[W == 0] <- Yobs[W == 0]
Y0[W == 1] <- Yobs[W == 1] - 1

Y1[W == 1] <- Yobs[W == 1]
Y1[W == 0] <- Yobs[W == 0] + 1

## recompute the distribution of Y1[w = 1] - Y0[w = 0] - 1
tau_rand <- apply(W_all, 1, function(w) mean(Y1[w==1]) - mean(Y0[w==0]))
p = mean(abs(tau_rand - 1) >= abs(tau_obs - 1))
print(p)

[1] 0.01785714

• Calculate the randomization \(p\)-value for a sequence of sharp null hypotheses of additive treatment effects of \(\{-20,-19, \ldots, 19,20\}\), and get the \(95\%\) confidence interval

# need to write a function that compute p-value for any tau
getPvalue <- function(tau, Yobs, W) {
  Y0 <- Y1 <- rep(NA, length(Yobs))
  Y0[W == 0] <- Yobs[W == 0]
  Y0[W == 1] <- Yobs[W == 1] - tau

  Y1[W == 1] <- Yobs[W == 1]
  Y1[W == 0] <- Yobs[W == 0] + tau

  ## recompute the distribution of Y1[w = 1] - Y0[w = 0] - 1
  tau_rand <- apply(W_all, 1, function(w) mean(Y1[w==1]) - mean(Y0[w==0]))
  p = mean(abs(tau_rand - tau) >= abs(tau_obs - tau))
  return(p)
}

## Compute the p-value for a seqeunce of tau
tau_seq <- seq(-20, 20, by = 1)
pvalue_seq <- sapply(tau_seq, getPvalue, Yobs, W)
names(pvalue_seq) <- tau_seq
print(pvalue_seq, digits = 3)

   -20    -19    -18    -17    -16    -15    -14    -13    -12    -11    -10     -9     -8     -7     -6 
0.0179 0.0179 0.0179 0.0179 0.0179 0.0179 0.0179 0.0179 0.0179 0.0179 0.0179 0.0179 0.0179 0.0179 0.0179 
    -5     -4     -3     -2     -1      0      1      2      3      4      5      6      7      8      9 
0.0179 0.0179 0.0179 0.0179 0.0179 0.0179 0.0179 0.0357 0.1786 0.3929 0.8571 0.6964 0.3571 0.1429 0.0536 
    10     11     12     13     14     15     16     17     18     19     20 
0.0357 0.0179 0.0179 0.0179 0.0179 0.0179 0.0179 0.0179 0.0179 0.0179 0.0179 

A rough 95% confidence interval is [2, 10]. You may want a more refined confidence interval by finding the cutoff that is close to 0.05.

## find a better lower bound
tau_seq1 <- seq(2, 3, by = 0.1)
pvalue_seq1 <- sapply(tau_seq1, getPvalue, Yobs, W)
names(pvalue_seq1) <- tau_seq1
print(pvalue_seq1, digits = 3)

     2    2.1    2.2    2.3    2.4    2.5    2.6    2.7    2.8    2.9      3 
0.0357 0.0357 0.0357 0.0357 0.0357 0.0357 0.0357 0.0357 0.0357 0.0357 0.1786 

## find a better upper bound

tau_seq2 <- seq(9, 10, by = 0.1)
pvalue_seq2 <- sapply(tau_seq2, getPvalue, Yobs, W)
names(pvalue_seq2) <- tau_seq2
print(pvalue_seq2, digits = 3)

     9    9.1    9.2    9.3    9.4    9.5    9.6    9.7    9.8    9.9     10 
0.0536 0.0357 0.0357 0.0357 0.0357 0.0357 0.0357 0.0357 0.0357 0.0357 0.0357 

A more refined 95% CI is [2.9, 9.1]

• When the sample size is large, it is impossible to enumerate all possible randomizations. How would you get an approximation for the randomization test \(p\)-value?

## Say we want M = 10000 random samples of the treatment assignments
### get M random samples of possible assignments
getPvalueMC <- function(tau, Yobs, W, M) {
  Y0 <- Y1 <- rep(NA, length(Yobs))
  Y0[W == 0] <- Yobs[W == 0]
  Y0[W == 1] <- Yobs[W == 1] - tau

  Y1[W == 1] <- Yobs[W == 1]
  Y1[W == 0] <- Yobs[W == 0] + tau

  ## randomly sample the treatment assignment vector following the complete randomization mechanism
  
  set.seed(0)
  tau_rand <- sapply(1:M, function(i) {
    w <- sample(W)
    return(mean(Y1[w==1]) - mean(Y0[w==0]))
  })
  p = mean(abs(tau_rand - tau) >= abs(tau_obs - tau))
  return(p)
}

getPvalueMC(0, Yobs, W, 10000)

[1] 0.0171