6.1 Selection on Observables and the Propensity Score

When we cannot randomise

Most policy data are observational: people chose (or were chosen for) the programme. If we observe every variable that drives both participation and outcomes, we can still compare like with like. Two assumptions make that work:

1. Conditional independence (unconfoundedness): given covariates X, treatment is as good as random: (Y(1), Y(0)) ⫫ D | X. In DAG terms, X satisfies the backdoor criterion (Module 5).
2. Overlap (common support): every type of person has some chance of being treated and untreated: 0 < P(D = 1 | X) < 1.

The propensity score e(X) = P(D = 1 | X) compresses all covariates into one number. Rosenbaum and Rubin showed that if unconfoundedness holds given X, it also holds given e(X), so we can match or weight on a single dimension.

!pip install -q numpy pandas statsmodels scikit-learn matplotlib seaborn causaldata
import numpy as np, pandas as pd
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt, seaborn as sns
from causaldata import nsw_mixtape, cps_mixtape
sns.set_theme(style="whitegrid")

LaLonde's challenge: experimental vs observational

LaLonde (1986) took the NSW experiment (Module 4: about $1,794 effect on 1978 earnings) and replaced the randomised control group with people from a national survey, the CPS. Could standard methods recover the experimental answer? Mostly they could not, and that finding launched modern causal inference.

nsw = nsw_mixtape.load_pandas().data
cps = cps_mixtape.load_pandas().data
obs = pd.concat([nsw[nsw.treat == 1], cps], ignore_index=True)     # NSW trainees + CPS non-trainees
obs["u74"] = (obs.re74 == 0).astype(int)                            # unemployed in 1974
obs["u75"] = (obs.re75 == 0).astype(int)
EXPERIMENTAL = 1_794

naive = obs.re78[obs.treat == 1].mean() - obs.re78[obs.treat == 0].mean()
ols = smf.ols("re78 ~ treat + age + educ + black + hisp + marr + nodegree + re74 + re75",
              data=obs).fit(cov_type="HC1").params["treat"]
print(f"Treated: {obs.treat.sum()} | CPS comparisons: {(obs.treat == 0).sum():,}")
print(f"Experimental benchmark: ${EXPERIMENTAL:,}")
print(f"Naive difference:       ${naive:,.0f}")
print(f"OLS with controls:      ${ols:,.0f}")
print(obs.groupby("treat")[["age", "educ", "black", "marr", "nodegree", "re74", "re75"]].mean().round(2).T)

The trainees were young, mostly Black, unmarried, without diplomas and with almost no earnings in 1974–75. The typical CPS respondent is nothing like them. The naive comparison is hugely negative and even regression, which extrapolates linearly across very different people, misses badly.

Estimating the propensity score

We fit a flexible logit, following Dehejia and Wahba (1999): squares, a cube in age, unemployment dummies and an interaction. The goal is not to "explain" participation but to balance covariates.

ps_formula = ("treat ~ age + I(age**2) + I(age**3) + educ + I(educ**2) + marr + nodegree"
              " + black + hisp + re74 + re75 + u74 + u75 + educ:re74")
ps_model = smf.logit(ps_formula, data=obs).fit(disp=0)
obs["ps"] = ps_model.predict(obs)
print(obs.groupby("treat")["ps"].describe().round(4))

Overlap and common support

Overlap is checked, not assumed. Plot the propensity score by group; where one group has no counterparts in the other, there is nothing to compare with, and any estimate there is pure extrapolation.

fig, axes = plt.subplots(1, 2, figsize=(13, 4))
sns.histplot(data=obs, x="ps", hue="treat", bins=50, stat="density", common_norm=False, ax=axes[0])
axes[0].set(title="Propensity score by group", yscale="log")
obs["logit_ps"] = np.log(obs.ps / (1 - obs.ps))
sns.kdeplot(data=obs, x="logit_ps", hue="treat", common_norm=False, ax=axes[1])
axes[1].set(title="Log-odds scale shows the overlap region")
plt.tight_layout(); plt.show()

Most CPS respondents have a propensity score near zero: they look nothing like trainees. A simple and widely used rule (Crump et al., 2009) keeps observations with scores between 0.1 and 0.9, which targets the effect for the overlap population. Here, where the target is the ATT, we only drop controls below the lowest treated score and treated units above the highest control score.

lo, hi = obs.loc[obs.treat == 1, "ps"].min(), obs.loc[obs.treat == 0, "ps"].max()
support = obs[(obs.ps >= lo) & (obs.ps <= hi)]
print(f"Common support: [{lo:.5f}, {hi:.3f}]")
print(f"Kept {len(support):,} of {len(obs):,} rows; treated kept: {support.treat.sum()} of {obs.treat.sum()}")

Exercise

simple = smf.logit("treat ~ age + educ + black + hisp + marr + nodegree + re74 + re75", data=obs).fit(disp=0)
obs["ps_simple"] = simple.predict(obs)
lo_s, hi_s = obs.loc[obs.treat == 1, "ps_simple"].min(), obs.loc[obs.treat == 0, "ps_simple"].max()
print(f"Treated outside support (simple spec): {((obs.treat == 1) & (obs.ps_simple > hi_s)).sum()}")
print(obs.groupby("treat")["ps_simple"].describe().round(4))
# The specification matters: the score must balance covariates. Lesson 6.2 shows how to check that.

Checklist and recap

• I can state unconfoundedness and overlap, and link unconfoundedness to the backdoor criterion.
• I can estimate a propensity score and plot overlap.
• I trim to common support and report how the estimand changes.