Simulate covariate adjustment

covariate_sim.R

# based on http://egap.org/methods-guides/10-things-know-about-covariate-adjustment
library(MASS)  # for mvrnorm()
library(tidyverse)

set.seed(1234567)
num.reps = 1000

# True treatment effect is 0 for every unit

adj.est = function(n, cov.matrix, treated) {
  Y.and.X  =  mvrnorm(n, mu = c(0, 0), Sigma = cov.matrix)
  Y   =  Y.and.X[, 1]
  X   =  Y.and.X[, 2]
  coef(lm(Y ~ treated + X))[2]
}

unadj.est = function(n, treated) {
  Y = rnorm(n)
  coef(lm(Y ~ treated))[2]
}

rmse = function(half.n, rho = 0, control = TRUE) {
  treated  =  rep(c(0, 1), half.n)
  n = 2 * half.n

  if (control) {
    cov.matrix  =  matrix(c(1, rho, rho, 1), nrow = 2, ncol = 2)
    return( sqrt(mean(replicate(num.reps, adj.est(n, cov.matrix, treated)) ^ 2)) )
  }
  else {
    return( sqrt(mean(replicate(num.reps, unadj.est(n, treated)) ^ 2)) )
  }
}

d <- expand.grid(half_n =  c(5, 10, 20, 50, 100, 200),
                 rho = seq(0, .8, .2)) %>%
  rowwise() %>%
  mutate(rmse = rmse(half_n, rho))

ggplot(d, aes(x = half_n, y = rmse, 
              col = factor(rho))) +
  geom_line() + 
  theme_classic() + 
  ggthemes::scale_color_solarized(name = "corr(X, Z)") + 
  xlab("N per group") +
  ylab("Root mean squared error") + 
  theme(legend.position = "bottom")