Robustness workflow¶
A publication-quality causal inference result is not one estimate. It
is a package: the point estimate + a stack of robustness checks that
convince a reviewer the conclusion doesn't hinge on a single modelling
choice. StatsPAI's .next_steps() gives you a priority-ordered list.
This guide expands the full recipe for each design.
The three-layer defence¶
Layer 1 — IDENTIFICATION: is my design valid?
DID: parallel trends
RD: continuity, no manipulation
IV: relevance, exclusion, monotonicity
SOO: unconfoundedness, overlap
Layer 2 — SPECIFICATION: does the result survive alternative choices?
Bandwidth, polynomial order, covariate set, kernel, sample restriction
Layer 3 — SENSITIVITY: how much unobserved confounding would
kill the conclusion?
E-value, Oster, sensemakr, Rambachan-Roth
Layer 1 — Identification checks¶
DID¶
# Placebo pre-trends (per-period joint)
sp.pretrends_test(r)
sp.bjs_pretrend_joint(r) # Borusyak-Jaravel-Spiess
# Honest inference under parallel-trends violations
sp.honest_did(r, m_grid=[0.0, 0.25, 0.5], method='smoothness')
# Decomposition: which 2x2 cells drive TWFE?
sp.bacon_decomposition(df, y='y', treat='d', time='t', id='i')
RD¶
# Density continuity (McCrary 2008 / Cattaneo-Jansson-Ma 2020)
sp.rddensity(df, x='running_var', c=0.0)
sp.mccrary_test(df, x='running_var', c=0.0)
# Covariate balance
sp.rdbalance(df, x='running_var', c=0.0, covs=['X1', 'X2'])
# Placebo cutoffs
sp.rdplacebo(df, y='y', x='running_var',
c=0.0, placebo_cutoffs=[-0.5, 0.5])
IV¶
# First-stage strength
r.diagnostics['First-stage F (d)'] # > 10 rule; > 30 comfortable
# Exclusion sensitivity
sp.iv(method='ltz', y='y', endog='d', instruments='z',
gamma_mean=0.0, gamma_var=0.01, data=df)
# Monotonicity (Imbens-Rubin / Kitagawa)
sp.kitagawa_test(df, y='y', treatment='d', instrument='z')
# Reduced form
sp.regress('y ~ z', data=df).tidy()
Selection-on-observables¶
# Overlap / common support
sp.overlap_plot(r)
sp.love_plot(r)
sp.ps_balance(r)
# E-value: how strong must unobserved U be to explain away effect?
sp.evalue_from_result(r)
# Oster delta bounds
sp.oster_bounds(r, delta=1.0)
# Sensemakr (Cinelli-Hazlett)
sp.sensemakr(df, y='y', treat='d', controls=['X1', 'X2'], benchmark=['X1'])
Layer 2 — Specification robustness¶
Specification curve / multiverse¶
spec = sp.spec_curve(
df, y='y', x='d',
controls=[
['X1'], ['X1', 'X2'], ['X1', 'X2', 'X3'],
['X1', 'X2', 'X3', 'X4'],
],
se_types=['nonrobust', 'HC1'],
)
spec.plot() # classic Simonsohn-Simmons-Nelson spec curve
Bandwidth / kernel sensitivity (RD)¶
Subsample stability¶
sp.robustness_report(df, formula='y ~ d + x1 + x2', x='d',
subsets={'male': 'sex == 1'},
winsor_levels=[0.01], trim_pct=0.01)
Layer 3 — Sensitivity to unobserved confounding¶
E-value (VanderWeele-Ding 2017)¶
For any causal estimate on observational data, the E-value is the minimum strength of association an unobserved confounder U would need to have with both treatment and outcome (on the risk-ratio scale) to fully explain away the observed effect.
# From a fitted causal result (interpreted as a standardised mean diff):
r = sp.dml(df, y='y', treat='d', covariates=[...])
ev = sp.evalue_from_result(r)
print(ev["evalue_estimate"]) # e.g. 2.3 -> U needs RR > 2.3 with both
print(ev["evalue_ci"]) # E-value for the CI limit nearest the null
# Or directly from a reported effect on any supported scale:
sp.evalue(estimate=1.8, ci=(1.2, 2.7), measure="OR") # rare=False by default
sp.evalue(estimate=1.5, ci=(1.1, 2.0), measure="HR", rare=False)
sp.evalue(estimate=0.5, se=0.1, sd=2.0, measure="OLS")
sp.evalue(estimate=2.5, ci=(1.8, 3.2), measure="RR", true=1.5) # non-null E-value
sp.evalue_rd(200, 150, 100, 250) # exact risk-difference E-value from a 2x2 table
sp.evalue returns a dict (evalue_estimate, evalue_ci, rr_estimate,
interpretation, ...) and reproduces the R EValue package to machine
precision across RR/OR/HR (rare and common outcomes), MD/SMD, OLS, and the
2x2-table risk difference. If the confidence interval already contains the
null (or the specified true value), its E-value is exactly 1.
Rule of thumb: E > 2 is "reasonably robust"; E > 3 is strong; E < 1.5 is fragile (any of the measured covariates already has that level of association).
Oster (2019) coefficient stability bounds¶
b = sp.oster_bounds(r, delta=1.0, r_max=1.3 * r.glance()['r_squared'].iloc[0])
print(b.lower, b.upper)
[lower, upper] is where the true effect must lie if
unobservables are as important as observables (δ=1). If 0 is in this
interval, the effect is not robust.
Sensemakr (Cinelli-Hazlett 2020)¶
s = sp.sensemakr(df, y='y', treat='d', controls=['X1', 'X2', 'X3'],
benchmark=['X1', 'X2'])
s.summary()
s.plot() # contour plot: partial R^2 of U with T and Y
Bounds (Manski-style)¶
When no assumption feels credible, compute worst-case bounds:
Full robustness report in one call¶
report = sp.robustness_report(
df, formula='y ~ d + x1 + x2', x='d',
cluster_var='state',
extra_controls=['x3'], drop_controls=['x2'],
winsor_levels=[0.01, 0.05],
)
report.summary() # every variant, side by side
report.results_df # tidy DataFrame, one row per specification
report.plot(); report.to_latex()
report.baseline_estimate / .baseline_se pin the reference specification and .n_checks counts the variants run.
This covers Layer 2 for a regression-style design. The Layer 1 and
Layer 3 checks above are separate calls — there is no single function
that returns a verdict.
Reporting template¶
For the final paper, present:
- Main table: point estimate, SE, 95% CI, N.
- Identification check: one test per identifying assumption (pre-trend plot, McCrary plot, first-stage F).
- Specification panel: point estimate varying key design choices (covariate sets, bandwidths, subsamples).
- Sensitivity panel: E-value, Oster bounds, sensemakr contour.
sp.regtable() auto-formats this for LaTeX/Word: