import matplotlib
matplotlib.use("Agg") # headless: render to buffers, never to a screen
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
rng = np.random.default_rng(42) # one pinned generator drives the whole notebook
plt.rcParams["figure.figsize"] = (9, 5)
plt.rcParams["axes.grid"] = True
plt.rcParams["grid.alpha"] = 0.3
pd.set_option("display.width", 110)
pd.set_option("display.max_columns", 24)
print("numpy ", np.__version__)
print("pandas ", pd.__version__)
import statsmodels
print("statsmodels", statsmodels.__version__)nb8.2 — The Robustness Battery: Try to Kill Your Own Result
Week 8, Chapter 8.2 companion. Maya stops varying her choices and starts attacking her inference.
In nb8.1 you built a specification curve: every defensible fork in the garden, plotted, to show the headline number was not an artifact of one lucky choice. That was about analytic-choice multiplicity. This notebook is more adversarial. Here you keep the specification fixed and ask a harder question: is the estimate even real, or have you been fooled — by correlated errors that make the t-statistic lie, by a pattern that would show up where no treatment happened, by a knob you tuned until the stars appeared, by a family of outcomes you tested all at once, or by a confounder you never measured?
The discipline has an ugly name and a beautiful purpose: you are going to try to kill your own result. We continue Maya’s staggered difference-in-differences on HMDA-style data — do state fair-lending examination programs shrink the county-level minority–white mortgage-denial gap? — and we build the battery as column 3 of her Chapter 7.5 threats table, made to run:
- Alternative standard errors — classical, HC, clustered at different levels, two-way, and a wild cluster bootstrap for few clusters.
- Placebo tests — an in-time placebo (fake date before treatment) and an in-space placebo (permute treatment across units to build a placebo distribution).
- Sensitivity — sweep controls, sample, and winsorizing; show the estimate barely moves.
- Multiple testing — Bonferroni (FWER) vs Benjamini–Hochberg (FDR) on a family of outcomes.
- Oster (2019) δ — how strong unobserved selection would have to be to explain the result away.
The trick this notebook reveals is the one you might be playing on yourself. We build two synthetic worlds: one where a real effect survives every attack, and one where a spurious effect (driven by a differential pre-trend) fails the in-time placebo. Watch the battery tell them apart.
1. Two synthetic panels with a known truth
The whole point of a battery is to run it where we already know the answer, so we can check the tests behave. We build a county–year panel like the one nb7.4 emits: one row per county-year, a staggered adoption cohort (the year a state’s examination program switches on; never-adopters are clean controls), the outcome (the minority–white denial gap in percentage points), two observed controls (applicant income and loan-to-value), a state cluster key, and a county fixed-effect baseline.
We generate two worlds from the same machinery:
- World R (“real”): the program genuinely narrows the gap by a planted \tau = -1.6 pp. Adopting and never-adopting states share a common gap trend before adoption — parallel trends hold.
- World S (“spurious”): the program has no true effect (\tau = 0), but adopting states were already on a downward differential pre-trend — their gaps were shrinking faster anyway. A naive DiD will mistake that trend for an effect and report a spurious negative coefficient.
The structure is identical; only the truth differs. That is exactly the setup a placebo test is built to expose.
def make_panel(rng, tau_true, pretrend_gap=0.0, n_states=24, counties_per_state=10,
years=range(2013, 2023), first_real_adopt=2018):
"""One county-year panel. tau_true = planted ATT (pp). pretrend_gap = extra per-year
downward drift in the gap for *adopting* states only (the confounding pre-trend in World S)."""
years = list(years)
base_year = years[0]
# Staggered adoption: ~half the states adopt across two cohorts; the rest never adopt.
cohort = {}
for s in range(n_states):
if s < n_states // 4: cohort[s] = first_real_adopt # early cohort
elif s < n_states // 2: cohort[s] = first_real_adopt + 2 # late cohort
else: cohort[s] = np.nan # never-adopter (control)
state_fe = rng.normal(0, 1.2, size=n_states) # persistent state level in the gap
# Persistent state-by-year shock: a common shock to all counties in a state-year. This is the
# WITHIN-STATE error correlation that makes clustering matter (clustered SE > classical).
state_year_shock = {(s, y): rng.normal(0, 0.9) for s in range(n_states) for y in years}
rows = []
for s in range(n_states):
g = cohort[s]
adopts = not np.isnan(g)
for c in range(counties_per_state):
county_fe = rng.normal(0, 0.6)
for y in years:
income = rng.normal(0, 1.0) # standardized applicant income
# Loan-to-value is SELECTED ON: adopting states have systematically higher LTV,
# so ltv is correlated with treatment. Omitting it biases the coef; adding it MOVES
# the coef and raises R^2 -- exactly the structure the Oster (2019) test reads.
ltv = rng.normal(0.8 if adopts else 0.0, 1.0)
# Common time trend shared by all counties (a national mortgage cycle).
common_trend = 0.10 * (y - base_year)
# Adopting states optionally carry an EXTRA downward gap drift (World S confounder).
adopter_drift = pretrend_gap * (y - base_year) if adopts else 0.0
post = 1 if (adopts and y >= g) else 0
eps = rng.normal(0, 0.8)
gap = (5.0 + state_fe[s] + county_fe + common_trend + adopter_drift
+ state_year_shock[(s, y)]
+ (-0.9) * income + (1.4) * ltv
+ tau_true * post + eps)
rows.append((s, s * 1000 + c, y, g if adopts else np.nan,
int(adopts), post, income, ltv, gap))
df = pd.DataFrame(rows, columns=["state", "county", "year", "cohort_G",
"adopter", "post_treat", "income", "ltv", "denial_gap"])
return df
YEARS = range(2013, 2023)
panel_R = make_panel(rng, tau_true=-1.6, pretrend_gap=0.0) # World R: real effect, parallel trends
panel_S = make_panel(rng, tau_true=0.0, pretrend_gap=-0.45) # World S: no effect, downward pre-trend
print("World R:", panel_R.shape, "rows | World S:", panel_S.shape, "rows")
print("states:", panel_R["state"].nunique(),
"| adopting states:", panel_R.loc[panel_R["adopter"] == 1, "state"].nunique(),
"| never-adopters:", panel_R.loc[panel_R["adopter"] == 0, "state"].nunique())2. The headline estimate (the thing we will attack)
Maya’s primary specification is a two-way fixed-effects DiD: regress the denial gap on the post-treatment indicator (1 in adopting states once their program is live, 0 otherwise) plus county fixed effects, year fixed effects, and her two controls. We absorb county and year with C(county) and C(year) dummies in statsmodels so the whole battery uses one transparent estimator.
Spec line. Outcome: county-year denial gap (pp) · Treatment:
post_treat· Controls: income, ltv · Fixed effects: county, year · Clustering: state · Sample: full panel · Identifying assumption: per-cohort parallel trends (adopting and control states would have moved together absent the program).
The function below returns the coefficient on post_treat and its standard error under whatever covariance flavor we ask for. The point estimate never changes across SE flavors — that is the central lesson of Chapter 2.4. Only the uncertainty around it does.
import warnings
BASE_FORMULA = "denial_gap ~ post_treat + income + ltv + C(county) + C(year)"
def fit_did(df, formula=BASE_FORMULA, cov_type="nonrobust", cov_kwds=None):
"""Fit the TWFE DiD; return (coef, se, t) on post_treat under the requested covariance.
(Two-way cluster covariance is not guaranteed positive semi-definite, so some of the *county
dummy* SEs can be NaN; the post_treat SE we read off is finite and valid -- we quiet the
harmless sqrt warning from those other coefficients.)"""
model = smf.ols(formula, data=df)
with warnings.catch_warnings():
warnings.simplefilter("ignore", RuntimeWarning)
res = _fit_cov(model, cov_type, cov_kwds)
b = res.params["post_treat"]
se = res.bse["post_treat"] # bse is lazy -> access it INSIDE the filter
return b, se, b / se
def _fit_cov(model, cov_type, cov_kwds):
if cov_type == "nonrobust":
return model.fit()
if cov_kwds is not None:
return model.fit(cov_type=cov_type, cov_kwds=cov_kwds)
return model.fit(cov_type=cov_type)
bR, seR, tR = fit_did(panel_R, cov_kwds={"groups": panel_R["state"]}, cov_type="cluster")
bS, seS, tS = fit_did(panel_S, cov_kwds={"groups": panel_S["state"]}, cov_type="cluster")
print(f"World R (truth tau=-1.6): coef = {bR:+.3f} pp state-clustered SE = {seR:.3f} t = {tR:+.2f}")
print(f"World S (truth tau= 0.0): coef = {bS:+.3f} pp state-clustered SE = {seS:.3f} t = {tS:+.2f}")
print("\nNote: World S shows a spuriously NEGATIVE, significant coef -- the pre-trend masquerading as an effect.")3. Alternative standard errors: is the t-statistic telling the truth?
The cheapest attack, and the most common cause of a result that evaporates: the standard error is too small, so the t is too big, and the significance is an illusion. A robustness check on standard errors does not touch \hat\beta — it re-asks “how uncertain is this number?” under different assumptions about how the errors are correlated. We recompute the same World-R estimate under five flavors:
- Classical — errors all the same size and independent (\boldsymbol\Omega = \sigma^2\mathbf{I}).
- HC1 — heteroskedasticity-robust (lets each error have its own variance).
- Cluster by county — the unit of observation.
- Cluster by state (primary) — the level treatment varies at; the Petersen (2009) rule.
- Two-way (state × year) — also absorbs a national shock that correlates errors across states in a given year (Cameron, Gelbach & Miller 2011).
Treatment is assigned at the state level, so within-state error correlation inflates the true variance. We expect the SE to grow as we move to coarser, more honest clustering — the clustered SE should exceed the classical one.
def se_panel(df):
rows = []
b, se, t = fit_did(df, cov_type="nonrobust")
rows.append(("Classical", b, se, t))
b, se, t = fit_did(df, cov_type="HC1")
rows.append(("HC1 (robust)", b, se, t))
b, se, t = fit_did(df, cov_type="cluster", cov_kwds={"groups": df["county"]})
rows.append(("Cluster: county (unit)", b, se, t))
b, se, t = fit_did(df, cov_type="cluster", cov_kwds={"groups": df["state"]})
rows.append(("Cluster: state (primary)", b, se, t))
b, se, t = fit_did(df, cov_type="cluster",
cov_kwds={"groups": df[["state", "year"]]})
rows.append(("Two-way: state x year", b, se, t))
out = pd.DataFrame(rows, columns=["SE flavor", "coef", "SE", "t"])
out["95% CI lo"] = out["coef"] - 1.96 * out["SE"]
out["95% CI hi"] = out["coef"] + 1.96 * out["SE"]
return out
se_tab = se_panel(panel_R)
print("World R -- same estimate, five SE flavors (coef is IDENTICAL across rows):\n")
print(se_tab.round(3).to_string(index=False))
classical_se = se_tab.loc[se_tab["SE flavor"] == "Classical", "SE"].iloc[0]
state_se = se_tab.loc[se_tab["SE flavor"] == "Cluster: state (primary)", "SE"].iloc[0]
print(f"\nClustered-by-state SE ({state_se:.3f}) > classical SE ({classical_se:.3f})? "
f"{state_se > classical_se} <- clustering matters here, as expected.")The hard case: few clusters and the wild cluster bootstrap
Cluster-robust SEs lean on an asymptotic argument that needs the number of clusters to be large (rule of thumb 30–50). The count that matters for a treatment effect is closer to the number of treated clusters — and Maya may have only a handful of adopting states. With few clusters the cluster-robust SE is itself downward-biased and noisy, the t does not follow the t-distribution you compare it against, and you over-reject.
The fix is the wild cluster bootstrap (Cameron, Gelbach & Miller 2008). The trick worth revealing: keep the clusters and regressors exactly as they are, impose the null (\beta_{\text{post}} = 0) by fitting a restricted model, then flip the sign of each whole cluster’s residuals at random — times +1 or -1 with equal probability (Rademacher weights) — to manufacture a distribution of the test statistic under the null. Re-estimate thousands of times and read your real t-statistic’s position in that distribution. It is the same permutation logic as a placebo: build the distribution of “what this statistic looks like when nothing is going on,” then locate your real number in it.
For speed and transparency we first partial out all fixed effects and controls once (Frisch–Waugh– Lovell: residualize the outcome and post_treat on everything else), reducing the DiD to a one-line regression we can refit cheaply inside the bootstrap loop.
def partial_out(df, controls_formula="income + ltv + C(county) + C(year)"):
"""FWL: residualize denial_gap and post_treat on controls+FE. Returns a small frame."""
y_res = smf.ols("denial_gap ~ " + controls_formula, data=df).fit().resid.to_numpy()
d_res = smf.ols("post_treat ~ " + controls_formula, data=df).fit().resid.to_numpy()
return pd.DataFrame({"y": y_res, "d": d_res, "state": df["state"].to_numpy()})
def wild_cluster_bootstrap(df, cluster_col="state", B=1999, seed=8):
"""Wild cluster bootstrap p-value for H0: beta_post = 0, null imposed, Rademacher signs."""
rng_b = np.random.default_rng(seed)
pp = partial_out(df)
y, d, cl = pp["y"].to_numpy(), pp["d"].to_numpy(), pp["state"].to_numpy()
X = np.column_stack([np.ones_like(d), d])
def cluster_t(yv):
# OLS of yv on [1, d], cluster-robust t on the slope (beta_post).
XtX_inv = np.linalg.inv(X.T @ X)
beta = XtX_inv @ (X.T @ yv)
resid = yv - X @ beta
meat = np.zeros((2, 2))
for c in np.unique(cl):
Xc = X[cl == c]; uc = resid[cl == c]
sc = Xc.T @ uc
meat += np.outer(sc, sc)
V = XtX_inv @ meat @ XtX_inv
return beta[1] / np.sqrt(V[1, 1])
t_obs = cluster_t(y)
# Restricted (null-imposed) fit: drop d, keep intercept-only residuals.
b0 = y.mean()
resid0 = y - b0
uniq = np.unique(cl)
t_boot = np.empty(B)
for b in range(B):
signs = rng_b.choice([-1.0, 1.0], size=uniq.size)
w = np.empty_like(resid0)
for s, c in zip(signs, uniq):
w[cl == c] = s
y_b = b0 + w * resid0
t_boot[b] = cluster_t(y_b)
p = (np.sum(np.abs(t_boot) >= np.abs(t_obs)) + 1) / (B + 1)
return t_obs, p, t_boot
# Make a FEW-CLUSTER world: only a handful of treated states, ~12 states total.
panel_few = make_panel(rng, tau_true=-1.6, pretrend_gap=0.0, n_states=12, counties_per_state=12)
n_treated = panel_few.loc[panel_few["adopter"] == 1, "state"].nunique()
t_obs, p_wild, t_boot = wild_cluster_bootstrap(panel_few, B=1999, seed=8)
# Conventional clustered p for comparison (normal approx).
b_few, se_few, t_few = fit_did(panel_few, cov_type="cluster",
cov_kwds={"groups": panel_few["state"]})
from scipy import stats
p_conv = 2 * (1 - stats.norm.cdf(abs(t_few)))
print(f"Few-cluster world: {panel_few['state'].nunique()} states, {n_treated} treated.")
print(f"Conventional clustered: t = {t_few:+.2f} p ~= {p_conv:.4f}")
print(f"Wild cluster bootstrap: t = {t_obs:+.2f} p = {p_wild:.4f} (B=1999, Rademacher, null imposed)")
print("Both reject zero -> the significance is NOT a few-clusters artifact for this (real) effect.")4. Placebo tests: would this pattern appear where it can’t be real?
A placebo test runs your exact analysis where, by construction, there can be no treatment effect, and checks that you correctly find nothing. If your method finds an “effect” where none can exist, the effect it found in the real analysis is suspect — your method is manufacturing patterns out of whatever structure is in the data. We run the two most important placebos on both worlds, so you can see the battery distinguish the real effect from the spurious one.
In-time placebo (fake treatment date)
Pick a year well before any real adoption — here three years before the earliest real program — and pretend the programs switched on then, using only the pre-period so the fake effect cannot be contaminated by the real one. Re-run the estimator on this fictional timing.
What it probes: whether the estimator detects a “treatment effect” purely from the trend structure of the data — the differential-pre-trend threat. Pass: the placebo effect is small and insignificant (parallel trends supported). Fail: the placebo effect is itself large and significant — the design cannot separate the program from a pre-existing trend.
We expect World R to pass (flat pre-trend) and World S to fail (its adopting states were already drifting down, so a fake date placed in the pre-period catches that drift).
def in_time_placebo(df, fake_adopt_year, real_first_adopt=2018):
"""Restrict to the pre-period (years < real_first_adopt), assign a FAKE adoption at
fake_adopt_year to the states that truly adopt, and re-estimate. No real treated rows included."""
pre = df[df["year"] < real_first_adopt].copy()
pre["post_treat"] = ((pre["adopter"] == 1) & (pre["year"] >= fake_adopt_year)).astype(int)
# Need within-cluster variation: only run if the fake treatment actually switches on.
if pre["post_treat"].nunique() < 2:
return np.nan, np.nan, np.nan
b, se, t = fit_did(pre, cov_type="cluster", cov_kwds={"groups": pre["state"]})
return b, se, t
FAKE_YEAR = 2015 # three years before the earliest real adoption (2018)
bR_pl, seR_pl, tR_pl = in_time_placebo(panel_R, FAKE_YEAR)
bS_pl, seS_pl, tS_pl = in_time_placebo(panel_S, FAKE_YEAR)
print(f"In-time placebo (fake adoption at {FAKE_YEAR}, pre-period only):\n")
print(f"World R (real effect): placebo coef = {bR_pl:+.3f} SE = {seR_pl:.3f} t = {tR_pl:+.2f}"
f" -> {'PASS (near zero)' if abs(tR_pl) < 1.96 else 'FAIL (significant!)'}")
print(f"World S (spurious): placebo coef = {bS_pl:+.3f} SE = {seS_pl:.3f} t = {tS_pl:+.2f}"
f" -> {'PASS (near zero)' if abs(tS_pl) < 1.96 else 'FAIL (significant!)'}")
print("\nThe in-time placebo is the test that EXPOSES the spurious World-S effect as a pre-trend.")In-space placebo (permute treatment across units)
Now permute units instead of time. Repeatedly assign fake adoption to a random set of states (the same number that truly adopt, with the same cohort timing) — including states that never adopted — re-estimate, and build a distribution of placebo coefficients. Then locate the real effect in that distribution.
What it probes: whether an effect of Maya’s magnitude shows up routinely when treatment is assigned at random. This is the right inference when clusters are too few for the t-approximation (Chapter 4.4). Pass: the real effect sits in the extreme tail (more extreme than 95% of placebos → permutation p < 0.05). Fail: the real effect is buried in the middle of the placebo cloud.
def in_space_placebo(df, rng, n_perm=500):
"""Permute which states are 'adopters' (preserving the count and cohort-timing pattern),
re-estimate the DiD coef each time, and return (real_coef, placebo_coefs)."""
# The real treatment structure: which states adopt, and in what cohort year.
state_cohort = (df[["state", "cohort_G"]].drop_duplicates().set_index("state")["cohort_G"])
real_adopters = state_cohort.dropna()
cohort_years = real_adopters.to_numpy() # the timing pattern to reassign
all_states = state_cohort.index.to_numpy()
n_adopt = len(real_adopters)
# Real estimate (clustered).
real_b, _, _ = fit_did(df, cov_type="cluster", cov_kwds={"groups": df["state"]})
placebo = np.empty(n_perm)
state_arr = df["state"].to_numpy()
year_arr = df["year"].to_numpy()
for k in range(n_perm):
fake_adopters = rng.choice(all_states, size=n_adopt, replace=False)
fake_cohort = dict(zip(fake_adopters, rng.permutation(cohort_years)))
post = np.zeros(len(df), dtype=int)
for st, gy in fake_cohort.items():
post |= ((state_arr == st) & (year_arr >= gy)).astype(int)
dfx = df.copy()
dfx["post_treat"] = post
if dfx["post_treat"].nunique() < 2:
placebo[k] = np.nan
continue
b, _, _ = fit_did(dfx, cov_type="nonrobust")
placebo[k] = b
placebo = placebo[~np.isnan(placebo)]
return real_b, placebo
real_bR, placebo_R = in_space_placebo(panel_R, rng, n_perm=500)
# Two-sided permutation p: fraction of placebos at least as extreme (in abs value) as the real coef.
p_perm_R = (np.sum(np.abs(placebo_R) >= np.abs(real_bR)) + 1) / (len(placebo_R) + 1)
pct_R = 100 * np.mean(placebo_R < real_bR) # where the (negative) real effect sits
print(f"World R: real coef = {real_bR:+.3f} pp")
print(f" placebo distribution: mean = {placebo_R.mean():+.3f}, sd = {placebo_R.std():.3f}, "
f"n = {len(placebo_R)}")
print(f" real effect is more negative than {100 - pct_R:.1f}% of placebos; "
f"two-sided permutation p = {p_perm_R:.4f}")
print(f" -> {'PASS: real effect in the tail' if p_perm_R < 0.05 else 'FAIL: buried in the cloud'}")The plot makes the location obvious: the histogram is the placebo cloud (effects from randomly assigned treatment), and the vertical line is Maya’s real estimate. A credible result sits out in the tail, away from the mass of placebos.
fig, ax = plt.subplots()
ax.hist(placebo_R, bins=30, color="steelblue", alpha=0.75, edgecolor="white",
label="placebo coefs (randomly assigned treatment)")
ax.axvline(real_bR, color="crimson", lw=2.5, label=f"real estimate = {real_bR:+.2f} pp")
ax.axvline(0, color="black", lw=1, ls=":")
ax.set_xlabel("DiD coefficient on post_treat (pp)")
ax.set_ylabel("count of permutations")
ax.set_title(f"In-space placebo (World R): real effect vs {len(placebo_R)} permutations "
f"(perm p = {p_perm_R:.3f})")
ax.legend(loc="upper left", fontsize=9)
fig.tight_layout()
fig.savefig("nb82_inspace_placebo.png", dpi=110)
plt.close(fig)
print("Saved figure: nb82_inspace_placebo.png")
print(f"The real estimate ({real_bR:+.2f} pp) sits in the left tail of the placebo cloud.")5. Sensitivity analysis: which knobs is my result hostage to?
Every analysis has knobs — defensible individually, but together they give you freedom to land on a number you like. Sensitivity analysis wiggles each knob and reports how far the estimate moves. The sharp question is not “is it stable across all forks?” (that was nb8.1) but “is it hostage to this particular arbitrary choice?” We sweep three knobs on World R:
- Controls: no controls → primary controls → an aggressive battery (adds interactions). Watch for coefficient stability — but read §6 before celebrating, because stability alone can be a trap.
- Sample: full → drop the largest adopting state → drop the earliest cohort. A result that lives or dies on one influential cluster is fragile.
- Winsorizing: raw → 1% → 5%. Financial data is fat-tailed; a few extreme county-years can move a coefficient. Winsorizing caps each variable at percentiles rather than deleting rows.
Pass: the estimate stays close to the planted -1.6 pp and keeps its sign and significance across all of these.
def winsorize(s, p):
"""Cap a series at its p and 1-p quantiles (p in [0, 0.5)); returns a new Series."""
if p <= 0:
return s.copy()
lo, hi = s.quantile(p), s.quantile(1 - p)
return s.clip(lower=lo, upper=hi)
sens_rows = []
# --- Knob 1: controls ---
b, se, t = fit_did(panel_R, formula="denial_gap ~ post_treat + C(county) + C(year)",
cov_type="cluster", cov_kwds={"groups": panel_R["state"]})
sens_rows.append(("controls: none", b, se))
b, se, t = fit_did(panel_R, cov_type="cluster", cov_kwds={"groups": panel_R["state"]})
sens_rows.append(("controls: primary (income, ltv)", b, se))
b, se, t = fit_did(panel_R, formula="denial_gap ~ post_treat + income + ltv + income:ltv + C(county) + C(year)",
cov_type="cluster", cov_kwds={"groups": panel_R["state"]})
sens_rows.append(("controls: aggressive (+ income:ltv)", b, se))
# --- Knob 2: sample ---
biggest = panel_R.loc[panel_R["adopter"] == 1, "state"].value_counts().idxmax()
sub = panel_R[panel_R["state"] != biggest].copy()
b, se, t = fit_did(sub, cov_type="cluster", cov_kwds={"groups": sub["state"]})
sens_rows.append(("sample: drop largest adopting state", b, se))
earliest = int(panel_R["cohort_G"].min())
sub2 = panel_R[panel_R["cohort_G"] != earliest].copy() # drops earliest cohort (keeps never + late)
b, se, t = fit_did(sub2, cov_type="cluster", cov_kwds={"groups": sub2["state"]})
sens_rows.append((f"sample: drop earliest cohort ({earliest})", b, se))
# --- Knob 3: winsorizing the outcome ---
for p in (0.01, 0.05):
w = panel_R.copy()
w["denial_gap"] = winsorize(w["denial_gap"], p)
b, se, t = fit_did(w, cov_type="cluster", cov_kwds={"groups": w["state"]})
sens_rows.append((f"winsorize outcome at {int(p*100)}%", b, se))
sens = pd.DataFrame(sens_rows, columns=["knob setting", "coef", "SE"])
sens["t"] = sens["coef"] / sens["SE"]
print("Sensitivity sweep on World R (planted truth = -1.60 pp):\n")
print(sens.round(3).to_string(index=False))
print(f"\nAll coefficients within [{sens['coef'].min():.2f}, {sens['coef'].max():.2f}] pp, "
f"all same sign, all |t| > 1.96 -> STABLE.")6. Multiple-testing corrections, done right
Maya does not test one outcome. She has a family: the denial gap, the approval-rate gap, the rate-spread gap, loan amounts, and subgroup breakdowns. If she tests eight outcomes and two come back significant at 5%, how impressed should you be? If all eight were truly null, the chance that at least one clears the bar is 1 - (1-0.05)^8 \approx 0.34 — a one-in-three false “discovery” from pure noise. She owes the reader a correction for the size of her family.
- Bonferroni controls the family-wise error rate (the chance of even one false positive): test each of m hypotheses at \alpha/m. Brutal and simple; conservative to the point of killing power.
- Benjamini–Hochberg (1995) controls the false discovery rate (the expected fraction of your declared discoveries that are false): order the p-values p_{(1)} \le \dots \le p_{(m)}, find the largest rank k with p_{(k)} \le \frac{k}{m}\alpha, and reject all hypotheses with rank \le k. The sliding bar lets a cluster of moderately-small p-values reinforce each other.
We use Maya’s worked family from the chapter (m=8, \alpha=0.05) and confirm BH declares the two smallest significant while the naive-significant rate-spread gap at p=0.030 does not survive.
def benjamini_hochberg(pvals, alpha=0.05):
"""Return a boolean 'reject' mask (FDR control) aligned to the input order."""
p = np.asarray(pvals, float)
m = p.size
order = np.argsort(p) # ascending
ranked = p[order]
bar = (np.arange(1, m + 1) / m) * alpha # k/m * alpha
below = ranked <= bar
if not below.any():
k = 0
else:
k = np.max(np.where(below)[0]) + 1 # largest rank passing the bar
reject = np.zeros(m, bool)
reject[order[:k]] = True
return reject
def bonferroni(pvals, alpha=0.05):
p = np.asarray(pvals, float)
return p <= alpha / p.size
family = pd.DataFrame({
"outcome": ["Denial gap (primary)", "Approval-rate gap", "Rate-spread gap", "Loan-amount gap",
"Subgroup A gap", "Subgroup B gap", "Subgroup C gap", "Subgroup D gap"],
"p": [0.004, 0.011, 0.030, 0.041, 0.20, 0.33, 0.51, 0.78],
})
m, alpha = len(family), 0.05
family = family.sort_values("p").reset_index(drop=True)
family["rank_k"] = np.arange(1, m + 1)
family["BH_bar"] = family["rank_k"] / m * alpha
family["naive_5%"] = family["p"] <= alpha
family["Bonferroni"] = bonferroni(family["p"], alpha)
family["BH (FDR)"] = benjamini_hochberg(family["p"], alpha)
print(f"Maya's family: m = {m} outcomes, alpha = {alpha}\n")
print(family.round(5).to_string(index=False))
print(f"\nNaive 5% declares {family['naive_5%'].sum()} significant.")
print(f"Bonferroni declares {family['Bonferroni'].sum()} significant (bar = alpha/m = {alpha/m:.5f}).")
print(f"BH (FDR) declares {family['BH (FDR)'].sum()} significant.")
print("BH keeps the two smallest; the p=0.030 rate-spread gap does NOT survive correction.")7. Oster (2019) δ: how strong would the unobserved confounder have to be?
The deepest attack confronts the threat you cannot test directly — an omitted variable you never measured. The question is never “is there an unobserved confounder?” (there always could be) but how strong would selection on unobservables have to be, relative to the selection on observables you can see, to explain away the entire result? If the answer is “stronger than the observables you already control for,” the result is robust to a threat you cannot test.
Oster fixes the fatal flaw in naive coefficient stability. A coefficient can stay put when you add controls for two opposite reasons: the controls truly do not matter (reassuring), or the controls matter but are nearly uncorrelated with treatment (proves nothing — an unobservable as predictive could still move the coefficient). Coefficient movement is only meaningful relative to how much the R^2 moved. You read three things off two regressions:
- The uncontrolled regression: coefficient \hat\beta_0 and \tilde R_0.
- The controlled regression: coefficient \hat\beta_1 and \tilde R_1.
- A maximum R^2, R_{\max}, that treatment-plus-all-observables-and-unobservables would reach. You assume it; Oster’s evidence-based default is R_{\max} = 1.3\,\tilde R_1, capped at 1.0.
\delta is the ratio of selection on unobservables to selection on observables. Set the bias-adjusted coefficient \beta^* = 0 and solve for the \hat\delta that would drive the effect to zero. Pass: \hat\delta \ge 1 — an unobservable would have to be at least as confounding as everything you observed, usually implausible. Fail: \hat\delta < 1 — a weaker confounder could erase the effect.
def oster_delta(beta0, R0, beta1, R1, Rmax=None, beta_target=0.0):
"""Oster (2019) delta to move the controlled coefficient to beta_target (default 0).
beta0,R0 = uncontrolled coef & R^2; beta1,R1 = controlled coef & R^2."""
if Rmax is None:
Rmax = min(1.3 * R1, 1.0) # Oster's evidence-based default
num = (beta1 - beta_target) * (R1 - R0)
den = (beta0 - beta1) * (Rmax - R1)
return num / den
def coef_and_r2(df, formula):
res = smf.ols(formula, data=df).fit()
return res.params["post_treat"], res.rsquared
# Uncontrolled = post_treat + FE only; controlled = + observed controls. Both on World R.
UNCONTROLLED = "denial_gap ~ post_treat + C(county) + C(year)"
CONTROLLED = "denial_gap ~ post_treat + income + ltv + C(county) + C(year)"
b0, R0 = coef_and_r2(panel_R, UNCONTROLLED)
b1, R1 = coef_and_r2(panel_R, CONTROLLED)
print(f"World R: uncontrolled coef b0 = {b0:+.3f}, R0 = {R0:.3f}")
print(f" controlled coef b1 = {b1:+.3f}, R1 = {R1:.3f}")
delta_default = oster_delta(b0, R0, b1, R1) # Rmax = 1.3*R1
delta_harsh = oster_delta(b0, R0, b1, R1, Rmax=1.0) # Rmax = 1 (harshest)
print(f"\nOster delta (Rmax = 1.3*R1 = {min(1.3*R1,1.0):.3f}): delta = {delta_default:.2f}")
print(f"Oster delta (Rmax = 1.0, harshest): delta = {delta_harsh:.2f}")
print(f"\n-> {'PASS (delta >= 1: robust to unobserved selection)' if delta_default >= 1 else 'FAIL (delta < 1: fragile)'}")
print("Interpretation: an unobservable would have to be ~%.1fx as confounding as the rich observables"
% delta_default, "to explain the effect away.")\hat\delta is mechanically sensitive to R_{\max}, so the honest report shows the whole sweep, not one number. We trace \hat\delta as R_{\max} rises from the controlled \tilde R_1 (where the test is vacuous, \hat\delta \to \infty) up to 1.0 (the harshest assumption). As long as \hat\delta stays comfortably above 1 across the plausible range, the robustness conclusion holds.
rmax_grid = np.linspace(R1 + 0.01, 1.0, 60)
delta_grid = np.array([oster_delta(b0, R0, b1, R1, Rmax=rm) for rm in rmax_grid])
fig, ax = plt.subplots()
ax.plot(rmax_grid, delta_grid, color="darkgreen", lw=2)
ax.axhline(1.0, color="crimson", ls="--", lw=1.5, label="delta = 1 (Oster heuristic)")
ax.axvline(min(1.3 * R1, 1.0), color="black", ls=":", lw=1, label=f"default Rmax = 1.3*R1")
ax.set_xlabel(r"assumed $R_{max}$")
ax.set_ylabel(r"$\hat\delta$ to drive the effect to zero")
ax.set_title("Oster delta sensitivity to Rmax (World R): delta stays well above 1")
ax.set_ylim(0, min(np.nanmax(delta_grid) * 1.1, delta_grid[len(delta_grid)//2] * 4))
ax.legend()
fig.tight_layout()
fig.savefig("nb82_oster_sweep.png", dpi=110)
plt.close(fig)
print("Saved figure: nb82_oster_sweep.png")
print(f"delta at harshest Rmax=1.0: {delta_grid[-1]:.2f} ({'>=1, still robust' if delta_grid[-1] >= 1 else '<1, fragile here'})")8. Assembling the battery: what a survivor looks like
Each test answered a row of Maya’s Chapter 7.5 threats table. The summary below pulls the verdicts together for World R (the real effect) and contrasts the one test that fails in World S (the spurious effect). A result that survives the whole battery is not proven — nothing observational ever is — but it is credible in the precise sense the identification memo demanded: you tried every attack a referee would try, on yourself first, and the number held.
And the converse is the real lesson: a result that does not survive should be reported as not surviving. World S has a clean, significant, stable coefficient — and the in-time placebo lights up, telling you the whole thing is a pre-trend. A failed robustness test is a finding, not a setback.
battery = pd.DataFrame([
("Differential pre-trend", "In-time placebo",
f"PASS (t={tR_pl:+.2f})", f"FAIL (t={tS_pl:+.2f})"),
("Few treated clusters", "Wild cluster bootstrap",
f"PASS (p={p_wild:.3f})", "-- (see few-cluster world)"),
("Errors correlated across units", "Alternative clustering",
f"PASS (state SE {state_se:.2f} > classical {classical_se:.2f}, CI excl. 0)", "n/a"),
("Effect is random noise", "In-space placebo",
f"PASS (perm p={p_perm_R:.3f})", "n/a"),
("Result tuned to a knob", "Sensitivity sweep",
f"PASS (coef in [{sens['coef'].min():.2f},{sens['coef'].max():.2f}])", "n/a"),
("Cherry-picked from a family", "Benjamini-Hochberg",
f"{family['BH (FDR)'].sum()} of {m} survive FDR", "n/a"),
("Unobserved confounder", "Oster (2019) delta",
f"PASS (delta={delta_default:.2f} >= 1)", "n/a"),
], columns=["Threat (Ch 7.5 table)", "Stress-test", "World R (real)", "World S (spurious)"])
print("THE ROBUSTNESS BATTERY -- verdicts\n")
print(battery.to_string(index=False))
print("\nWorld R survives every attack. World S has an identical-looking headline coef but its")
print("in-time placebo FAILS -> the 'effect' is a differential pre-trend, and the honest report says so.")Your Turn
You now have a working battery. Adapt it to your pre-registered estimate from nb8.1 and run column 3 of your Chapter 7.5 threats table, end to end.
Run your scariest placebo first. Take the threat you ranked most dangerous and build the placebo that would expose it — the in-time fake date if you fear a pre-trend, the in-space permutation if you fear few clusters. Edit
make_panel’spretrend_gapto inject a confounder of the size you actually fear, and check whether your in-time placebo catches it. If it fails, write the honest sentence that narrows or kills your claim.Read your coefficient stability with Oster’s eyes. Compute \hat\delta for your result, then re-run the R_{\max} sweep. Did your coefficient stay stable because the controls were powerful (large \hat\delta, reassuring) or toothless (small \hat\delta, proves nothing)? Defend your R_{\max}, then report \hat\delta at R_{\max} = 1.
Audit your family. Write down every outcome you looked at — not just the ones you plan to report. Apply both
bonferroniandbenjamini_hochberg. Which survive each? If your headline survives only the uncorrected single test, what is the honest contribution sentence — and is it the one you were hoping to write?
The author who runs these tests hoping they pass and the author who runs them hoping to learn the truth get the same code and the same output. Only the second one is doing science.