The idea of reweighting
A targeted programme enrols mostly poor households. To make treated and control groups comparable, IPW up-weights people who received an unlikely treatment: a well-off household that enrolled represents many similar households that did not. Weights:
| Target | Treated weight | Control weight |
|---|---|---|
| ATE | 1 / e(X) | 1 / (1 − e(X)) |
| ATT | 1 | e(X) / (1 − e(X)) |
In the weighted pseudo-population, treatment is independent of X, so a weighted difference in means is causal under unconfoundedness and overlap.
!pip install -q numpy pandas statsmodels matplotlib seaborn causaldata
import numpy as np, pandas as pd
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt, seaborn as sns
rng = np.random.default_rng(42)
sns.set_theme(style="whitegrid")
A targeted tutoring programme
A supplemental tutoring programme targets disadvantaged pupils: a higher poverty index raises the chance of enrolment, and poverty also lowers test scores. The true effect is +10 points.
N = 2_000
poverty = rng.normal(0, 1, N)
p_true = 1 / (1 + np.exp(-(-0.5 + 1.2 * poverty)))
treat = rng.binomial(1, p_true)
score = 50 + 10 * treat - 15 * poverty + rng.normal(0, 5, N)
df = pd.DataFrame({"poverty": poverty, "treat": treat, "score": score})
df["ps"] = smf.logit("treat ~ poverty", data=df).fit(disp=0).predict(df)
w_ate = np.where(df.treat == 1, 1 / df.ps, 1 / (1 - df.ps))
def weighted_diff(y, d, w):
"""Hájek (normalised) estimator: weighted mean of treated minus weighted mean of controls."""
return np.average(y[d == 1], weights=w[d == 1]) - np.average(y[d == 0], weights=w[d == 0])
naive = df.score[df.treat == 1].mean() - df.score[df.treat == 0].mean()
ipw = weighted_diff(df.score.to_numpy(), df.treat.to_numpy(), w_ate)
wls = smf.wls("score ~ treat", data=df, weights=w_ate).fit(cov_type="HC1")
print(f"Naive: {naive:6.2f} | IPW (Hájek): {ipw:6.2f} | WLS: {wls.params['treat']:6.2f} (robust SE {wls.bse['treat']:.2f}) | truth 10")
The naive estimate is negative: tutored pupils are poorer, and poverty depresses scores. IPW recovers the truth. Weighted least squares with a treatment dummy gives the same Hájek estimate plus a robust standard error (conservative, because it ignores that the score was estimated).
Weighted balance check. After weighting, the covariate means should match:
for g in (1, 0):
m = df.treat == g
print(f"treat={g}: raw mean poverty {df.poverty[m].mean():+.3f} | weighted {np.average(df.poverty[m], weights=w_ate[m]):+.3f}")
Extreme weights and trimming
When e(X) is close to 0 or 1, weights explode: one observation can dominate the estimate. Inspect the distribution of weights and the effective sample size (ESS = (Σw)² / Σw²).
def ess(w):
return w.sum() ** 2 / (w ** 2).sum()
# A harsher DGP: strong targeting -> poor overlap
p_hard = 1 / (1 + np.exp(-(-0.5 + 3.0 * poverty)))
treat_h = rng.binomial(1, p_hard)
score_h = 50 + 10 * treat_h - 15 * poverty + rng.normal(0, 5, N)
h = pd.DataFrame({"poverty": poverty, "treat": treat_h, "score": score_h})
h["ps"] = smf.logit("treat ~ poverty", data=h).fit(disp=0).predict(h)
w_h = np.where(h.treat == 1, 1 / h.ps, 1 / (1 - h.ps))
print(f"Max weight {w_h.max():,.0f} | ESS {ess(w_h):,.0f} of {N:,}")
keep = (h.ps > 0.05) & (h.ps < 0.95)
est_all = weighted_diff(h.score.to_numpy(), h.treat.to_numpy(), w_h)
est_trim = weighted_diff(h.score[keep].to_numpy(), h.treat[keep].to_numpy(), w_h[keep])
print(f"IPW all: {est_all:.2f} | IPW trimmed to [0.05, 0.95]: {est_trim:.2f} (kept {keep.sum():,})")
sns.histplot(w_h, bins=80, log_scale=True); plt.title("IPW weights (log scale)"); plt.show()
Trimming trades a little bias (the estimand becomes the overlap population) for much lower variance. In Module 7, doubly robust estimators reduce the damage from imperfect weights further, by combining weighting with an outcome model.
ATT weights on the LaLonde data
from causaldata import nsw_mixtape, cps_mixtape
nsw = nsw_mixtape.load_pandas().data
obs = pd.concat([nsw[nsw.treat == 1], cps_mixtape.load_pandas().data], ignore_index=True)
obs["u74"] = (obs.re74 == 0).astype(int); obs["u75"] = (obs.re75 == 0).astype(int)
obs["ps"] = smf.logit("treat ~ age + I(age**2) + I(age**3) + educ + I(educ**2) + marr + nodegree"
" + black + hisp + re74 + re75 + u74 + u75 + educ:re74", data=obs).fit(disp=0).predict(obs)
w_att = np.where(obs.treat == 1, 1.0, obs.ps / (1 - obs.ps))
att = weighted_diff(obs.re78.to_numpy(), obs.treat.to_numpy(), w_att)
print(f"IPW ATT: ${att:,.0f} (experimental $1,794)")
print(f"Control ESS: {ess(w_att[obs.treat == 0]):,.0f} of {(obs.treat == 0).sum():,} CPS respondents")
The ESS shows that the 16,000 CPS respondents are worth only a few hundred "effective" controls: the ones who resemble trainees. That is the information really available for the comparison.
Exercise
Write a function ipw_bootstrap(data, reps=200) that re-estimates the propensity score and the IPW ATT in
each bootstrap sample of the LaLonde data and returns a 95% percentile interval. Does it contain $1,794?
Show a solution
F = ("treat ~ age + I(age**2) + I(age**3) + educ + I(educ**2) + marr + nodegree"
" + black + hisp + re74 + re75 + u74 + u75 + educ:re74")
def ipw_bootstrap(data, reps=200):
out = []
for s in range(reps):
b = data.sample(frac=1, replace=True, random_state=s)
ps = smf.logit(F, data=b).fit(disp=0).predict(b)
w = np.where(b.treat == 1, 1.0, ps / (1 - ps))
out.append(weighted_diff(b.re78.to_numpy(), b.treat.to_numpy(), w))
return np.percentile(out, [2.5, 97.5])
print(ipw_bootstrap(obs, reps=100).round(0))
Checklist and recap
- I can write ATE and ATT weights and explain the pseudo-population.
- I check weighted balance, the largest weights and the effective sample size.
- I trim when overlap is poor and say how the estimand changes.