Causal Inference I(b) — Covariate Adjustment in Experiments (R companion)¶
estimatr::lm_lin — Lin's reference estimator on two real randomized experiments¶
estimatr (Blair, Cooper, Coppock, Humphreys) is the design-based inference toolkit for experiments. Its lm_lin is the reference implementation of Lin's (2013) covariate-adjusted estimator — center the covariates, interact them with treatment, return heteroskedasticity-robust (HC2) standard errors — and lm_robust / difference_in_means give the unadjusted and ANCOVA comparisons. This companion reproduces the Python notebook on the same two real experiments: The Electric Company (strong pre-test baseline → ~79% variance reduction) and Social Pressure GOTV (weak baseline → ~3%), confirming that adjustment buys precision in exact proportion to $\rho^2$.
1. Precision, not bias — the Electric Company experiment¶
192 classrooms (grades 1–4) randomized in matched pairs to view The Electric Company or not, with reading measured before and after. The pre-test strongly predicts the post-test (r ≈ 0.88), so most outcome variation is baseline ability, not the show. difference_in_means gives the unbiased-but-imprecise estimate.
suppressMessages(library(estimatr))
ec<-read.csv("electric.csv"); ec$X<-NULL
cat(sprintf("Electric Company: %d classrooms (%d treated), corr(pre,post) = %.3f\n", nrow(ec), sum(ec$treatment), cor(ec$pre_test,ec$post_test)))
dm<-difference_in_means(post_test~treatment, data=ec)
cat(sprintf("Unadjusted difference in means: %+.2f points (SE %.2f, 95%% CI [%.2f, %.2f])\n",
dm$coefficients[1], dm$std.error[1], dm$conf.low[1], dm$conf.high[1]))
options(repr.plot.width=13, repr.plot.height=4.6); par(mfrow=c(1,2))
cols<-ifelse(ec$treatment==1,"#2f855a","#2b6cb0")
plot(ec$pre_test, ec$post_test, col=cols, pch=19, xlab="pre-test", ylab="post-test",
main=sprintf("Post-test tracks pre-test (r=%.2f)",cor(ec$pre_test,ec$post_test))); abline(0,1,lty=3,col="grey")
legend("topleft",c("treated","control"),col=c("#2f855a","#2b6cb0"),pch=19,bty="n")
hist(ec$post_test[ec$treatment==0],breaks=18,col=rgb(43,108,176,120,max=255),xlim=range(ec$post_test),
main="Raw outcomes overlap heavily",xlab="post-test"); hist(ec$post_test[ec$treatment==1],breaks=18,col=rgb(47,133,90,120,max=255),add=TRUE)
par(mfrow=c(1,1))
Warning message: "package 'estimatr' was built under R version 4.6.1"
Electric Company: 192 classrooms (96 treated), corr(pre,post) = 0.883
Unadjusted difference in means: +5.66 points (SE 2.54, 95% CI [0.65, 10.66])
2. Regression adjustment (ANCOVA)¶
lm_robust(post_test ~ treatment + pre_test) adds the covariate. Randomization keeps $\hat\tau$ essentially unchanged (unbiased) while the pre-test soaks up residual variance, cutting the standard error sharply — a ~79% reduction in sampling variance.
m0<-lm_robust(post_test~treatment, data=ec)
m1<-lm_robust(post_test~treatment+pre_test, data=ec)
vr<-1-(m1$std.error["treatment"]/m0$std.error["treatment"])^2
cat("Effect of The Electric Company on post-test reading:\n")
cat(sprintf(" unadjusted %+.2f SE %.3f CI width %.2f\n", m0$coefficients["treatment"], m0$std.error["treatment"], m0$conf.high["treatment"]-m0$conf.low["treatment"]))
cat(sprintf(" ANCOVA(pre) %+.2f SE %.3f CI width %.2f\n", m1$coefficients["treatment"], m1$std.error["treatment"], m1$conf.high["treatment"]-m1$conf.low["treatment"]))
cat(sprintf(" --> estimate barely moves, sampling VARIANCE down %.0f%%\n", 100*vr))
options(repr.plot.width=7.5, repr.plot.height=5)
plot(ec$pre_test, ec$post_test, col=ifelse(ec$treatment==1,"#2f855a","#2b6cb0"), pch=19, xlab="pre-test", ylab="post-test",
main=sprintf("ANCOVA: parallel lines, gap = tau = %+.1f", m1$coefficients["treatment"]))
xg<-seq(min(ec$pre_test),max(ec$pre_test),length=50); b<-m1$coefficients
lines(xg, b["(Intercept)"]+b["treatment"]+b["pre_test"]*xg, col="#2f855a", lwd=2.5)
lines(xg, b["(Intercept)"]+b["pre_test"]*xg, col="#2b6cb0", lwd=2.5)
legend("topleft",c("treated","control"),col=c("#2f855a","#2b6cb0"),lwd=2.5,bty="n")
Effect of The Electric Company on post-test reading:
unadjusted +5.66 SE 2.537 CI width 10.01
ANCOVA(pre) +4.73 SE 1.172 CI width 4.63
--> estimate barely moves, sampling VARIANCE down 79%
3. CUPED and the $\rho^2$ law — two real experiments¶
CUPED (Deng et al. 2013) is ANCOVA in A/B-testing dress: $Y^{\text{cuped}}=Y-\theta(X-\bar X)$ with $\theta=\operatorname{Cov}(Y,X)/\operatorname{Var}(X)$, and its variance reduction is exactly $\rho^2$. We confirm it on both experiments — the Electric Company ($\rho^2\approx0.78$) and Social Pressure GOTV ($\rho^2\approx0.03$) — bracketing the law from ~79% down to ~3%.
cuped<-function(df,y,x,t){
th<-cov(df[[y]],df[[x]])/var(df[[x]]); df$yc<-df[[y]]-th*(df[[x]]-mean(df[[x]]))
a<-lm_robust(as.formula(paste(y,"~",t)),data=df); c<-lm_robust(yc~get(t),data=df)
rho<-as.numeric(cor(df[[y]],df[[x]])); vr<-as.numeric(1-(c$std.error[2]/a$std.error[2])^2)
c(rho=rho, rho2=rho^2, vr=vr, eff_u=as.numeric(a$coefficients[2]), eff_c=as.numeric(c$coefficients[2])) }
re<-cuped(ec,"post_test","pre_test","treatment")
sp<-read.csv("social.csv"); sp<-sp[sp$messages %in% c("Control","Neighbors"),]
sp$voted<-as.numeric(sp$primary2006); sp$past<-as.numeric(sp$primary2004); sp$treatment<-as.numeric(sp$messages=="Neighbors")
rs<-cuped(sp,"voted","past","treatment")
cat("CUPED variance reduction tracks rho^2:\n")
cat(sprintf(" Electric Company : rho=%.3f rho^2=%.3f -> var-reduction %.0f%% (effect %+.2f -> %+.2f)\n", re["rho"],re["rho2"],100*re["vr"],re["eff_u"],re["eff_c"]))
cat(sprintf(" Social Pressure : rho=%.3f rho^2=%.3f -> var-reduction %.0f%% (effect %+.3f -> %+.3f)\n", rs["rho"],rs["rho2"],100*rs["vr"],rs["eff_u"],rs["eff_c"]))
options(repr.plot.width=8, repr.plot.height=4.6)
M<-rbind(`rho^2 (predicted)`=c(100*re["rho2"],100*rs["rho2"]), `actual var-reduction`=c(100*re["vr"],100*rs["vr"]))
colnames(M)<-c("Electric\n(strong)","Social Pressure\n(weak)")
bp<-barplot(M,beside=TRUE,col=c("#a0aec0","#2f855a"),ylab="% variance reduction",
main="CUPED variance reduction = rho^2, on real data",legend.text=TRUE,args.legend=list(x="topright",bty="n"))
text(bp, as.vector(M)+2, sprintf("%.0f",as.vector(M)), cex=.8)
CUPED variance reduction tracks rho^2:
Electric Company : rho=0.883 rho^2=0.780 -> var-reduction 79% (effect +5.66 -> +4.73)
Social Pressure : rho=0.163 rho^2=0.026 -> var-reduction 3% (effect +0.081 -> +0.080)
4. The Freedman critique and Lin's fix — lm_lin¶
Freedman (2008): plain ANCOVA can be slightly biased in finite samples and can hurt precision when the covariate–outcome slope differs across arms. Lin (2013) centers the covariates and fully interacts them with treatment, guaranteeing the adjusted estimator is asymptotically no worse than the unadjusted difference. estimatr::lm_lin is its reference implementation. In the huge Social Pressure sample Lin ≈ ANCOVA; in the small Electric Company sample we also respect the blocked (grade) design by adjusting within grade, where per-grade slopes genuinely differ.
le<-lm_lin(post_test~treatment, covariates=~pre_test, data=ec)
ls<-lm_lin(voted~treatment, covariates=~past, data=sp)
cat("Unadjusted vs ANCOVA vs Lin (lm_lin), estimate [SE]:\n")
cat(sprintf(" Electric Company : unadj %+.2f [%.3f] ANCOVA %+.2f [%.3f] Lin %+.2f [%.3f]\n",
m0$coefficients["treatment"],m0$std.error["treatment"], m1$coefficients["treatment"],m1$std.error["treatment"],
le$coefficients["treatment"],le$std.error["treatment"]))
s0<-lm_robust(voted~treatment,sp); s1<-lm_robust(voted~treatment+past,sp)
cat(sprintf(" Social Pressure : unadj %+.3f [%.4f] ANCOVA %+.3f [%.4f] Lin %+.3f [%.4f]\n",
s0$coefficients["treatment"],s0$std.error["treatment"], s1$coefficients["treatment"],s1$std.error["treatment"],
ls$coefficients["treatment"],ls$std.error["treatment"]))
# grade-blocked (Electric's true design)
gr<-do.call(rbind,lapply(sort(unique(ec$grade)),function(g){d<-ec[ec$grade==g,]
a<-lm_robust(post_test~treatment,d); b<-lm_robust(post_test~treatment+pre_test,d)
data.frame(grade=g,unadj=a$coefficients["treatment"],se_u=a$std.error["treatment"],anc=b$coefficients["treatment"],se_a=b$std.error["treatment"])}))
cat("\nElectric Company within grade (blocked design):\n"); print(round(gr,2),row.names=FALSE)
options(repr.plot.width=8, repr.plot.height=4.4)
plot(gr$grade-0.08, gr$unadj, ylim=range(c(gr$unadj-2*gr$se_u,gr$anc+2*gr$se_a)), pch=19, col="#a0aec0",
xlab="grade", ylab="effect on post-test (95% CI)", xaxt="n", main="Adjustment tightens the CI in every grade"); axis(1,at=1:4)
arrows(gr$grade-0.08,gr$unadj-1.96*gr$se_u,gr$grade-0.08,gr$unadj+1.96*gr$se_u,angle=90,code=3,length=.05,col="#a0aec0")
points(gr$grade+0.08, gr$anc, pch=19, col="#2f855a")
arrows(gr$grade+0.08,gr$anc-1.96*gr$se_a,gr$grade+0.08,gr$anc+1.96*gr$se_a,angle=90,code=3,length=.05,col="#2f855a")
abline(h=0,lty=3); legend("topright",c("unadjusted","pre-test adjusted"),col=c("#a0aec0","#2f855a"),pch=19,bty="n")
Unadjusted vs ANCOVA vs Lin (lm_lin), estimate [SE]:
Electric Company : unadj +5.66 [2.537] ANCOVA +4.73 [1.172] Lin +4.73 [1.164]
Social Pressure : unadj +0.081 [0.0027] ANCOVA +0.080 [0.0027] Lin +0.080 [0.0026]
Electric Company within grade (blocked design):
grade unadj se_u anc se_a
1 8.30 4.62 8.79 2.61
2 8.36 2.70 4.27 1.37
3 0.33 2.35 1.91 0.77
4 3.71 1.84 1.70 0.71
5. Summary¶
estimatr reproduced the precision story on two real experiments: difference_in_means gives the unbiased-but-imprecise estimate, lm_robust adds ANCOVA, and lm_lin implements Lin's centered, treatment-interacted estimator with robust SEs. The Electric Company's strong pre-test cut sampling variance ~79% while barely moving the point estimate; Social Pressure GOTV's weak baseline gave ~3% — both exactly $\rho^2$, confirming CUPED's law from real data at both extremes.
Guidance unchanged from the Python notebook: pre-specify strongly predictive baseline covariates (a pre-period outcome is ideal), adjust with lm_lin, report robust/clustered SEs — keeping randomization's unbiasedness while gaining power for free. Cross-links: precision layer on RCT foundations (subsection 1a); the pre-period-covariate idea recurs as the baseline in DiD (subsection 6) and as residualization in Double ML (subsection 9).