Causal Inference II — Potential Outcomes & Matching¶
Recovering an experiment from observational data — the LaLonde job-training problem¶
In the previous notebook randomization did the identifying work for free: treatment was independent of the potential outcomes by design. Almost no real policy question comes with that luxury. When we only observe who happened to be treated, the treated and control groups differ systematically, and the naive comparison confounds the treatment effect with those pre-existing differences. This notebook introduces the assumption that stands in for randomization — unconfoundedness — and the family of methods built on it: propensity scores, matching, inverse-probability weighting, and doubly-robust estimation.
We use the most famous test case in the field, and one with a built-in answer key. Robert LaLonde (1986) asked a devastating question: if we take a program whose true effect we know from a randomized experiment, throw away the experimental controls, and try to recover that effect from observational comparison data, do our econometric methods succeed? His answer was largely "no." Dehejia & Wahba (1999) revived the debate by showing that propensity-score methods could recover the experimental benchmark after all. We replay that whole argument here:
- the experimental benchmark — the true effect, from the randomized NSW trial;
- the observational trap — the same treated units compared to a survey control group, where the naive estimate is badly wrong;
- the fixes — propensity scores, nearest-neighbor and Mahalanobis matching, covariate-balance diagnostics, IPW, and doubly-robust AIPW — and how close each gets to the truth.
Python-lead (from-scratch estimators + scikit-learn); the R companion uses MatchIt, cobalt, and WeightIt.
1. The data, the benchmark, and the confounding problem¶
The National Supported Work Demonstration (NSW) was a mid-1970s randomized program that gave 6–18 months of guaranteed employment and training to severely disadvantaged workers — long-term welfare recipients, ex-addicts, ex-offenders, and high-school dropouts. The outcome is real annual earnings in 1978 (re78), after the program. The covariates are pre-treatment: age, education, race (black/hispanic), marital status, a no-high-school-degree indicator, and real earnings in 1974 and 1975 (re74,re75).
Because assignment was randomized, the experimental sample (185 treated + 260 experimental controls, the Dehejia-Wahba subset) identifies the true effect by a simple difference in means — exactly the logic of the previous notebook. That number is our answer key: an average treatment effect on the treated (ATT) of about USD 1,794.
LaLonde's challenge: discard the 260 experimental controls and instead compare the 185 treated to 429 controls drawn from the Current Population Survey (CPS) — a general survey sample. These CPS respondents are nothing like the NSW participants (they are older, better-educated, and far higher-earning), so the naive difference in means is badly confounded. We show that raw gap and the covariate imbalance that causes it. The estimand throughout is the ATT — the effect for those who actually took the training — which is the policy-relevant quantity and what the experimental benchmark measures.
import numpy as np, pandas as pd, matplotlib.pyplot as plt, warnings
warnings.filterwarnings("ignore")
from sklearn.linear_model import LogisticRegression, LinearRegression
from scipy.spatial.distance import cdist
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
cov=["age","educ","black","hispan","married","nodegree","re74","re75"]
exp=pd.read_csv("lalonde_exp.csv"); obs=pd.read_csv("lalonde_obs.csv")
bench=exp.re78[exp.treat==1].mean()-exp.re78[exp.treat==0].mean() # experimental ATT (randomized)
W=obs.treat.values; Y=obs.re78.values; X=obs[cov].values.astype(float)
naive=Y[W==1].mean()-Y[W==0].mean()
print(f"EXPERIMENTAL benchmark (randomized NSW, n={len(exp)}: {(exp.treat==1).sum()} treated, "
f"{(exp.treat==0).sum()} experimental controls): ATT = ${bench:,.0f} <-- the truth")
print(f"OBSERVATIONAL sample (185 NSW treated + 429 CPS controls): naive diff = ${naive:,.0f} <-- confounded")
print(f" -> the naive observational comparison suggests the program REDUCED earnings; bias = ${naive-bench:,.0f}\n")
def smd(col): # standardized mean diff (treated sd)
t=X[W==1,col]; c=X[W==0,col]; return (t.mean()-c.mean())/t.std()
tab=pd.DataFrame({"treated (NSW)":X[W==1].mean(0),"control (CPS)":X[W==0].mean(0),
"std.mean.diff":[smd(j) for j in range(len(cov))]},index=cov)
print("Covariate means — NSW treated vs CPS controls:"); print(tab.round(2).to_string())
fig,ax=plt.subplots(1,2,figsize=(13,4.4))
ax[0].barh(cov,[smd(j) for j in range(len(cov))],color=[RED if abs(smd(j))>0.1 else GREY for j in range(len(cov))])
ax[0].axvline(0,color="k",lw=.7); ax[0].axvline(-0.1,color=GREY,ls=":"); ax[0].axvline(0.1,color=GREY,ls=":")
ax[0].set_xlabel("standardized mean difference (treated − control)"); ax[0].set_title("Severe covariate imbalance (|SMD|>0.1 = red)")
ax[1].hist(Y[W==0],bins=40,alpha=.6,color=BLUE,density=True,label="CPS controls"); ax[1].hist(Y[W==1],bins=40,alpha=.6,color=RED,density=True,label="NSW treated")
ax[1].set_xlabel("1978 earnings (re78, USD)"); ax[1].set_title("CPS controls earn far more — the confounding"); ax[1].legend()
plt.tight_layout(); plt.show()
print("The CPS controls are older, more educated, far more likely married, and earned thousands more BEFORE the program.")
print("Comparing NSW trainees to them without adjustment is the confounding the rest of this notebook must undo.")
EXPERIMENTAL benchmark (randomized NSW, n=445: 185 treated, 260 experimental controls): ATT = $1,794 <-- the truth
OBSERVATIONAL sample (185 NSW treated + 429 CPS controls): naive diff = $-635 <-- confounded
-> the naive observational comparison suggests the program REDUCED earnings; bias = $-2,429
Covariate means — NSW treated vs CPS controls:
treated (NSW) control (CPS) std.mean.diff
age 25.82 28.03 -0.31
educ 10.35 10.24 0.06
black 0.84 0.20 1.76
hispan 0.06 0.14 -0.35
married 0.19 0.51 -0.83
nodegree 0.71 0.60 0.24
re74 2095.57 5619.24 -0.72
re75 1532.06 2466.48 -0.29
The CPS controls are older, more educated, far more likely married, and earned thousands more BEFORE the program. Comparing NSW trainees to them without adjustment is the confounding the rest of this notebook must undo.
2. Unconfoundedness and the propensity score¶
To learn the effect from observational data we invoke unconfoundedness (a.k.a. conditional independence, selection on observables): conditional on the covariates $X$, treatment is as good as random, $$\big(Y(1),Y(0)\big)\ \perp\ W \ \mid\ X,$$ together with overlap ($0<e(X)<1$ for all $X$). This is the assumption that replaces randomization — and, unlike randomization, it is not testable: it asserts there is no unobserved confounder. All we can check is balance on the covariates we do have.
Conditioning on the full covariate vector is hard when $X$ is high-dimensional. Rosenbaum & Rubin (1983) proved a remarkable simplification: the scalar propensity score $e(X)=P(W=1\mid X)$ is a balancing score — if unconfoundedness holds given $X$, it also holds given $e(X)$ alone. So we can match or weight on one number instead of eight. We estimate $e(X)$ by logistic regression and inspect overlap: the region where treated and control propensity distributions coincide. The CPS controls pile up near $e\approx0$ (they look nothing like trainees), so only a subset provides valid comparisons — the rest are off the common support.
Which covariates enter the propensity model is itself a variable-selection problem — include every confounder, avoid post-treatment variables and instruments — connecting directly to the Variable Selection arc.
Xs=(X-X.mean(0))/X.std(0) # standardize for the PS model
ps=LogisticRegression(penalty=None,max_iter=5000).fit(Xs,W).predict_proba(Xs)[:,1]
obs["ps"]=ps
fig,ax=plt.subplots(1,2,figsize=(13,4.4))
ax[0].hist(ps[W==0],bins=40,alpha=.6,color=BLUE,density=True,label="CPS controls"); ax[0].hist(ps[W==1],bins=40,alpha=.6,color=RED,density=True,label="NSW treated")
ax[0].set_xlabel("estimated propensity score e(X)"); ax[0].set_ylabel("density"); ax[0].set_title("Propensity overlap: most CPS controls look nothing like trainees"); ax[0].legend()
# common support region
lo,hi=ps[W==1].min(),ps[W==1].max()
ax[1].scatter(ps[W==0],np.random.uniform(0,0.4,(W==0).sum()),s=6,color=BLUE,alpha=.3,label="CPS controls")
ax[1].scatter(ps[W==1],np.random.uniform(0.6,1.0,(W==1).sum()),s=8,color=RED,alpha=.5,label="NSW treated")
ax[1].axvspan(0.05,0.95,color=GREEN,alpha=.08); ax[1].set_xlabel("propensity score"); ax[1].set_yticks([]); ax[1].set_title("Only where the two overlap can effects be estimated"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
n_supp=((ps[W==0]>=0.05)&(ps[W==0]<=0.95)).sum()
print(f"Of {(W==0).sum()} CPS controls, only {n_supp} fall in the [0.05, 0.95] common-support band with the trainees.")
print("The propensity score reduces eight covariates to one balancing score (Rosenbaum-Rubin) and exposes the overlap problem.")
Of 429 CPS controls, only 273 fall in the [0.05, 0.95] common-support band with the trainees. The propensity score reduces eight covariates to one balancing score (Rosenbaum-Rubin) and exposes the overlap problem.
3. Matching, and checking balance¶
Matching builds a synthetic control group by pairing each treated unit with the most similar control(s), so the matched groups look alike on $X$ — mimicking the experiment we could not run. We implement two classic schemes from scratch, both estimating the ATT with replacement:
- Nearest-neighbor on the propensity score — pair each treated unit with the control whose $e(X)$ is closest;
- Mahalanobis matching — pair on the full covariate vector using Mahalanobis distance (scale- and correlation-aware).
The decisive diagnostic is covariate balance: the standardized mean differences should collapse toward zero after matching, conventionally inside a ±0.1 band. The Love plot shows the severe pre-matching imbalance collapsing — from 7 of 8 covariates outside the band to 2 under nearest-neighbour propensity matching, though age and married do not quite make it in.
Then we read off the ATT, and the two schemes disagree by a factor of two and a half. The reason is the interesting part, and it runs against the usual moral: Mahalanobis produces the tidier balance table and the worse estimate. Which covariates a scheme chooses to balance turns out to matter more than how well it balances them on average.
ti=np.where(W==1)[0]; ci=np.where(W==0)[0]
# 1-NN propensity matching (with replacement)
mps=ci[np.abs(ps[ti][:,None]-ps[ci][None,:]).argmin(1)]
att_ps=(Y[ti]-Y[mps]).mean()
# Mahalanobis matching (with replacement) on covariates
VI=np.linalg.pinv(np.cov(X[W==0].T))
mmh=ci[cdist(X[ti],X[ci],"mahalanobis",VI=VI).argmin(1)]
att_mah=(Y[ti]-Y[mmh]).mean()
def smd_after(match): return [ (X[ti,j].mean()-X[match,j].mean())/X[ti,j].std() for j in range(len(cov)) ]
before=[smd(j) for j in range(len(cov))]; aps=smd_after(mps); amh=smd_after(mmh)
fig,ax=plt.subplots(1,2,figsize=(13,4.6))
yv=np.arange(len(cov))
ax[0].scatter(before,yv,color=RED,s=55,label="before matching",zorder=3)
ax[0].scatter(aps,yv,color=GREEN,s=55,label="after NN-propensity",zorder=3)
ax[0].scatter(amh,yv,color=BLUE,marker="D",s=40,label="after Mahalanobis",zorder=3)
for y in yv: ax[0].plot([before[y],aps[y]],[y,y],color=GREY,lw=.8,zorder=1)
ax[0].axvline(0,color="k",lw=.7); ax[0].axvline(-0.1,color=GREY,ls=":"); ax[0].axvline(0.1,color=GREY,ls=":")
ax[0].set_yticks(yv); ax[0].set_yticklabels(cov); ax[0].set_xlabel("standardized mean difference"); ax[0].set_title("Love plot: which covariates each scheme balances"); ax[0].legend(fontsize=8)
est={"naive (CPS)":naive,"Mahalanobis":att_mah,"NN-propensity":att_ps}
names=list(est); vals=[est[k] for k in names]
ax[1].barh(names,vals,color=[GREY,BLUE,GREEN]); ax[1].axvline(bench,color=RED,lw=2,ls="--",label=f"experimental truth ${bench:,.0f}"); ax[1].axvline(0,color="k",lw=.6)
lo_,hi_=min(vals+[0]),max(vals+[bench]); pad_=0.03*(hi_-lo_)
ax[1].set_xlim(lo_-7*pad_, hi_+5*pad_) # room for the labels on both sides
for i,v in enumerate(vals):
ax[1].text(v+pad_ if v>=0 else v-pad_, i, f"${v:,.0f}", va="center",
ha="left" if v>=0 else "right", fontsize=8)
ax[1].set_xlabel("estimated ATT (USD)"); ax[1].set_title("ATT vs the experimental benchmark"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
print(f"NN-propensity matching ATT = ${att_ps:,.0f} (experimental truth ${bench:,.0f})")
print(f"Mahalanobis matching ATT = ${att_mah:,.0f}")
bal = pd.DataFrame({"before":before, "after NN-ps":aps, "after Mahalanobis":amh}, index=cov)
print("\nStandardized mean differences:"); print(bal.round(3).to_string())
mx_ps, mx_mh = max(np.abs(aps)), max(np.abs(amh))
out_ps, out_mh = int((np.abs(aps)>0.1).sum()), int((np.abs(amh)>0.1).sum())
print(f"\n max |SMD|: before {max(np.abs(before)):.3f} NN-ps {mx_ps:.3f} Mahalanobis {mx_mh:.3f}")
print(f" outside +/-0.1: before {int((np.abs(before)>0.1).sum())}/8 NN-ps {out_ps}/8 Mahalanobis {out_mh}/8")
print()
print("Read that table carefully, because it does not say what a balance table is usually assumed to say.")
print(f"MAHALANOBIS IS THE BETTER-BALANCED MATCH on every aggregate summary -- max |SMD| {mx_mh:.3f} against")
print(f"{mx_ps:.3f}, and {out_mh} of 8 covariates outside the +/-0.1 band against {out_ps}. And it is the estimator that")
print(f"MISSES, at ${att_mah:,.0f} against a truth of ${bench:,.0f}, while the worse-balanced match lands on ${att_ps:,.0f}.")
print()
print("What separates them is WHICH covariates each one balances. Lagged earnings are the strongly")
print(f"prognostic ones here, and on re75 NN-propensity gets |SMD| {abs(aps[7]):.3f} against Mahalanobis's {abs(amh[7]):.3f};")
print(f"on re74, {abs(aps[6]):.3f} against {abs(amh[6]):.3f}. Mahalanobis buys its tidy aggregate by balancing age and the")
print("demographic dummies almost perfectly -- variables that matter far less for 1978 earnings -- and pays")
print("for it where it counts. Nothing in an overall balance summary would have told you that.")
print()
print("So 'check balance, not just the estimate' is necessary and not sufficient: checked the usual way,")
print("aggregate balance would have SELECTED THE WRONG ESTIMATOR here. Balance has to be judged on the")
print("covariates that predict the outcome, which is a claim about the outcome model, not about the")
print("matching -- and it is exactly the knowledge unconfoundedness already assumes you have.")
NN-propensity matching ATT = $1,792 (experimental truth $1,794)
Mahalanobis matching ATT = $719
Standardized mean differences:
before after NN-ps after Mahalanobis
age -0.310 0.228 0.028
educ 0.055 -0.065 -0.035
black 1.762 0.015 0.000
hispan -0.350 -0.023 0.000
married -0.826 0.152 0.028
nodegree 0.245 0.059 0.000
re74 -0.723 -0.056 0.061
re75 -0.291 0.008 0.138
max |SMD|: before 1.762 NN-ps 0.228 Mahalanobis 0.138
outside +/-0.1: before 7/8 NN-ps 2/8 Mahalanobis 1/8
Read that table carefully, because it does not say what a balance table is usually assumed to say.
MAHALANOBIS IS THE BETTER-BALANCED MATCH on every aggregate summary -- max |SMD| 0.138 against
0.228, and 1 of 8 covariates outside the +/-0.1 band against 2. And it is the estimator that
MISSES, at $719 against a truth of $1,794, while the worse-balanced match lands on $1,792.
What separates them is WHICH covariates each one balances. Lagged earnings are the strongly
prognostic ones here, and on re75 NN-propensity gets |SMD| 0.008 against Mahalanobis's 0.138;
on re74, 0.056 against 0.061. Mahalanobis buys its tidy aggregate by balancing age and the
demographic dummies almost perfectly -- variables that matter far less for 1978 earnings -- and pays
for it where it counts. Nothing in an overall balance summary would have told you that.
So 'check balance, not just the estimate' is necessary and not sufficient: checked the usual way,
aggregate balance would have SELECTED THE WRONG ESTIMATOR here. Balance has to be judged on the
covariates that predict the outcome, which is a claim about the outcome model, not about the
matching -- and it is exactly the knowledge unconfoundedness already assumes you have.
4. Weighting and doubly-robust estimation¶
Matching discards most of the control sample (each treated unit uses one control). Inverse-probability weighting (IPW) instead uses all units, weighting controls by $e(X)/(1-e(X))$ so that the reweighted control group resembles the treated group (Horvitz-Thompson logic). It is consistent if the propensity model is correct.
Doubly-robust AIPW (Augmented IPW) combines the propensity model with an outcome regression $m_0(X)=E[Y\mid X,W=0]$, and has the celebrated double-robustness property: it is consistent if either the propensity model or the outcome model is correct — two chances to get it right. We implement both from scratch (trimming to the common-support band so a near-zero denominator cannot explode a weight), and place every estimator against the benchmark. Three of the four land between USD 1,200 and USD 1,800 — a world away from the naive −USD 635 — while Mahalanobis matching sits well below that band, and every one of them falls at or below the experimental truth rather than straddling it.
keep=(ps>0.05)&(ps<0.95); Wk,Yk,pk=W[keep],Y[keep],ps[keep]
# IPW ATT: treated get weight 1, controls get e/(1-e)
wc=pk/(1-pk)
att_ipw=Yk[Wk==1].mean()-np.sum(wc[Wk==0]*Yk[Wk==0])/np.sum(wc[Wk==0])
# AIPW ATT (doubly robust): outcome model on controls + propensity augmentation
m0=LinearRegression().fit(X[W==0],Y[W==0]).predict(X); N1=(W==1).sum()
att_aipw=(np.sum((W==1)*(Y-m0))-np.sum((W==0)*(ps/(1-ps))*(Y-m0)))/N1
allest={"naive (CPS)":naive,"Mahalanobis\nmatching":att_mah,"IPW":att_ipw,"AIPW\n(doubly robust)":att_aipw,"NN-propensity\nmatching":att_ps}
nm=list(allest); vv=[allest[k] for k in nm]
fig,ax=plt.subplots(figsize=(9,4.8))
cols=[GREY,BLUE,ORANGE,PURP,GREEN]
ax.barh(nm,vv,color=cols); ax.axvline(bench,color=RED,lw=2.5,ls="--"); ax.axvline(0,color="k",lw=.7)
lo_,hi_=min(vv+[0]),max(vv+[bench]); pad_=0.025*(hi_-lo_)
ax.set_xlim(lo_-7*pad_, hi_+7*pad_)
for i,v in enumerate(vv):
ax.text(v+pad_ if v>=0 else v-pad_, i, f"${v:,.0f}", va="center",
ha="left" if v>=0 else "right", fontsize=9)
ax.invert_yaxis(); ax.set_xlabel("estimated ATT (USD)")
ax.set_title(f"Every adjustment method vs the randomized benchmark (dashed = ${bench:,.0f})")
plt.tight_layout(); plt.show()
print(f"IPW ATT = ${att_ipw:,.0f}")
print(f"AIPW ATT = ${att_aipw:,.0f} (doubly robust)")
print(f"Experimental truth ${bench:,.0f}; naive observational ${naive:,.0f}.")
print()
print(f"Every method moves the estimate off the naive ${naive:,.0f} and onto the correct side of zero, which is")
print("the headline. The spread among them is the caveat, and it is wide:")
print(f" NN-propensity ${att_ps:>6,.0f} essentially exact")
print(f" IPW ${att_ipw:>6,.0f} {100*(att_ipw-bench)/bench:>+5.0f}% against the benchmark")
print(f" AIPW ${att_aipw:>6,.0f} {100*(att_aipw-bench)/bench:>+5.0f}%")
print(f" Mahalanobis ${att_mah:>6,.0f} {100*(att_mah-bench)/bench:>+5.0f}%")
print()
print("Three of the four sit between $1,200 and $1,800; Mahalanobis sits well below that band, and all four")
print(f"are at or BELOW the benchmark -- none of them brackets it from above. A practitioner without the answer")
print("key would have four defensible numbers spanning a factor of two and no way to choose between them,")
print("which is LaLonde's original complaint restated. The adjustment is doing real work -- IF")
print("unconfoundedness holds, and nothing here can check that.")
IPW ATT = $1,316 AIPW ATT = $1,226 (doubly robust) Experimental truth $1,794; naive observational $-635. Every method moves the estimate off the naive $-635 and onto the correct side of zero, which is the headline. The spread among them is the caveat, and it is wide: NN-propensity $ 1,792 essentially exact IPW $ 1,316 -27% against the benchmark AIPW $ 1,226 -32% Mahalanobis $ 719 -60% Three of the four sit between $1,200 and $1,800; Mahalanobis sits well below that band, and all four are at or BELOW the benchmark -- none of them brackets it from above. A practitioner without the answer key would have four defensible numbers spanning a factor of two and no way to choose between them, which is LaLonde's original complaint restated. The adjustment is doing real work -- IF unconfoundedness holds, and nothing here can check that.
5. Summary¶
Starting from an observational comparison that said the job-training program reduced earnings by USD 635, covariate adjustment under unconfoundedness moved every estimator onto the correct side of zero, and nearest-neighbour propensity matching landed on the randomized truth of +USD 1,794 almost exactly. The machinery worked in the sense that matters most: the propensity score reduced eight covariates to one balancing score (Rosenbaum-Rubin), and matching and weighting rebuilt a control group out of survey respondents who initially looked nothing like the trainees.
It worked less cleanly than the headline suggests. The four adjusted estimates span USD 719 to USD 1,792, a factor of two and a half, and all four sit at or below the benchmark rather than scattering around it. Without the answer key there would be no way to choose among them.
The honest caveats are the real lesson, and they are why this dataset is famous:
- Balance is testable; unconfoundedness is not. We can show the covariates are balanced after matching, but the whole edifice rests on the untestable claim that no unobserved confounder remains. LaLonde's (1986) original point was that observational methods often fail to recover experimental benchmarks; Dehejia-Wahba (1999) showed careful propensity methods can — and later work showed the answer is sensitive to specification and subsample. The experiment of the previous notebook needed none of this faith.
- Better balance is not a better estimate. This is the finding that surprises. Mahalanobis matching produces the tidier balance table — max |SMD| of 0.138 against nearest-neighbour's 0.228, one covariate outside the ±0.1 band against two — and lands at USD 719, more than a thousand dollars further from the truth. It achieves its tidy aggregate on age and the demographic dummies while balancing lagged earnings less well, and lagged earnings are what predict 1978 earnings. So "always check balance" is necessary and not sufficient: applied to the aggregate, it would have chosen the wrong estimator here. Balance must be judged on the covariates that drive the outcome — which is knowledge about the outcome model, not about the matching.
- Standard errors need care. Matching SEs from a naive bootstrap are invalid (Abadie-Imbens); the R companion uses the correct variance estimators.
Cross-links. Choosing the covariates for the propensity model is the confounder-selection face of the Variable Selection arc; the doubly-robust idea returns, generalized with machine-learning nuisance models and cross-fitting, in the Double/Debiased ML subsection (and already underlies the DML causal forest). Next: Instrumental Variables, which confronts the case unconfoundedness cannot handle — an unobserved confounder — by finding a source of variation in treatment that is as good as randomly assigned.