Causal Inference I(b) — Covariate Adjustment in Experiments¶

Precision, not bias: ANCOVA, CUPED, and the Freedman–Lin debate — on two real randomized experiments¶

Randomization already gives an unbiased treatment effect from a simple difference in means (subsection 1a). So why adjust for covariates at all? Not to remove bias — there is none to remove — but to buy precision: a pre-treatment variable that predicts the outcome soaks up noise, shrinking the standard error and tightening the confidence interval without changing what the estimate targets. This is the everyday workhorse of modern experimentation, from clinical trials (ANCOVA) to tech A/B testing (CUPED).

We work with two real randomized experiments that bracket the phenomenon:

  • The Electric Company (1971) — a classic education RCT: 192 classrooms (grades 1–4) randomized to watch an educational TV show or not, with a reading pre-test and post-test. The pre-test is a strong predictor of the post-test — the ideal case for adjustment.
  • Social Pressure GOTV (Gerber, Green & Larimer 2008) — a 305,866-voter field experiment mailing get-out-the-vote messages; here the only baseline covariate (voted in the 2004 primary) is a weak predictor of 2006 turnout.

The two together demonstrate the governing law — variance reduction ≈ $\rho^2$, the squared outcome–covariate correlation — at both ends: a dramatic ~79% reduction where the baseline is strong, a modest ~3% where it is weak. Python leads (from-scratch OLS/ANCOVA/CUPED/Lin); the R companion uses estimatr's lm_lin. This is the precision companion to the RCT-foundations notebook.

1. Precision, not bias — the Electric Company experiment¶

192 classrooms were randomized (in matched pairs within each grade) to either view The Electric Company or serve as controls; reading was measured before and after the year. The simple difference in mean post-test scores is already an unbiased estimate of the treatment effect. But the pre-test score is strongly correlated with the post-test (r ≈ 0.88): most of the variation in outcomes is explained by where a class started, not by the show. That predictable variation is noise for the purpose of estimating the treatment effect — and a covariate that captures it can be removed.

In [1]:
import numpy as np, pandas as pd, matplotlib.pyplot as plt, statsmodels.formula.api as smf, warnings
warnings.filterwarnings("ignore")
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
ec=pd.read_csv("electric.csv", index_col=0)
print(f"Electric Company: {len(ec)} classrooms, {int(ec.treatment.sum())} treated / {int((1-ec.treatment).sum())} control, grades {sorted(ec.grade.unique())}")
print(f"corr(pre-test, post-test) = {ec.pre_test.corr(ec.post_test):.3f}   (a strong baseline predictor)")
m0=smf.ols("post_test~treatment", data=ec).fit(cov_type="HC1")
print(f"\nUnadjusted difference in means: {m0.params['treatment']:+.2f} points  (SE {m0.bse['treatment']:.2f}, 95% CI [{m0.conf_int().loc['treatment',0]:.2f}, {m0.conf_int().loc['treatment',1]:.2f}])")
fig,ax=plt.subplots(1,2,figsize=(13,4.6))
for t,c,lab in [(1,GREEN,"treated (watched show)"),(0,BLUE,"control")]:
    g=ec[ec.treatment==t]; ax[0].scatter(g.pre_test,g.post_test,c=c,alpha=.6,s=28,label=lab)
