Code example for post `Randomization-based inference: the Neymanian approach'

cc_neymanian-basic-example.R

library(magrittr)
library(purrr)
library(gtools)


## Estimand function and estimator function
tau.fn <- function(science) mean(science$y1) - mean(science$y0)
hat.tau.fn <- function(w, y) mean(y[w==1]) - mean(y[w==0])


## ----------------------------
## The science table + estimand
## ----------------------------

science <- data.frame(y0 = c(2, 1, 2, 0, 3, 0),
                      y1 = c(3, 1, 4, 1, 3, 2))  

tau <- tau.fn(science)

## ----------------------------
## observed data
## ----------------------------

W.obs  <- c(1, 1, 0, 1, 0, 0)                  
Y.obs  <- science$y1 * W.obs + science$y0 * (1-W.obs)
hat.tau.obs <-  hat.tau.fn(W.obs, Y.obs)


## ----------------------------
## bias of estimator
## ----------------------------

W.ls <- permutations(6, 6, 1:6) %>% array_tree(1) %>% map(~ W.obs[.]) %>% unique

hat.tau.ls <- vector('numeric', length=length(W.ls))

for(i in seq_along(W.ls)) {
    W <- W.ls[[i]]
    Y <- science$y1 * W + science$y0 * (1-W)
    hat.tau.ls[i] <- hat.tau.fn(W, Y)
    print(hat.tau.fn(W, Y))
}

bias <- mean(hat.tau.ls) - tau