Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Simulate covariate adjustment
# 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")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.