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

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

## 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]))

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