When does a difference in means not capture the true treatment effects vs a regression with pre-treatment controls? A question from Gelman - Regression & Other Stories...
In answer to my own question...my understanding is that a difference in means should not capture the treatment effect when there are pre-treatment variables that are highly predictive(associate) with the outcome, i.e., confounding, or in a situation with selection bias.
However, I am somewhat struggling to produce this situation via simulation (as opposed to drawing as the question asks). See my code below. I don't understand why there isn't a bigger difference in the bias of the OLS estimate and diff in means. The inclusion of the pre-treatment variable in creating the outcome suggests that a difference in means should be more biased...
'Linearity assumptions and causal inference: Consider a study with an outcome, y, a treatment indicator, z, and a single confounding covariate, x. Draw a scatterplot of treatment and control observations that demonstrates each of the following:
(a) A scenario where the difference in means estimate would not capture the true treatment
effect but a regression of y on x and z would yield the correct estimate.'
library(dplyr)

experiment = function(N, T){
  n = N
  x1 = rnorm(n, 65, 3 )
  t = T
  z = sample(rep(c(1,0), each = n/2), replace = FALSE)
  y0 = 10 +1.3*x1  + rnorm(n, 0, 1)
  y1 = 10 +1.3*x1  + t + rnorm(n, 0, 1)
  y = (1-z)*y0 + z*y1
  data = data_frame(x1, z,y0,y1, y) 
  fit = lm(y ~ z  + x1 , data = data)
  dm_res = mean(data[data$z == 1 ,]$y) - mean(data[data$z==0,]$y)
  ols = fit$coefficients[[2]]
  results = c('ols' = ols, 'meandiff' = dm_res)
  return(results)
}

vals = seq(from = .1, to = 5,by = .1)
l = list()
for (val in vals){

  #replicate experiment 100 times and take average value of estimates
  res = replicate(100, experiment(2000, val)) %>% rowMeans()
  #subtract true value from estimated value
  l[[paste(val)]] = res - val
}

l %>% bind_rows() %>% colMeans()



There is no confounding in your example because the treatment was generated completely independently of the predictor of the outcome. You basically implemented a randomized trial. In such a scenario, both the difference in means and covariate-adjusted estimate are unbiased for the treatment effect (and therefore for each other). To induce confounding, you need the covariate to cause selection into treatment. Something like the following would do that:
z <- rbinom(n, 1, plogis(.6 * x1))

This creates a dependence between the covariate and treatment that would yield confounding.