ax[0].plot([0,120],[0,120],ls=":",c=GREY); ax[0].set_xlabel("pre-test score"); ax[0].set_ylabel("post-test score")
ax[0].set_title(f"Post-test is mostly predicted by pre-test (r={ec.pre_test.corr(ec.post_test):.2f})"); ax[0].legend(fontsize=8)
ax[1].hist(ec[ec.treatment==0].post_test,bins=18,alpha=.55,color=BLUE,label="control"); ax[1].hist(ec[ec.treatment==1].post_test,bins=18,alpha=.55,color=GREEN,label="treated")
ax[1].axvline(ec[ec.treatment==0].post_test.mean(),color=BLUE,lw=2); ax[1].axvline(ec[ec.treatment==1].post_test.mean(),color=GREEN,lw=2)
ax[1].set_xlabel("post-test score"); ax[1].set_title("Raw outcome distributions overlap heavily"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
print("The unadjusted effect is unbiased but IMPRECISE -- the treated/control post-test distributions overlap because most")
print("of the spread is baseline reading ability, not the show. Adjusting for the pre-test will remove exactly that spread.")

Electric Company: 192 classrooms, 96 treated / 96 control, grades [np.int64(1), np.int64(2), np.int64(3), np.int64(4)]
corr(pre-test, post-test) = 0.883   (a strong baseline predictor)

Unadjusted difference in means: +5.66 points  (SE 2.54, 95% CI [0.68, 10.63])
No description has been provided for this image 
The unadjusted effect is unbiased but IMPRECISE -- the treated/control post-test distributions overlap because most
of the spread is baseline reading ability, not the show. Adjusting for the pre-test will remove exactly that spread.

2. Regression adjustment (ANCOVA)¶

The standard fix is analysis of covariance: regress the outcome on treatment and the pre-treatment covariate, $$Y_i=\alpha+\tau\,T_i+\beta\,X_i+\varepsilon_i.$$ Because randomization makes $X$ (approximately) orthogonal to $T$, adding $X$ does not move $\hat\tau$ (it stays unbiased) but it removes $\beta X$ worth of residual variance, shrinking $\operatorname{SE}(\hat\tau)$. Here the point estimate barely changes while the standard error more than halves — a ~79% reduction in sampling variance. The small shift that does occur corrects a slight pre-test imbalance across grades, which is the other thing adjustment does: it removes chance imbalance in the covariate.

In [2]:
m1=smf.ols("post_test~treatment+pre_test", data=ec).fit(cov_type="HC1")
vr=1-(m1.bse['treatment']/m0.bse['treatment'])**2
print("Effect of The Electric Company on post-test reading:")
print(f"  unadjusted (diff in means)   {m0.params['treatment']:+.2f}   SE {m0.bse['treatment']:.3f}   CI width {m0.conf_int().loc['treatment',1]-m0.conf_int().loc['treatment',0]:.2f}")
print(f"  ANCOVA (adjust for pre-test) {m1.params['treatment']:+.2f}   SE {m1.bse['treatment']:.3f}   CI width {m1.conf_int().loc['treatment',1]-m1.conf_int().loc['treatment',0]:.2f}")
print(f"  --> point estimate barely moves (unbiased), SE cut by {100*(1-m1.bse['treatment']/m0.bse['treatment']):.0f}%, sampling VARIANCE down {100*vr:.0f}%")
# bootstrap the sampling distributions to make 'precision' visible
rng=np.random.default_rng(0); bu=[]; ba=[]
for _ in range(2000):
    s=ec.sample(len(ec),replace=True,random_state=rng.integers(1e9))
    bu.append(smf.ols("post_test~treatment",data=s).fit().params['treatment'])
    ba.append(smf.ols("post_test~treatment+pre_test",data=s).fit().params['treatment'])
fig,ax=plt.subplots(1,2,figsize=(13,4.4))
xg=np.linspace(min(ec.pre_test),max(ec.pre_test),50)
for t,c in [(1,GREEN),(0,BLUE)]:
    g=ec[ec.treatment==t]; ax[0].scatter(g.pre_test,g.post_test,c=c,alpha=.5,s=22)
    ax[0].plot(xg, m1.params['Intercept']+m1.params['treatment']*t+m1.params['pre_test']*xg, c=c, lw=2.5)
ax[0].set_xlabel("pre-test"); ax[0].set_ylabel("post-test"); ax[0].set_title(f"ANCOVA: two parallel lines, gap = tau = {m1.params['treatment']:+.1f}")
ax[1].hist(bu,bins=40,alpha=.55,color=GREY,label=f"unadjusted (SD {np.std(bu):.2f})")
ax[1].hist(ba,bins=40,alpha=.65,color=GREEN,label=f"ANCOVA (SD {np.std(ba):.2f})")
ax[1].axvline(0,color="k",lw=.6); ax[1].set_xlabel("estimated effect (bootstrap)"); ax[1].set_title("Adjusted estimator has a far tighter sampling distribution"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
print(f"Bootstrap SD: unadjusted {np.std(bu):.2f} vs ANCOVA {np.std(ba):.2f} -- the adjusted estimate is ~{np.std(bu)/np.std(ba):.1f}x more precise, same target.")

Effect of The Electric Company on post-test reading:
  unadjusted (diff in means)   +5.66   SE 2.537   CI width 9.94
  ANCOVA (adjust for pre-test) +4.73   SE 1.169   CI width 4.58
  --> point estimate barely moves (unbiased), SE cut by 54%, sampling VARIANCE down 79%
No description has been provided for this image 
Bootstrap SD: unadjusted 2.47 vs ANCOVA 1.15 -- the adjusted estimate is ~2.1x more precise, same target.

3. CUPED and the $\rho^2$ law — confirmed on two real experiments¶

Tech A/B testing rediscovered ANCOVA as CUPED (Controlled-experiment Using Pre-Experiment Data, Deng et al. 2013): replace the outcome with a covariate-adjusted version $$Y_i^{\text{cuped}}=Y_i-\theta\,(X_i-\bar X),\qquad \theta=\frac{\operatorname{Cov}(Y,X)}{\operatorname{Var}(X)},$$ then take the simple difference in means of $Y^{\text{cuped}}$. Because $\theta$ is the regression slope, CUPED equals ANCOVA — and its variance reduction is exactly $$1-\frac{\operatorname{Var}(\hat\tau_{\text{cuped}})}{\operatorname{Var}(\hat\tau)}\approx\rho^2,$$ the squared outcome–covariate correlation. We confirm this on both experiments: the Electric Company ($\rho\approx0.88$, so $\rho^2\approx0.78$) yields ~79% reduction, while Social Pressure GOTV ($\rho\approx0.16$, $\rho^2\approx0.026$) yields ~3%. The lesson is precise and general: adjustment helps exactly as much as the covariate predicts the outcome — no more, no less.

In [3]:
def cuped(df,y,x,t):
    th=np.cov(df[y],df[x])[0,1]/np.var(df[x],ddof=1)
    dd=df.copy(); dd["yc"]=dd[y]-th*(dd[x]-dd[x].mean())
    a=smf.ols(f"{y}~{t}",data=dd).fit(cov_type="HC1"); c=smf.ols(f"yc~{t}",data=dd).fit(cov_type="HC1")
    rho=df[y].corr(df[x]); vr=1-(c.bse[t]/a.bse[t])**2
    return th,rho,vr,a.params[t],c.params[t],a.bse[t],c.bse[t]
# Electric Company
th_e,rho_e,vr_e,tu_e,tc_e,su_e,sc_e=cuped(ec,"post_test","pre_test","treatment")
# Social Pressure GOTV -- Neighbors vs Control, covariate = voted in 2004 primary
sp=pd.read_csv("social.csv"); sp["voted"]=sp.primary2006.astype(float); sp["past"]=sp.primary2004.astype(float)
sp=sp[sp.messages.isin(["Control","Neighbors"])].copy(); sp["treatment"]=(sp.messages=="Neighbors").astype(float)
th_s,rho_s,vr_s,tu_s,tc_s,su_s,sc_s=cuped(sp,"voted","past","treatment")
print("CUPED = regression adjustment; variance reduction tracks rho^2:")
print(f"  Electric Company : rho={rho_e:.3f}  rho^2={rho_e**2:.3f}  ->  variance reduction {100*vr_e:.0f}%   (effect {tu_e:+.2f} -> {tc_e:+.2f})")
print(f"  Social Pressure  : rho={rho_s:.3f}  rho^2={rho_s**2:.3f}  ->  variance reduction {100*vr_s:.0f}%   (effect {tu_s:+.3f} -> {tc_s:+.3f})")
fig,ax=plt.subplots(1,2,figsize=(13,4.4))
labs=["Electric Company\n(strong baseline)","Social Pressure GOTV\n(weak baseline)"]
x=np.arange(2); w=0.36
ax[0].bar(x-w/2,[100*rho_e**2,100*rho_s**2],w,color=GREY,label="rho^2 (predicted)")
ax[0].bar(x+w/2,[100*vr_e,100*vr_s],w,color=GREEN,label="actual variance reduction")
ax[0].set_xticks(x); ax[0].set_xticklabels(labs,fontsize=8); ax[0].set_ylabel("% variance reduction"); ax[0].set_title("CUPED variance reduction = rho^2, on real data"); ax[0].legend(fontsize=8)
for i,(a_,b_) in enumerate(zip([rho_e**2,rho_s**2],[vr_e,vr_s])): ax[0].text(i-w/2,100*a_+1,f"{100*a_:.0f}",ha="center",fontsize=8); ax[0].text(i+w/2,100*b_+1,f"{100*b_:.0f}",ha="center",fontsize=8)
rr=np.linspace(0,1,100); ax[1].plot(rr,rr,color=GREEN,lw=2)
ax[1].scatter([rho_e**2,rho_s**2],[vr_e,vr_s],c=[BLUE,ORANGE],s=90,zorder=3)
ax[1].annotate("Electric",(rho_e**2,vr_e),textcoords="offset points",xytext=(-55,0),fontsize=8,color=BLUE)
ax[1].annotate("Social Pressure",(rho_s**2,vr_s),textcoords="offset points",xytext=(10,-4),fontsize=8,color=ORANGE)
ax[1].set_xlabel("rho^2"); ax[1].set_ylabel("variance reduction"); ax[1].set_title("Both experiments sit on the 45-degree line"); ax[1].set_xlim(-.03,1); ax[1].set_ylim(-.03,1)
plt.tight_layout(); plt.show()
print("Same law at both extremes: a strong baseline (Electric) buys ~79%, a weak one (Social Pressure) ~3% -- exactly rho^2.")

CUPED = regression adjustment; variance reduction tracks rho^2:
  Electric Company : rho=0.883  rho^2=0.780  ->  variance reduction 79%   (effect +5.66 -> +4.73)
  Social Pressure  : rho=0.163  rho^2=0.026  ->  variance reduction 3%   (effect +0.081 -> +0.080)
No description has been provided for this image 
Same law at both extremes: a strong baseline (Electric) buys ~79%, a weak one (Social Pressure) ~3% -- exactly rho^2.

4. The Freedman critique and Lin's fix¶

Freedman (2008) pointed out an uncomfortable fact: plain ANCOVA can be slightly biased in finite samples and can even hurt precision if the covariate's relationship to the outcome differs across arms, because a single common slope $\beta$ is imposed on both groups. Lin (2013) resolved it with a small change — center the covariates and fully interact them with treatment: $$Y_i=\alpha+\tau\,T_i+\beta\,(X_i-\bar X)+\gamma\,T_i(X_i-\bar X)+\varepsilon_i.$$ Fitting a separate slope per arm guarantees Lin's estimator is at least as precise as the unadjusted difference asymptotically and removes Freedman's bias, at the cost of a few extra parameters. In large samples the interaction rarely matters (Social Pressure: Lin ≈ ANCOVA to the fourth decimal); in small samples with heterogeneous slopes it is the safe default. Here we also respect the Electric Company's design by adjusting within grade (its randomization was blocked by grade), where per-grade pre/post slopes genuinely differ.

In [4]:
def lin(df,y,x,t):
    dd=df.copy(); dd["xc"]=dd[x]-dd[x].mean()
    return smf.ols(f"{y}~{t}+xc+{t}:xc",data=dd).fit(cov_type="HC1")
le=lin(ec,"post_test","pre_test","treatment"); ls=lin(sp,"voted","past","treatment")
print("Unadjusted vs ANCOVA vs Lin (interacted), point estimate [SE]:")
print(f"  Electric Company : unadj {tu_e:+.2f} [{su_e:.3f}]   ANCOVA {smf.ols('post_test~treatment+pre_test',ec).fit(cov_type='HC1').params['treatment']:+.2f} [{m1.bse['treatment']:.3f}]   Lin {le.params['treatment']:+.2f} [{le.bse['treatment']:.3f}]")
print(f"  Social Pressure  : unadj {tu_s:+.3f} [{su_s:.4f}]   ANCOVA {smf.ols('voted~treatment+past',sp).fit(cov_type='HC1').params['treatment']:+.3f} [{smf.ols('voted~treatment+past',sp).fit(cov_type='HC1').bse['treatment']:.4f}]   Lin {ls.params['treatment']:+.3f} [{ls.bse['treatment']:.4f}]")
# grade-blocked adjustment for Electric (its true design)
rows=[]
for g in sorted(ec.grade.unique()):
    dg=ec[ec.grade==g]; a=smf.ols("post_test~treatment",dg).fit(cov_type="HC1"); b=smf.ols("post_test~treatment+pre_test",dg).fit(cov_type="HC1")
    rows.append((g,a.params['treatment'],a.bse['treatment'],b.params['treatment'],b.bse['treatment']))
gr=pd.DataFrame(rows,columns=["grade","unadj","se_u","anc","se_a"])
print("\nElectric Company within grade (its blocked design) -- adjustment tightens every grade:")
print(gr.round(2).to_string(index=False))
fig,ax=plt.subplots(figsize=(8,4.2))
ax.errorbar(gr.grade-0.08,gr.unadj,yerr=1.96*gr.se_u,fmt="o",color=GREY,capsize=4,label="unadjusted")
ax.errorbar(gr.grade+0.08,gr.anc,yerr=1.96*gr.se_a,fmt="o",color=GREEN,capsize=4,label="pre-test adjusted")
ax.axhline(0,color="k",lw=.6); ax.set_xlabel("grade"); ax.set_ylabel("effect on post-test (95% CI)"); ax.set_xticks([1,2,3,4])
ax.set_title("Adjustment shrinks the CI in every grade (largest effect in grade 1)"); ax.legend(fontsize=8)
plt.tight_layout(); plt.show()
print("Lin = center + interact: guaranteed no worse than unadjusted, fixes Freedman's small-sample bias. Use it by default.")

Unadjusted vs ANCOVA vs Lin (interacted), point estimate [SE]:
  Electric Company : unadj +5.66 [2.537]   ANCOVA +4.73 [1.169]   Lin +4.73 [1.157]

  Social Pressure  : unadj +0.081 [0.0027]   ANCOVA +0.080 [0.0027]   Lin +0.080 [0.0026]

Electric Company within grade (its blocked design) -- adjustment tightens every grade:
 grade  unadj  se_u  anc  se_a
     1   8.30  4.62 8.79  2.61
     2   8.36  2.70 4.27  1.36
     3   0.33  2.35 1.91  0.77
     4   3.71  1.84 1.70  0.71
No description has been provided for this image 
Lin = center + interact: guaranteed no worse than unadjusted, fixes Freedman's small-sample bias. Use it by default.

5. Summary¶

Covariate adjustment in a randomized experiment is about precision, not identification — the effect is already unbiased, and adjustment simply removes predictable outcome variance:

  • ANCOVA (regress $Y$ on $T$ and a pre-treatment $X$) left the Electric Company estimate essentially unchanged (+5.7 → +4.7 points) while cutting the sampling variance ~79% — the treated/control comparison went from imprecise to sharp because the pre-test explained most of the outcome spread.
  • CUPED is the same estimator (the A/B-testing formulation), and its variance reduction is exactly $\rho^2$ — confirmed at both extremes on real data: ~79% where the baseline is strong (Electric, $\rho^2=0.78$), ~3% where it is weak (Social Pressure, $\rho^2=0.03$). Adjustment pays off in exact proportion to how well the covariate predicts the outcome.
  • Freedman's critique + Lin's fix: plain ANCOVA can be slightly biased in small samples; Lin's centered, treatment-interacted regression is guaranteed asymptotically no worse than the unadjusted estimate and is the safe default.

Practical guidance: pre-specify a small set of strongly predictive baseline covariates (a pre-period measure of the outcome is best), adjust with Lin's interacted estimator, and report cluster/heteroskedasticity-robust SEs. You keep randomization's unbiasedness and gain power for free. Cross-links: this is the precision layer on the RCT foundations (subsection 1a); the same pre-period-covariate logic reappears as the baseline period in Difference-in-Differences (subsection 6) and as the residualization idea behind Double Machine Learning (subsection 9); the Electric Company's blocked (grade) design connects to stratified/matched-pair randomization. The R companion runs both experiments through estimatr's lm_lin (Lin's reference implementation).