Use when the user asks to run a full empirical / causal analysis in Python — by default in the style of an applied economics paper (AER / QJE / JPE / ReStud / AEJ) with DID / RD / IV / SCM / DML / matching, written-out estimating equation + identifying assumption, Table 1 / Table 2 / event-study figure / robustness gauntlet — OR in epidemiology / public health style (target-trial emulation, IPTW + g-formula + TMLE triplet, Mendelian randomization, KM/AFT survival, E-value sensitivity, STROBE/TRIPOD reporting) — OR in ML causal inference style (DML, S/T/X/R/DR meta-learners, causal forest, Dragonnet/TARNet/CEVAE, BCF, CATE distribution, policy learning, conformal causal, fairness audit, causal discovery) — OR in distributional / gap-decomposition style (Oaxaca–Blinder `sp.oaxaca`, Kitagawa `sp.kitagawa_decompose`, DiNardo–Fortin–Lemieux `sp.dfl_decompose`, Gelbach `sp.gelbach`, Fairlie `sp.fairlie`, RIF / FFL `sp.rif_decomposition`, all reachable through the `sp.decompose` dispatcher). Also covers exporting mu
Install with Codex or Claude Copy this prompt, paste it into Codex, Claude, or another assistant, and let it review the skill page and install it for you.
A direct command skips the review prompt. Inspect the source before running it.
StatsPAI is a validation-tiered Python package for causal inference and applied econometrics: one import statspai as sp, 1,100+ registered functions behind a self-describing API, and mature estimator result objects that commonly export to LaTeX / Word / Excel / BibTeX.
This skill drives StatsPAI through the canonical pipeline of an applied AER empirical paper. Each step emits a paper-ready artifact (Table 1, event-study figure, Table 2 main results, robustness panel, replication stamp).
Install: pip install "statspai[fixest,plotting]" (API surface re-validated against statspai 1.19.0 — every sp.* reference, signature, and result-object attribute claim in this skill is checked by validate_api_claims.py in this folder). The bare pip install statspai is not enough for the default pipeline — see the dependency matrix below.
Paper: JOSS submission under review; JSS materials in Paper-JSS/README.md and docs/jss_source_audit_dossier.md
Install the right extras or the documented calls will raise ImportError. Several core functions live behind optional dependency groups (verified from pyproject.toml):
You use…
Needs extra
Install
Symptom if missing
sp.feols / sp.fepois / sp.feglm (high-dim FE — the default for any y ~ x | fe regression)
fixest (pyfixest)
pip install "statspai[fixest]"
ImportError: pyfixest is required …
Any figure (sp.coefplot, sp.binscatter, event-study/RD/SCM plots, .plot())
plotting (matplotlib/seaborn)
pip install "statspai[plotting]"
ImportError on first plot
This SKILL.md is very large, so SkillsMP previews the first section here.View on GitHub
A one-shot install covering the whole skill: pip install "statspai[fixest,plotting,neural,text]". sp.regtable / sp.collect / Word+Excel+LaTeX export, sp.regress, IV, RD, DID (callaway_santanna), matching, DML, meta-learners, causal forest, BCF, TMLE, and the epi stack work on the base install.
Verified skeleton (copy, then swap in your columns)
This minimal pipeline runs start-to-finish against statspai 1.19.0 (every call below was executed). It is the golden path — adapt column names / design, keep the call shapes and the unpack-then-save figure idiom. The full playbook (§−1 → §8) expands each step.
import numpy as np, pandas as pd, statspai as sp
# df has: wage, training(0/1), worker_id, firm_id, year, first_treat_year, age, edu, tenure, ...# §1 Table 1 → Word/Excel/LaTeX
mc = sp.mean_comparison(df, ["age","edu","tenure"], group="training", test="ttest",
title="Table 1. Summary statistics")
mc.to_word("tables/table1.docx"); mc.to_excel("tables/table1.xlsx")
# §2 Estimand-first plan (freeze BEFORE estimating)
q = sp.causal_question(treatment="training", outcome="wage", data=df, estimand="ATT",
design="did", time_structure="panel", time="year", id="worker_id",
covariates=["age","edu","tenure"])
plan = q.identify(); print(plan.summary())
# §3 Identification figure — from a CS/SA result (NOT event_study()); plotters return (fig, ax)
cs = sp.callaway_santanna(df, y="wage", g="first_treat_year", t="year", i="worker_id", x=["age","edu"])
fig, ax = sp.enhanced_event_study_plot(cs, shade_pre=True); fig.savefig("figures/fig2a.png", dpi=300)
# §4 Main table — mix sp.regress (no FE) + sp.feols (HDFE, needs statspai[fixest]) in ONE regtable
M1 = sp.regress("wage ~ training", df, cluster="firm_id")
M2 = sp.feols("wage ~ training + age + edu + tenure | industry + year", df, vcov={"CRV1":"firm_id"})
rt = sp.regtable(M1, M2, template="aer", coef_labels={"training":"Job training"},
model_labels=["(1) OLS","(2) FE"], stats=["N","R2","Cluster","FE"],
title="Table 2. Effect of training on wages")
rt.to_word("tables/table2.docx"); rt.to_excel("tables/table2.xlsx")
open("tables/table2.tex","w").write(rt.to_latex())
# §5 Heterogeneity — per-row CATE at result.model_info["cate"] (there is NO .cate_estimates)
ml = sp.metalearner(df, y="wage", treat="training", covariates=["age","edu","tenure"], learner="dr")
fig, ax = sp.cate_plot(ml, kind="hist"); fig.savefig("figures/fig4.png", dpi=300)
# §7 Robustness — Oster + E-value + honest-DID sensitivity figure
sp.oster_bounds(data=df, y="wage", treat="training", controls=["age","edu","tenure"], r_max=1.3)
sp.evalue(estimate=M2.params["training"], ci=tuple(M2.conf_int().loc["training"]), measure="RR")
fig, ax = sp.sensitivity_plot(sp.honest_did(cs, method="smoothness"),
original_estimate=cs.estimate, original_ci=cs.ci)
fig.savefig("figures/fig6.png", dpi=300)
# §8 One-file replication bundle (Word/Excel/LaTeX/Markdown from one source)
c = sp.collect("Replication", template="aer")
c.add_summary(df, vars=["wage","age","edu","tenure"], stats=["mean","sd","n"], title="Table 1")
c.add_regression(M1, M2, model_labels=["(1)","(2)"], stats=["N","R2"], title="Table 2")
for ext in ("docx","xlsx","tex","md"): c.save(f"replication/paper.{ext}")
Epi (§A) and ML-causal (§B) reuse this exact scaffolding — only the §4 estimator stack changes (TMLE/g-formula/MR for epi; DML/meta-learner/causal-forest for ML), and every estimator still returns a result that drops into sp.regtable / sp.collect.
Why for Agents
Self-describing: sp.list_functions() / sp.describe_function(name) / sp.function_schema(name) — registered symbols are discoverable without doc lookup.
Structured results: mature estimators return result objects with methods such as .summary(), .plot(), .diagnostics, .to_latex(), .to_word(), .cite() when supported.
One import, full pipeline: data contract → Table 1 → estimand-first DSL → identification graphs → main table → heterogeneity → mechanisms → robustness → replication package.
Estimand-first: sp.causal_question(...).identify() forces the "DID vs RD vs IV?" decision before estimation, with the identifying assumption written down — the way a referee expects to read it.
SkillOpt-derived operating loop (read before the playbook)
SkillOpt's useful lesson for this skill is procedural, not cosmetic: a skill is a
bounded decision policy that should improve from rollout evidence while preserving
verified behavior. Treat every StatsPAI request as a mini rollout:
Route the mode first: choose Default/AER, Mode A/epi, Mode B/ML-causal, or
a narrow export-only path from the user's words. Do not run the full paper
pipeline when the request is only "make Table 1" or "export this regression".
Freeze the contract before estimating: name y, treatment/exposure,
unit/time ids, estimand, design, required artifacts, and install extras. If any
field is missing, infer only when the column names make the choice obvious;
otherwise produce a short blocking checklist instead of hallucinating columns.
Start from the smallest verified call shape: prefer the skeleton and the
relevant section-specific snippet over ad hoc API guesses. For an unfamiliar
function, call sp.describe_function(name) / sp.function_schema(name) before
writing code.
Widen one block at a time: data contract → plan → diagnostic figure →
main estimate → robustness/export. After each block, read warnings and object
attributes before passing the result downstream.
Gate the answer on artifacts, not intentions: final responses should list
the files produced, the identifying assumption, the estimator class, and any
failed or skipped gate. Never claim "paper-ready" if Word/Excel/LaTeX exports
or required diagnostics were not actually generated.
Turn failures into bounded corrections: if a call raises, fix the smallest
wrong rule (signature, result type, optional extra, plot return shape) and
continue from the last verified artifact. Do not rewrite the pipeline wholesale.
SkillOpt-style execution gate (task-local card)
Before generating or revising StatsPAI analysis code, compress the request into a task-local best_skill card:
best_skill: <mode + design + artifact target>
train_signal: <current failure, user goal, or missing evidence>
selection_split: <focal dataset/spec/output used to judge the candidate>
heldout_gate: <checks the patch must pass beyond the focal example>
accepted_patterns: <rules to reuse after validation>
rejected_patterns: <failed shortcuts not to retry without new evidence>
patch_scope: <one estimator/sample/export/robustness change>
reject_if: <conditions that force rollback to the last passing spec>
Route card: record the mode (econ, epi, or ml-causal), estimand, identification design, focal outcome/treatment, StatsPAI install extras, and required artifacts.
Bounded edit: change one decision at a time (sample rule, estimator, optional extra, plot return shape, export format, or robustness check). Prefer the smallest patch that can pass validation.
Selection split discipline: treat the user's immediate failure or requested artifact as the selection split. Reserve at least one alternate outcome, sample window, estimator family, or export target as the held-out gate.
Held-out gate: define checks before running code: row counts, key uniqueness, treatment support, missingness thresholds, expected table/figure files, and one non-focal robustness/specification that the change must not break.
Reject buffer: if a candidate spec fails the gate, log the failure, code diff, and gate output in analysis_log.md; revert to the last passing spec and do not retry the same unchecked pattern.
Slow/meta update: at the end of the task, write down accepted_patterns and rejected_patterns from the trajectory. Do not widen the canonical project template from a single passing run.
Promote only after validation: only turn a one-off fix into reusable project boilerplate after it passes the current data and at least one alternate outcome/sample/specification.
Acceptance gates by request type
Request type
Minimum gates before final answer
Export-only / outreg2 equivalent
At least one RegtableResult or Collection object is created; requested .docx / .xlsx / .tex paths are written or the exact missing optional dependency is reported
AER DID / event study
sp.causal_question(...).identify() saved or printed; CS/SA result used for the event-study figure; numerical pre-trends checked separately with sp.event_study(...) or equivalent; Table 2 and at least one robustness/sensitivity artifact produced
IV
First-stage F and instrument story reported before the 2SLS coefficient; no | fe formula is passed to sp.ivreg; FE-IV needs explicit dummy construction or a stated limitation
RD
McCrary/manipulation check plus RD plot are produced before the treatment-effect table; bandwidth/kernel sensitivity is in the robustness block
Matching / weighting
Balance or love plot is produced before outcome estimation; weights are carried into the Table 1 / balance export when applicable
Epi / target-trial
Target-trial protocol is written before modeling; positivity/overlap is checked; IPTW/g-formula/TMLE estimates are compared when data support them; E-value or equivalent sensitivity is reported
ML causal / CATE
Train/holdout split and nuisance learners are explicit; per-row CATE source is valid (model_info["cate"] for meta-learners or cf.effect(X) for forests); policy/OPE claims use holdout data
Stata/R migration
Use StatsPAI's self-description or translator surface first; preserve semantic notes for unsupported options instead of silently pretending full parity
Maintenance rule for future skill edits
When improving this skill itself, follow a SkillOpt-style accept rule: propose a
small add/delete/replace edit, then accept it only if it helps a concrete failure
case and does not regress the verified skeleton, export cookbook, or Common
Mistakes table. Use EVALS.md as the held-out gate set for future skill edits.
Keep reusable fixes near the earliest section where an agent will need them; keep
rare API traps in Common Mistakes.
The AER-style empirical pipeline
The skill mirrors the canonical sections of an applied AER / QJE / AEJ paper. Each step below is one paper section and one set of artifacts on disk.
All code blocks below share one running example (training → wage, with worker_id / firm_id / year / age / edu / tenure) purely for readability. Column names, population, estimand, and design values are illustrative — substitute the user's actual columns and research question. Only sp.* function names and argument shapes are normative.
The default playbook above is AER-style applied econometrics — the AEA convention: written-out estimating equation, identifying assumption table, design horse-race, full robustness gauntlet. The skill also ships two parallel sub-pipelines for the other two big causal-inference traditions, each reusing the same export stack (sp.regtable / sp.collect / sp.paper_tables) and result objects:
"Mix" (e.g. "estimate DID + then ML CATE on the heterogeneity")
Default + Mode B in sequence — every estimator returns the same CausalResult, drop them all into one sp.regtable(...) for the horse-race column
The three modes share the same export stack, the same CausalResult interface, and the same sp.causal_question(...).identify() estimand-first DSL — switching modes only changes which Step 4 estimators you reach for, not the surrounding scaffolding. If you only want descriptive stats / Table 1 / a balance check, the AER sp.sumstats / sp.mean_comparison / sp.collect calls work in all three modes.
Paper-ready figure & table inventory (what to produce by section)
A modern AER paper has 5–7 figures and 3–5 main tables + an appendix robustness table. Every step below should leave at least one numbered artifact on disk. Default file names assume parallel .tex / .docx / .xlsx exports (the agent should produce all three so co-authors can edit in Word / Excel and the build system can use LaTeX):
sp.regtable(rob1...robN, panel_labels=[...]).to_word/.to_excel — or sp.paper_tables(robustness=[...]).to_docx()
tables/tableA1_robustness.{tex,docx,xlsx}
§7
Figure 5: spec curve
sp.spec_curve(...).plot()
figures/fig5_spec_curve.png
§7
Figure 6: honest-DID sensitivity plot (+ text dashboard)
sp.sensitivity_plot(sp.honest_did(cs, ...)) for the figure; print(sp.sensitivity_dashboard(result).summary()) for the Cinelli–Hazlett/Oster/E-value numbers (text, not a figure)
figures/fig6_sensitivity.png
§8
Replication bundle: all tables in one Word/Excel/LaTeX file
sp.collect("Paper").add_summary(...).add_regression(...)...save("paper.{docx,xlsx,tex}") — or sp.paper_tables(main=, heterogeneity=, robustness=, placebo=).to_docx/.to_xlsx
replication/paper.{docx,xlsx,tex}
Every CausalResult and OLS model can be passed straight into sp.regtable(...), sp.coefplot(...), and sp.collect(). Don't hand-roll LaTeX, and don't render Word/Excel from pandas — the export functions apply book-tab borders, AER-style stars, and the right SE label automatically.
Export cookbook — Word / Excel / LaTeX in one line
StatsPAI's export stack is the agent-native equivalent of Stata's outreg2 / esttab / collect and R's modelsummary / gtsummary. Three tiers, picked by scope of what you're exporting:
Tier
Use when
API
Hot kwargs
1. Single multi-column table (the outreg2 / summary_col equivalent)
Exporting one Table 2 / Table 3 / Table A1 with progressive columns
Naming gotcha: sp.regtable(..., output="docx") is invalid — the enum is {"text", "latex", "tex", "html", "markdown", "md", "qmd", "quarto", "word", "excel"}. Use output="word" / "excel", or — simpler — drop output= and call .to_word(filename) / .to_excel(filename) on the result.
Notebook setup — CJK fonts + retina DPI
Run once at the top of every analysis script / notebook, before any matplotlib-backed plot (sp.regtable.to_* exporters do not need this — only .savefig / sp.coefplot / sp.binscatter / sp.cate_plot / etc.). Two failures it fixes in one shot:
CJK labels render as ▢▢▢ tofu — the matplotlib default DejaVu Sans carries no Chinese / Japanese / Korean glyphs, so ax.set_title("教育回报") silently degrades into squares.
Plots look fuzzy on hi-DPI displays — matplotlib's default figure.dpi=100 is half the density of a Retina / 4K screen.
Drop-in snippet
import matplotlib as mpl
import matplotlib.pyplot as plt
defsetup_plot(retina: bool = True) -> None:
"""One-shot matplotlib boilerplate: CJK font fallback + retina DPI.
Idempotent — safe to call multiple times. Call BEFORE any plotting.
"""# 1. CJK font fallback chain — covers macOS / Windows / Linux in one list.# matplotlib uses the first available font; later names are fallbacks,# so listing all three platforms is harmless on any single host.
mpl.rcParams["font.sans-serif"] = [
"PingFang SC", "Heiti SC", "Hiragino Sans GB", # macOS"Microsoft YaHei", "SimHei", "SimSun", # Windows"Noto Sans CJK SC", "Source Han Sans SC", # Linux / Adobe"WenQuanYi Micro Hei", # Linux fallback"Arial Unicode MS", # universal fallback"DejaVu Sans", # last-resort Latin
]
mpl.rcParams["axes.unicode_minus"] = False# 修复中文字体下负号渲染成 □# 2. Retina-grade DPI. figure.dpi controls on-screen / inline rendering;# savefig.dpi controls .png exports. Set both — they are independent.if retina:
mpl.rcParams["figure.dpi"] = 144# 2× default — sharp on Retina/HiDPI
mpl.rcParams["savefig.dpi"] = 300# manuscript/export PNG (AER house norm)# Jupyter inline retina backend (no-op outside IPython):try:
from IPython import get_ipython
ipy = get_ipython()
if ipy isnotNone:
ipy.run_line_magic("config", "InlineBackend.figure_format = 'retina'")
except Exception:
pass
setup_plot() # call once at the top
Smoke test (5 seconds, run once after setup_plot())
If the saved PNG shows Chinese characters cleanly and the y-axis tick -1 is a real minus sign (not a square), the setup is good. Otherwise see troubleshooting below.
Saving figures — the (fig, ax) idiom (READ THIS)
Every StatsPAI plotter and every result .plot() returns a (fig, ax) tuple — NOT a bare Figure. So sp.parallel_trends_plot(...).savefig(...) raises AttributeError: 'tuple' object has no attribute 'savefig'. Always unpack, then save the figure:
sp.kaplan_meier(...).plot() returns a bare Axes (it is a KMResult, not a CausalResult) — save via ax = km.plot(); ax.figure.savefig(...).
This applies uniformly to coefplot, binscatter, rdplot, rddensity().plot(), bacon_plot, enhanced_event_study_plot, did_summary_plot/ggdid/group_time_plot, synthdid_plot, cate_plot, cate_group_plot, dose_response().plot(), sensitivity_plot, match().plot(), synth().plot(), and a generic result.plot(). The code blocks below all use the unpack-then-save form.
Troubleshooting
Symptom
Fix
Title still shows ▢▢▢ tofu after setup_plot()
Host has none of the listed fonts. Install one — macOS: pre-installed (no action). Linux: sudo apt install fonts-noto-cjk (Debian/Ubuntu) or sudo dnf install google-noto-sans-cjk-fonts (Fedora/RHEL). Windows: pre-installed. Then clear matplotlib's font cache: rm -rf ~/.cache/matplotlib (Linux/macOS) / %LOCALAPPDATA%\matplotlib (Windows), and restart the Python / Jupyter kernel.
Negative numbers render as ▢
axes.unicode_minus = False was overridden by a later plt.style.use(...) or mpl.rcParams.update(...). Re-call setup_plot() after any style change.
Plot blurry inside VSCode .ipynb
VSCode's notebook UI ignores figure.dpi for inline rendering. Either switch the cell output to "Open in Image Viewer", or use %matplotlib inlinebeforesetup_plot(). The saved .png (driven by savefig.dpi=300) is sharp regardless.
sp.<plot>(...) output still shows tofu
The sp.* plotters honor global rcParams, so this only happens when setup_plot() was called after the plot was drawn. Move the call to the very top of the script.
Need to verify which font matplotlib picked
mpl.font_manager.findfont(mpl.font_manager.FontProperties(family=mpl.rcParams["font.sans-serif"])) returns the resolved file path — if it ends in DejaVuSans.ttf despite Chinese labels, no CJK font is installed.
Persist as project default (optional)
Drop the same rcParams into a project-level matplotlibrc next to pyproject.toml so co-authors and CI runners pick it up without calling setup_plot():
# matplotlibrc — committed to the repo
font.sans-serif: PingFang SC, Heiti SC, Microsoft YaHei, SimHei, Noto Sans CJK SC, Arial Unicode MS, DejaVu Sans
axes.unicode_minus: False
figure.dpi: 144
savefig.dpi: 300
The setup_plot() function above is the in-script fallback when a project matplotlibrc is not present.
Step −1 — Pre-Analysis Plan (pre-data; AEA RCT Registry style)
sp.power(design, n=..., effect_size=..., power_target=...) is a unified dispatcher — leave one argument None to solve for it (sample size, MDE, or power). Convenience wrappers: sp.power_rct, sp.power_did, sp.power_rd, sp.power_iv, sp.power_cluster_rct, sp.power_ols.
# Always go through the dispatcher when you want auto-solve. The# `sp.power_<design>` wrappers (power_rct / power_did / power_rd /# power_iv / power_cluster_rct / power_ols) accept *only* the design's# native arguments — they will NOT solve for power_target / n / effect# unless you go via `sp.power(design, ..., power_target=...)`.
sp.power("rct", effect_size=0.3, power_target=0.80) # → PowerResult(n=349, power=0.80)
sp.power("did", n=200, effect_size=0.15, power_target=0.80,
n_periods=4, n_treated_periods=2) # DID: solves MDE / n / power
sp.power("cluster_rct", cluster_size=50, icc=0.05,
effect_size=0.2, power_target=0.80) # Cluster RCT: solves n_clusters# Roth (2022) pre-trends power is a POST-estimation diagnostic — it needs an estimated# event-study result, so run it in §3 once you have `es = sp.event_study(...)`:# sp.pretrends_power(es)
Persist the PowerResult next to data_contract.json and empirical_strategy.md — a referee will ask whether the design was powered before data collection, not after.
Step 0 — Sample construction & data contract (Section "Data")
An AER §1 Data section has three jobs: (a) describe sources, (b) document every sample restriction (the "footnote 4" sample log), (c) lock the panel structure. StatsPAI assumes an analysis-ready DataFrame — do ETL (imputation, type coercion, merges, transforms) in pandas first, then run the 5-check contract.
Paste this log verbatim as footnote 4 of your paper. AER reviewers use it to reconstruct the analysis sample.
0.2 Five-check data contract (go / no-go gate)
import pandas as pd, numpy as np, statspai as sp
defdata_contract(df, *, y, treatment, id=None, time=None, covariates=()):
"""Return a go/no-go dict. Stop the pipeline if any required check fails."""
keys = [y, treatment] + ([id, time] ifidand time else []) + list(covariates)
c = {
"n_obs": len(df), # 1. shape"dtypes": df[keys].dtypes.astype(str).to_dict(), # 2. dtypes on keys"n_missing": df[keys].isna().sum().to_dict(), # 3. missing pattern"n_dupes_on_keys": 0,
"panel_balanced": None,
"cohort_sizes": None,
}
ifidand time:
c["n_dupes_on_keys"] = int(df.duplicated([id, time]).sum()) # 4. duplicate (id,time)
balanced = sp.balance_panel(df, entity=id, time=time) # 5. panel balance
c["panel_balanced"] = len(balanced) == len(df)
c["n_dropped_by_balance"] = len(df) - len(balanced)
if"first_treat_year"in df.columns: # staggered cohorts
c["cohort_sizes"] = (
df.drop_duplicates(id).groupby("first_treat_year").size().to_dict()
)
c["y_range"] = (float(df[y].min()), float(df[y].max()))
c["treatment_share"] = float(df[treatment].mean())
# Missingness mechanism hint (Rubin): compare covariate means between# rows missing-on-y vs observed. Any p < 0.05 ⇒ NOT MCAR → use MI / IPW,# not listwise deletion.from scipy import stats
miss_y = df[y].isna()
c["mcar_hint"] = "likely MCAR (listwise OK)"if miss_y.any() and (~miss_y).any():
for cov in covariates:
if df[cov].dtype.kind in"fi":
_, p = stats.ttest_ind(df.loc[miss_y, cov].dropna(),
df.loc[~miss_y, cov].dropna(),
equal_var=False)
if p < 0.05:
c["mcar_hint"] = f"NOT MCAR (y-miss differs on {cov}, p={p:.3f}) → use MI / IPW"breakreturn c
contract = data_contract(df, y="wage", treatment="training",
id="worker_id", time="year",
covariates=["age", "edu", "tenure"])
assert contract["n_dupes_on_keys"] == 0, "duplicate (id, time) — fix before panel methods"assertall(v == 0for v in contract["n_missing"].values()), \
f"NaNs on keys: {contract['n_missing']}"
If any assertion fires, stop and fix it in pandas — StatsPAI estimators silently drop NaN rows, the most common source of "mysterious sample-size shrinkage" bugs. Persist:
The signature AER Table 1 has three column blocks plus a difference column:
| | (1) Full | (2) Treated | (3) Control | (4) Δ (t-test) |
The Imbens–Rubin rule of thumb: a normalized difference |Δ| / √((s²₁+s²₀)/2) > 0.25 flags substantive imbalance and should trigger matching / reweighting before you trust an OLS comparison.
# Quick text/LaTeX preview (use sumstats `output=` for a string-only render).# When `by=` is binary 0/1 and you don't pass `by_labels=`, sumstats auto-fills# the panel headers as **Control / Treated** so the academic Table 1 reads# correctly out of the box. For non-0/1 codings or different wording, pass# `by_labels={0:"Untrained", 1:"Trained"}` (or `{"A":"Control","B":"Treated"}`).print(sp.sumstats(df, vars=["wage","edu","exp","tenure","age"],
by="training", output="text"))
# AER-style balance table → Word + Excel + LaTeX in three lines.# `mean_comparison` returns a MeanComparisonResult that exposes the full# export chain (.to_word / .to_excel / .to_latex / .to_markdown / .to_html).
mc = sp.mean_comparison(df,
["age","edu","tenure","firm_size"],
group="training",
test="ttest",
title="Table 1. Summary statistics by treatment status")
mc.to_word ("tables/table1_summary.docx") # editable in Word
mc.to_excel("tables/table1_summary.xlsx") # editable in Excelopen("tables/table1_summary.tex", "w").write(mc.to_latex())
sp.describe(df).to_markdown("references/codebook.md") # auto-codebook
1.1 Multi-panel Table 1 (AER convention)
Group rows into Panel A: Outcomes, Panel B: Treatment intensity, Panel C: Controls, Panel D: Sample composition. The cleanest path is to push each panel into a sp.collect() bundle — one .save("file.docx") call then writes the whole multi-panel Table 1 with AER book-tab borders, in Word and Excel and LaTeX from one source.
panels = {
"A. Outcomes": ["wage", "log_wage", "weeks_employed"],
"B. Treatment": ["training", "training_hours"],
"C. Demographic controls": ["age", "edu", "female", "married"],
"D. Labor market": ["tenure", "firm_size", "industry_id"],
}
c1 = sp.collect("Table 1. Summary statistics", template="aer")
for label, vs in panels.items():
c1.add_heading(f"Panel {label}", level=2)
c1.add_summary(df, vars=vs, stats=["mean", "sd", "n"])
c1.save("tables/table1_summary.docx") # editable Word, AER book-tab borders
c1.save("tables/table1_summary.xlsx") # one sheet per panel (heading drives the sheet name)
c1.save("tables/table1_summary.tex") # multi-panel LaTeX# Plain-text alternative (no Collection): one `sp.sumstats` per panel, concat strings.# Useful when you only need the .tex preview without a binary export.import io; buf = io.StringIO()
for label, vs in panels.items():
buf.write(f"\n% Panel {label}\n")
buf.write(sp.sumstats(df, vars=vs, by="training",
stats=["mean", "sd", "n"], output="latex"))
open("tables/table1_summary_flat.tex", "w").write(buf.getvalue())
1.2 Figure 1 — raw trends / treatment rollout
For DID / event-study designs, the first figure of an applied paper is almost always either (a) raw treated-vs-control means over time, or (b) the staggered rollout heat-strip showing which units are treated when. Both are one-liners:
# (a) Raw trends with vertical line at treatment start (DID Figure 1 style)
fig, ax = sp.parallel_trends_plot(df, y="wage", time="year", treat="training",
treat_time=2015, ci=True,
labels={"treated":"Trained", "control":"Untrained"})
fig.savefig("figures/fig1a_raw_trends.png", dpi=300)
# (b) Treatment rollout heatmap (staggered DID convention; Goodman-Bacon-friendly)
fig, ax = sp.treatment_rollout_plot(df, time="year", treat="training", id="worker_id",
sort_by="first_treat_year",
title="Figure 1. Treatment timing")
fig.savefig("figures/fig1b_rollout.png", dpi=300)
For matching designs, also produce a love plot of standardized differences pre/post matching (Step 3.4).
This is the heart of an AER paper. Before any code, write down the equation explicitly and state the identifying assumption. Vague identification language is the single most common reason a referee rejects an applied paper.
When design="auto" is too opaque, use this decision tree:
┌─ running var + cutoff ───────────────── RDD (sp.rdrobust)
│
├─ exogenous instrument Z ─────────────── IV (sp.ivreg, sp.dml)
data + question ─┤
├─ pre/post × treat/control ─┬ 2 periods ── 2×2 DID (sp.did)
│ └ staggered ── CS / SA (sp.callaway_santanna)
│
├─ 1 treated unit + donor pool + long pre ── SCM (sp.synth, sp.sdid)
│
├─ high-dim X, selection-on-observables ── DML / Causal Forest
│
└─ none of the above ──────────────────── matching + E-value (sp.match, sp.evalue)
2.3 Estimand-first DSL = pre-registration
sp.causal_question declares the five-tuple (population, treatment, outcome, estimand, design) and .identify() picks the estimator with its assumptions written down. Treat the IdentificationPlan as your pre-registration artifact — freeze it before running q.estimate() so the analysis plan is a dated document, not a post-hoc rationalization.
q = sp.causal_question(
treatment="training", outcome="wage", data=df,
population="manufacturing workers, 2010–2020",
estimand="ATT",
design="auto", # 'auto' | 'did' | 'event_study' | 'regression_discontinuity'# | 'iv' | 'rct' | 'selection_on_observables'# | 'synthetic_control' | 'natural_experiment'# | 'policy_shock' | 'longitudinal_observational'
time_structure="panel", time="year", id="worker_id",
covariates=["age", "edu", "tenure"],
)
plan = q.identify() # IdentificationPlan: estimator + assumptions + fallbacksprint(plan.summary()) # human-readable Methods paragraphprint(plan.identification_story) # narrative of why this estimator identifies the estimand# FREEZE the plan to disk BEFORE estimating — this is your pre-registration.# `q` (CausalQuestion) carries the question (population / treatment / outcome).# `plan` (IdentificationPlan) carries the strategy (estimator / story /# assumptions / fallbacks / warnings). The estimating equation is *your*# job to write down — paste it from the §2.1 table that matches plan.estimator.from pathlib import Path
bullets = lambda xs: "\n".join(f"- {x}"for x in xs) if xs else"- (none)"
Path("artifacts/empirical_strategy.md").write_text(
f"# Empirical Strategy (pre-registration)\n\n"f"**Population**: {q.population}\n"f"**Treatment**: `{q.treatment}` **Outcome**: `{q.outcome}`\n"f"**Estimand**: {plan.estimand}\n"f"**Estimator**: `sp.{plan.estimator}`\n\n"f"## Estimating equation (paste from §2.1 row matching `{plan.estimator}`)\n"f"```\n<paste here>\n```\n\n"f"## Identification story\n{plan.identification_story}\n\n"f"## Identifying assumptions (must defend in §2)\n{bullets(plan.assumptions)}\n\n"f"## Auto-flagged warnings\n{bullets(plan.warnings)}\n\n"f"## Fallback estimators (Step 7 robustness)\n{bullets(plan.fallback_estimators)}\n"
)
# Machine-readable sidecar (full question, replayable):
Path("artifacts/causal_question.yaml").write_text(q.to_yaml())
result = q.estimate() # run only after the plan is committed to disk / git
2.5 (Optional) LLM-assisted DAG addendum
Useful when the user wants an explicit DAG to defend in §2 or §7. Pipe the discovered DAG into sp.causal(..., dag=...).
proposal = sp.llm_dag_propose(
variables=df.columns.tolist(),
domain="labor economics: training, wages, tenure",
client=my_llm_client, # .complete(prompt) -> str; None = heuristic
)
validation = sp.llm_dag_validate(proposal, df, alpha=0.05) # (dag, data) positionalprint(validation.edge_evidence)
discovered = sp.llm_dag_constrained(
df,
descriptions={"wage": "monthly wage USD", "training": "0/1 program"},
oracle=my_llm_client.suggest_edges, # optional; falls back to plain PC
max_iter=3,
)
# The result is an LLMConstrainedDAGResult — it has NO `.dag` attribute. Get a DAG with# `.to_dag()` (or inspect `.final_edges`). Pass into Step 4 as:# sp.causal(..., dag=discovered.to_dag())
AER convention: the identification figure precedes the regression table. The reader should see graphical evidence that PT holds / first stage is strong / RD jumps cleanly before you ask them to trust your point estimate.
3.1 Event-study plot + numerical pre-trends test (DID identification)
Pre-period coefficients ≈ 0 (with the −1 reference period normalized to zero) is the visual evidence for parallel trends. Pair the figure with a numerical pre-trends test so reviewers don't have to eyeball it.
# --- The event-study FIGURE comes from a Callaway–Sant'Anna (or sun_abraham)# result, NOT from sp.event_study(). The figure plotters# (enhanced_event_study_plot / cs.plot() / ggdid / group_time_plot) consume a# CS/SA result; feeding them sp.event_study() output raises KeyError('att').# Use `x=` for covariates (NOT `covariates=` — that kwarg does not exist on CS).
cs = sp.callaway_santanna(df, y="wage", g="first_treat_year",
t="year", i="worker_id",
x=["age", "edu"])
# Figure 2a — dynamic ATT / event-study coefficient plot. Plotters return (fig, ax).
fig, ax = sp.enhanced_event_study_plot(
cs, shade_pre=True,
title="Figure 2a. Event-study coefficients (95% CI; ref. period = −1)")
fig.savefig("figures/fig2a_event_study.png", dpi=300)
# (equivalently: `fig, ax = cs.plot()` or `fig, ax = sp.ggdid(cs)` /# `fig, ax = sp.group_time_plot(cs)` — all consume the CS result and return (fig, ax).)# Numerical pre-trends test (Roth 2022 power) for the table footnote. THIS is what# sp.event_study() is for — the coefficient/pre-trend numerics, not the figure.
es = sp.event_study(df, y="wage", treat_time="first_treat_year",
time="year", unit="worker_id",
window=(-4, 4), ref_period=-1,
covariates=["age", "edu"])
print(sp.pretrends_summary(es)) # F-stat, p-value, max-PT bound# es.model_info["pretrend_test"] holds the same numbers machine-readably.# Bacon decomposition figure for staggered DID (Figure 2a-bis)
bd = sp.bacon_decomposition(df, y="wage", treat="training",
time="year", id="worker_id")
fig, ax = sp.bacon_plot(bd, title="Figure 2a-bis. Goodman-Bacon weights")
fig.savefig("figures/fig2a2_bacon.png", dpi=300)
# Borusyak–Jaravel–Spiess joint pre-trends test — needs the CS/SA result# AND the underlying panel (NOT the event_study() output):
sp.bjs_pretrend_joint(cs, df, y="wage", group="first_treat_year",
time="year", first_treat="first_treat_year",
controls=["age", "edu"])
3.3 RD: McCrary density + canonical RD plot + binscatter
The signature RD figure is sp.rdplot (CCT-style binned scatter with local-polynomial fit on each side), paired with the McCrary manipulation test. Together they answer: (a) is there a visual jump? (b) is the density continuous at the cutoff?
# Figure 2b — canonical RD plot (binned means + local poly fit on each side)
fig, ax = sp.rdplot(df, y="y", x="running_var", c=0,
p=4, kernel="triangular", binselect="esmv",
shade_ci=True, ci_level=0.95)
fig.savefig("figures/fig2b_rdplot.png", dpi=300)
# Figure 2b-bis — McCrary density (manipulation test). .plot() → (fig, ax)
fig, ax = sp.rddensity(df, x="running_var", c=0).plot()
fig.savefig("figures/fig2b2_mccrary.png", dpi=300)
# Optional: covariate-adjusted binscatter (continuity in covariates is also testable)
fig, ax, _ = sp.binscatter(df, y="age", x="running_var", n_bins=40, ci=True)
fig.savefig("figures/fig2b3_cov_binscatter.png", dpi=300)
3.4 Matching: love plot (standardized differences)
m = sp.match(df, y="wage", treat="training",
covariates=["age", "edu", "tenure"], method="nearest")
fig, ax = m.plot() # |std diff| pre vs post; target |Δ|<0.1
fig.savefig("figures/fig2c_love_plot.png", dpi=300)
3.5 SCM: synthetic-control trajectory + gap plot
For synthetic-control designs the canonical Figure 2 is the treated-vs-synthetic time-series with treatment time annotated. synthdid_plot does this in one line.
Identification-specific checks (PT for DID, weak-IV F, density for RD, common support for matching) are also auto-run inside sp.causal(...) in Step 4 — don't duplicate the numerics here, but DO produce the figures: a referee scans the figures first.
Step 4 — Main results (multi-regression tables, AER style)
This is the densest section of an applied paper. A modern AER §4 typically contains 2–3 multi-regression tables and one coefficient plot:
Table 2 (main): progressive controls, 4–6 columns
Table 2-bis (design horse race): same coefficient under OLS / 2SLS / DID / DML
Table 2-ter (multi-outcome): same treatment, several outcomes side-by-side
Figure 3 (coefplot): visual summary of β̂ and 95% CI across specs
Estimator routing (memorize this — getting it wrong silently produces nonsense):
No FE → sp.regress("y ~ x1 + x2", df, cluster="firm_id")
DID / event-study → sp.callaway_santanna(...) / sp.sun_abraham(...)
Never write sp.regress("y ~ x | firm_id") — sp.regress does not parse | and silently treats x | firm_id as a single variable name. Use sp.feols for any formula containing |.
sp.regtable(*models, ...) is the workhorse. Useful kwargs:
keep : list of coef names to display (e.g. ["training"])
drop : list of coef names to suppress (controls)
model_labels : column labels ["(1) Baseline", "(2) +Demog", ...]
dep_var_labels : dep-var-row labels (for multi-outcome tables)
panel_labels : panel-A / panel-B layout for stacked tables
coef_labels : pretty-print names for coefficients
stars : "aer" → * 0.10 ** 0.05 *** 0.01 (or "default", "none")
stats : footer rows ["N","R2","Cluster","FE","DV mean", ...]
output : "latex" | "html" | "markdown" | "text"
filename : path to write the table
4.1 Pattern A — Progressive controls (the canonical Table 2)
Stable β̂ across columns ⇒ less concern that selection on observables is driving the estimate (Oster 2019 selection-stability logic; quantified in Step 7.5). sp.regtable(*models) is the StatsPAI equivalent of Stata outreg2 / esttab and R modelsummary::msummary / summary_col — it consolidates N models into ONE table with one column per model.
(1) Baseline
(2) +Demographics
(3) +Labor-market
(4) +Region×Industry FE
(5) +Worker FE
Controls
none
age, edu
+ tenure, firm_size
high-dim FE
individual FE
# RULE: pure OLS → sp.regress; high-dim FE absorption → sp.feols# (sp.regress does NOT parse `|` as FE — it's a thin OLS wrapper. Use# `sp.feols("y ~ x | fe1 + fe2", df, vcov={"CRV1":"firm_id"})` for FE.)
M1 = sp.regress("wage ~ training", df, cluster="firm_id")
M2 = sp.regress("wage ~ training + age + edu", df, cluster="firm_id")
M3 = sp.regress("wage ~ training + age + edu + tenure + firm_size", df, cluster="firm_id")
M4 = sp.feols ("wage ~ training + age + edu + tenure + firm_size | region + industry + year",
df, vcov={"CRV1": "firm_id"})
M5 = sp.feols ("wage ~ training + age + edu + tenure + firm_size | worker_id + year",
df, vcov={"CRV1": "firm_id"})
# Consolidate 5 models into ONE table (= Stata `outreg2 [M1..M5] using ..., replace`).# **Default = show ALL coefficients verbatim — controls AND the intercept**# (AER convention; readers verify the full spec). Pass NO `keep=`/`drop=` and# `regtable` will surface every estimated parameter. Add `drop=["Intercept"]`# only if you want to suppress the constant for paper aesthetics; add# `keep=[focal]` only when a focal-coefficient-only table is intentional.
rt = sp.regtable(M1, M2, M3, M4, M5,
template="aer", # auto-applies SE label, star levels, font
coef_labels={"training": "Job training"},
model_labels=["(1) Baseline", "(2) +Demog.", "(3) +Labor-mkt",
"(4) Region×Ind. FE", "(5) Worker FE"],
stats=["N", "R2", "Cluster", "FE", "DV mean"],
title="Table 2. Effect of training on wages")
# Variants (all opt-in — the default above is preferred):# • drop intercept only: sp.regtable(..., drop=["Intercept"])# • focal-coefficient only: sp.regtable(..., keep=["training"])# • mixed-magnitude table: sp.regtable(..., fmt="auto")# Use whenever a single table mixes dollar-magnitude coefficients# (e.g. earnings ≈ 1500) with elasticity-magnitude coefficients# (e.g. log-earnings ≈ 0.09). The default fmt="%.3f" pads the dollar# side; a fixed fmt="%.0f" rounds the elasticity side to "0" while# significance stars survive — the silent LaLonde-style precision# trap. fmt="auto" picks per-value precision: thousands separator# for |β|≥1000, integer for ≥100, 1 dp for ≥10, 2 dp for ≥1, 3 dp# below — so neither magnitude is killed.# Export to ALL THREE in three lines — Word for co-authors, Excel for editors, LaTeX for build:
rt.to_word ("tables/table2_main.docx")
rt.to_excel("tables/table2_main.xlsx")
open("tables/table2_main.tex", "w").write(rt.to_latex())
4.2 Pattern B — Design horse race (Table 2-bis)
Show the same coefficient of interest under multiple identification strategies. This is the AER credibility move: convergent evidence across designs each making different identifying assumptions.
4.3 Pattern C — Multi-outcome table (same X, several Y's)
A single treatment, several outcomes. Use dep_var_labels so each column carries the Y name.
ys = ["wage", "log_wage", "weeks_employed", "left_firm", "promoted"]
multi_y = [sp.feols(f"{y} ~ training + age + edu + tenure | industry + year",
df, vcov={"CRV1": "firm_id"})
for y in ys]
rt = sp.regtable(*multi_y,
template="aer",
dep_var_labels=ys, # column header: dep var
model_labels=["(1)","(2)","(3)","(4)","(5)"],
stats=["N","R2","DV mean","Cluster"],
title="Table 2-ter. Effect of training on multiple outcomes")
rt.to_word ("tables/table2c_multi_outcome.docx")
rt.to_excel("tables/table2c_multi_outcome.xlsx")
open("tables/table2c_multi_outcome.tex", "w").write(rt.to_latex())
4.4 Pattern D — Stacked Panel A / Panel B table
Same model family, two horizons (short-run / long-run) or two samples (pre-2015 / post-2015) stacked vertically. Use panel_labels.
panelA = [sp.feols("wage_t1 ~ training + X | industry + year", df, vcov={"CRV1":"firm_id"}),
sp.feols("wage_t1 ~ training + X | worker_id + year", df, vcov={"CRV1":"firm_id"})]
panelB = [sp.feols("wage_t5 ~ training + X | industry + year", df, vcov={"CRV1":"firm_id"}),
sp.feols("wage_t5 ~ training + X | worker_id + year", df, vcov={"CRV1":"firm_id"})]
rt = sp.regtable(*panelA, *panelB,
template="aer",
panel_labels=["Panel A. Short-run (1 year)",
"Panel A. Short-run (1 year)",
"Panel B. Long-run (5 years)",
"Panel B. Long-run (5 years)"],
model_labels=["(1) Industry FE","(2) Worker FE"]*2,
stats=["N","R2"],
title="Table 2-quater. Short- vs long-run effects")
rt.to_word ("tables/table2d_horizons.docx")
rt.to_excel("tables/table2d_horizons.xlsx")
open("tables/table2d_horizons.tex", "w").write(rt.to_latex())
4.5 Pattern E — IV reporting triplet (first-stage / reduced-form / 2SLS)
The textbook AER IV table presents the first stage, the reduced form, and the 2SLS in three columns so the reader can verify Wald-ratio = RF / FS.
Trap:sp.ivreg does not absorb | fe and does not parse C(fe) — it silently drops a | industry + year term (identical β̂ with or without it), so a 2SLS column written that way would not control for the FE the first-stage/reduced-form columns absorb. Keep the IV triplet on the same low-dim control set in all three columns; to control for fixed effects in a 2SLS, pre-build dummy columns in pandas and add them explicitly, or partial the FE out first.
fs = sp.feols("training ~ Z + age + edu", df, vcov={"CRV1":"firm_id"}) # 1st stage
rf = sp.feols("wage ~ Z + age + edu", df, vcov={"CRV1":"firm_id"}) # reduced form
iv = sp.ivreg("wage ~ (training ~ Z) + age + edu", df, cluster="firm_id") # 2SLS (same controls)
rt = sp.regtable(fs, rf, iv,
template="aer",
keep=["Z", "training"], # IV triplet is intentionally focal:# show only Z + endog so the reader can# eyeball Wald-ratio = RF / FS. For the# full coef list, drop the kwarg entirely.
dep_var_labels=["training", "wage", "wage"],
model_labels=["(1) First stage", "(2) Reduced form", "(3) 2SLS"],
stats=["First-stage F", "N", "R2", "Cluster"],
title="Table 2-quinto. IV reporting triplet")
rt.to_word ("tables/table2e_iv_triplet.docx")
rt.to_excel("tables/table2e_iv_triplet.xlsx")
open("tables/table2e_iv_triplet.tex", "w").write(rt.to_latex())
4.6 Pattern F — Causal-design main via sp.causal(...)
For DID / IV / RD / SCM mains, the sp.causal(...) orchestrator returns a CausalResult plus diagnostics and an automatic robustness preview. Pipe .result into regtable:
w = sp.causal(df, y="wage", treatment="training",
id="worker_id", time="year", design="did",
covariates=["age", "edu", "tenure"],
dag=discovered.to_dag()) # optional (LLMConstrainedDAGResult.to_dag())print(w.diagnostics) # PT verdict + warningsprint(w.recommendation) # which estimator + whyprint(w.result.summary()) # point estimate + cluster-robust SE + CIprint(w.robustness_findings) # automated robustness battery preview
4.7 Figure 3 — coefficient plot of the main table
Replace one of the wall-of-numbers tables with a coefplot in the body, push the table to the appendix. Modern AER papers increasingly do this.
fig, ax = sp.coefplot(M1, M2, M3, M4, M5,
model_names=["(1)","(2)","(3)","(4)","(5)"],
variables=["training"],
title="Figure 3. β̂ on training across specifications (95% CI)",
alpha=0.05)
fig.savefig("figures/fig3_coefplot.png", dpi=300)
Reporting checklist for the Table 2 footnote (AER house style)
Fixed-effects absorbed — regtable auto-adds one footer row per FE name (e.g. Industry FE: Yes / Year FE: Yes / Worker_id FE: No) whenever any column comes from sp.feols(... | fe1 + fe2 ...). Don't hand-roll these rows.
Sample size and number of clusters
Estimator (OLS / 2SLS / CS-DID / SCM / DML)
Stars convention * 0.10 ** 0.05 *** 0.01
Mean of dependent variable in the estimation sample (so β̂ can be read as a % of the base rate)
Step 5 — Heterogeneity (Table 3 + Figure 4)
The AER §5 Heterogeneity combines (a) a subgroup regression table with one column per subgroup (binary moderators + interaction terms), and (b) a CATE / dose-response figure for continuous moderators. Both should appear; they answer different questions.
5.1 Pattern G — Subgroup regtable (Table 3)
One column per subgroup, with the same specification re-run on each slice. Clean, easy to read, expected by referees.
slices = {
"(1) All": df,
"(2) Female": df[df["female"] == 1],
"(3) Male": df[df["female"] == 0],
"(4) Low skill": df[df["skill_quartile"].isin([1, 2])],
"(5) High skill": df[df["skill_quartile"].isin([3, 4])],
"(6) Small firm": df[df["firm_size"] < 100],
"(7) Large firm": df[df["firm_size"] >= 100],
}
gmodels = [sp.feols("wage ~ training + age + edu + tenure | industry + year",
d, vcov={"CRV1": "firm_id"}) for d in slices.values()]
rt = sp.regtable(*gmodels,
template="aer",
coef_labels={"training": "Training"},
model_labels=list(slices),
stats=["N","R2","DV mean"],
title="Table 3. Heterogeneous effects of training")
rt.to_word ("tables/table3_heterogeneity.docx")
rt.to_excel("tables/table3_heterogeneity.xlsx")
open("tables/table3_heterogeneity.tex", "w").write(rt.to_latex())
Test moderation formally with interaction terms — referees often ask whether the gap between subgroups is statistically significant, which requires the interaction p-value.
H1 = sp.feols("wage ~ training*female + age + edu + tenure | industry + year",
df, vcov={"CRV1": "firm_id"})
H2 = sp.feols("wage ~ training*C(skill_quartile) + age + edu + tenure | industry + year",
df, vcov={"CRV1": "firm_id"})
H3 = sp.feols("wage ~ training*log_firm_size + age + edu + tenure | industry + year",
df, vcov={"CRV1": "firm_id"})
rt = sp.regtable(H1, H2, H3,
template="aer",
keep=["training", "training:female", # interaction-form heterogeneity"training:C(skill_quartile)[T.2]", # is intentionally focal:"training:C(skill_quartile)[T.3]", # only the main effect + interactions"training:C(skill_quartile)[T.4]", # are reported. Drop this kwarg"training:log_firm_size"], # entirely to show full controls.
model_labels=["(1) ×Female", "(2) ×Skill quartile", "(3) ×log(Firm size)"],
stats=["N","R2"],
title="Table 3-bis. Interaction-form heterogeneity")
rt.to_word ("tables/table3b_interactions.docx")
rt.to_excel("tables/table3b_interactions.xlsx")
open("tables/table3b_interactions.tex", "w").write(rt.to_latex())
5.4 Figure 4-bis — CATE distribution (DR-Learner / causal forest)
The CATE plotters read per-row conditional effects out of the result's
model_info["cate"] array. There is no .cate_estimates attribute — the raw
per-row CATE vector lives at ml.model_info["cate"] (an ndarray of length n),
and summary stats at model_info["cate_mean"] / cate_q25 / cate_q75 / ....
sp.causal_forest returns a summary result that does not populate
model_info["cate"], so for the CATE histogram and grouped bar chart use a
meta-learner (or any DR-/X-/R-learner) and pass its result to the plotters.
ml = sp.metalearner(df, y="wage", treat="training",
covariates=["age","edu","tenure","firm_size"], learner="dr")
# Raw per-row CATE vector (if you need the numbers, not just the figure):
cate_i = ml.model_info["cate"] # ndarray, length n (NOT ml.cate_estimates)
fig, ax = sp.cate_plot(ml, kind="hist",
title="Figure 4b. Distribution of conditional ATE")
fig.savefig("figures/fig4b_cate_hist.png", dpi=300)
# CATE by group bar chart: first compute the group-level table, THEN plot it.# `cate_group_plot` takes a DataFrame (from cate_by_group), not the result object.
g = sp.cate_by_group(ml, df, by="skill_quartile", n_groups=4)
fig, ax = sp.cate_group_plot(g, title="Figure 4c. CATE by skill quartile")
fig.savefig("figures/fig4c_cate_by_group.png", dpi=300)
# Tabular summary for the appendixprint(sp.cate_summary(ml))
print(g) # group-level CATE table
5.5 Subgroup-analysis dispatcher (one-liner)
sp.subgroup_analysis(df, formula="wage ~ training + age + edu + tenure",
x="training",
by={"gender": "female", "skill": "skill_quartile"},
robust="hc1") # quick subgroup β̂ table (HC1 by default; no cluster arg)
For continuous moderators or many subgroups, prefer:
sp.causal_forest(formula="wage ~ training | X", data=df) — CATE summary only (does not populate model_info["cate"]; use a meta-learner for per-row CATEs)
Step 7 — Robustness gauntlet (the AER referee gauntlet)
The seven canonical robustness blocks of an applied paper. A modern AER paper expects most of these in the body or appendix — assemble a Table A1-style robustness panel from the outputs.
7.1 Placebo tests
sp.rdplacebo(df, y="y", x="running_var", c=0,
placebo_cutoffs=[-2, -1, 1, 2]) # RD: fake cutoffs
sp.synth_time_placebo(df, outcome="y", unit="unit", time="time",
treated_unit=1, treatment_time=2000,
n_placebo_times=10) # SCM in-time placebo
sp.synthdid_placebo(...) # SDID placebo# For DID: re-run with a fake treat year before actual treatment and confirm β̂ ≈ 0.
Cluster-level choice is itself a robustness check — show the result is not driven by an over-narrow cluster.
# For statsmodels-backed sp.regress / sp.ivreg results:
sp.twoway_cluster(M3, df, cluster1="firm_id", cluster2="year") # two-way clustering
sp.conley(M3, df, lat="lat", lon="lon",
dist_cutoff=100, kernel="uniform") # spatial HAC (Conley 1999)# For pyfixest-backed sp.feols results, set 2-way cluster directly in `vcov`:
sp.feols("y ~ x | firm_id + year", df,
vcov={"CRV1": "firm_id+year"}) # 2-way: firm × year
7.5 Oster (2019) selection bound
"How big would unobserved selection have to be for β to flip sign / vanish?" The Oster δ tells you whether the bound on selection on unobservables, relative to selection on observables, has to exceed an implausible value to overturn the result.
7.6 Honest DID — Rambachan–Roth (2023) PT sensitivity
honest_did only consumes a CS / SA / did_multiplegt event-study result
(or aggte(result, type='dynamic')). Pass the cs object built in §3.1,
not a generic OLS/FE main-table result:
sp.honest_did(cs, method="smoothness") # bound β under bounded PT violation
7.11 Pattern H — Robustness master table (Table A1, one row per check)
The canonical AER appendix Table A1 stacks every robustness specification next to the baseline so reviewers see at a glance that β̂ survives. sp.regtable accepts any mix of EconometricResults / CausalResult, so build the list dynamically:
baseline = sp.feols("wage ~ training + age + edu + tenure | industry + year",
df, vcov={"CRV1": "firm_id"})
rob = {
"(1) Baseline": baseline,
"(2) Drop top 1% wage": sp.feols("wage ~ training + age + edu + tenure | industry + year",
df.query("wage < wage.quantile(0.99)"),
vcov={"CRV1": "firm_id"}),
"(3) Balanced panel": sp.feols("wage ~ training + age + edu + tenure | industry + year",
sp.balance_panel(df, entity="worker_id", time="year"),
vcov={"CRV1": "firm_id"}),
"(4) Drop early cohorts": sp.feols("wage ~ training + age + edu + tenure | industry + year",
df.query("first_treat_year > 2008"),
vcov={"CRV1": "firm_id"}),
"(5) Worker FE": sp.feols("wage ~ training + age + edu + tenure | worker_id + year",
df, vcov={"CRV1": "firm_id"}),
"(6) 2-way cluster": sp.feols("wage ~ training + age + edu + tenure | industry + year",
df, vcov={"CRV1": "firm_id+year"}), # 2-way: firm × year# sp.conley needs a STATSMODELS-backed result (sp.regress/sp.ivreg) — it raises# KeyError on a pyfixest feols result. Re-fit the spec via sp.regress for this row."(7) Conley spatial SE": sp.conley(sp.regress("wage ~ training + age + edu + tenure",
df, cluster="firm_id"),
df, lat="lat", lon="lon", dist_cutoff=100),
"(8) Log outcome": sp.feols("log_wage ~ training + age + edu + tenure | industry + year",
df, vcov={"CRV1": "firm_id"}),
"(9) IHS outcome": sp.feols("ihs_wage ~ training + age + edu + tenure | industry + year",
df, vcov={"CRV1": "firm_id"}),
"(10) PSM-weighted": sp.match(df, y="wage", treat="training",
covariates=["age","edu","tenure","firm_size"],
method="nearest"),
"(11) Entropy balance": sp.ebalance(df, y="wage", treat="training",
covariates=["age","edu","tenure","firm_size"]),
"(12) DML-PLR": sp.dml(df, y="wage", treat="training",
covariates=["age","edu","tenure","firm_size"], model="plr"),
}
# Robustness master = AER Table A1 — readers MUST see every coefficient# across every spec to verify nothing is hiding behind `keep=`. Default to# the full coef table (intercept included); only switch to# `keep=["training"]` if a referee has explicitly asked for a focal-only# summary, or add `drop=["Intercept"]` if you want the constant suppressed.
rt = sp.regtable(*rob.values(),
template="aer",
coef_labels={"training": "Training (β̂)"},
model_labels=list(rob),
stats=["N", "R2", "Cluster", "FE"],
title="Table A1. Robustness of the main estimate")
rt.to_word ("tables/tableA1_robustness.docx")
rt.to_excel("tables/tableA1_robustness.xlsx")
open("tables/tableA1_robustness.tex", "w").write(rt.to_latex())
# Equivalent one-shot via the paper-format multi-panel API — produces a# single .docx / .xlsx that you can hand a co-author, with main + robustness# (+ heterogeneity / placebo if you have them) auto-laid-out per AER style:
sp.paper_tables(main=[M1, M2, M3, M4, M5],
robustness=list(rob.values()),
template="aer",
coef_labels={"training": "Training"},
model_labels_main=["(1)","(2)","(3)","(4)","(5)"],
model_labels_robustness=list(rob),
# paper_tables only accepts `keep=`, not `drop=`. Omit both to# show every coefficient (AER convention). Pass `keep=["training"]`# only when a focal-only summary is desired.
).to_docx("tables/paper_tables.docx")
7.12 Figure 5 — coefficient forest plot of all robustness specs
