simmons_replication.R

library(tidyverse)
# Replication of Simmons ####
set.seed(42)
reps <- 1500
n <- 40
sim <- matrix(NA, nrow = reps, ncol = 4)
colnames(sim) <- c("vanilla_est", "vanilla_p", "covariate_est", "covariate_p")

for (i in 1:reps){
  d = data_frame(
    treatment = sample(0:1, size = n, replace = TRUE),
    gender = sample(c(-.5,.5), size = n, replace = TRUE),
    out = rnorm(n = n) # + .5*gender
  )
  m_vanilla = lm(out ~ treatment, data = d)
  sim[i, 1:2] = summary(m_vanilla) %>% `[`("coefficients") %>% `[[`(1) %>% `[`(2,c(1,4))
  m_cov = lm(out ~ treatment*gender, data = d)
  sim[i, 3] = summary(m_cov) %>% `[`("coefficients") %>% `[[`(1) %>% `[`(2,1)
  tmp = summary(m_cov) %>% `[`("coefficients") %>% `[[`(1) 
  sim[i, 4] = min(tmp[c(2,4),4])
}

sim %>% as_data_frame() %>%
  summarise(
    eff_vanilla = sum(vanilla_p < .05)/reps*100,
    eff_covariate = sum(covariate_p < .05)/reps*100,
    diff = mean(vanilla_est - covariate_est)
  )

plot(density(sim[,1]-sim[,3]), type='l')

# Extension of Simmons: Ommitted Variable Bias example  ####
set.seed(42)
reps <- 1500
n <- 40
sim <- matrix(NA, nrow = reps, ncol = 5)
colnames(sim) <- c("vanilla_est", "vanilla_p", "covariate_est", "covariate_p_no_int","covariate_p_int")

for (i in 1:reps){
  d = data_frame(
    gender = sample(c(-.5,.5), size = n, replace = TRUE),
    treatment = ifelse(gender == .5, sample(0:1, size = n, replace = TRUE, prob = c(.6,.4)),
                       sample(0:1, size = n, replace = TRUE, prob = c(.4,.6))
                       ),
    out = rnorm(n = n)  + .8*gender
  )
  m_vanilla = lm(out ~ treatment, data = d)
  sim[i, 1:2] = summary(m_vanilla) %>% `[`("coefficients") %>% `[[`(1) %>% `[`(2,c(1,4))
  m_cov = lm(out ~ treatment*gender, data = d)
  sim[i, 3:4] = summary(m_cov) %>% `[`("coefficients") %>% `[[`(1) %>% `[`(2,c(1,4))
  tmp = summary(m_cov) %>% `[`("coefficients") %>% `[[`(1) 
  sim[i, 5] = min(tmp[c(2,4),4])
}

sim %>% as_data_frame() %>%
  summarise(
    eff_vanilla = sum(vanilla_p < .05)/reps*100,
    eff_covariate_no_int = sum(covariate_p_no_int < .05)/reps*100,
    eff_covariate_int = sum(covariate_p_int < .05)/reps*100,
    diff = mean(vanilla_est - covariate_est) 
  )


plot(density(sim[,1]-sim[,3]), type='l')