# StatsPAI_skill > 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 - Skill: `brycewang-stanford/statspai-skill` (Agent Skill, multi-file: 2 files) - Install (CLI): `npx skillmds add brycewang-stanford/statspai-skill` - Raw SKILL.md: https://api.skillmd.com/api/skills/brycewang-stanford/statspai-skill/raw - Safety review: pending (external: skill-scanner PASS, skillspector WARNING) - Works with: Claude Code, Claude.ai, OpenAI Codex - Category: Product & Planning - Author: brycewang-stanford (https://skillmd.com/u/brycewang-stanford) - Updated: 2026-08-19 - Page: https://skillmd.com/skills/brycewang-stanford/statspai-skill --- # StatsPAI: Agent-Native Causal Inference & AER-Style Empirical Workflow 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). - **Source**: https://github.com/brycewang-stanford/StatsPAI - **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 | > | `sp.dragonnet` / `sp.tarnet` / `sp.cfrnet` / `sp.cevae` (neural causal) | `neural` (torch) | `pip install "statspai[neural]"` | `ImportError: PyTorch is required …` | > | `sp.causal_text.*` (text-as-treatment) | `text` (sentence-transformers) | `pip install "statspai[text]"` | `ImportError` on embed | > > 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. ```python 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 1. **Self-describing**: `sp.list_functions()` / `sp.describe_function(name)` / `sp.function_schema(name)` — registered symbols are discoverable without doc lookup. 2. **Structured results**: mature estimators return result objects with methods such as `.summary()`, `.plot()`, `.diagnostics`, `.to_latex()`, `.to_word()`, `.cite()` when supported. 3. **One import, full pipeline**: data contract → Table 1 → estimand-first DSL → identification graphs → main table → heterogeneity → mechanisms → robustness → replication package. 4. **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: 1. **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". 2. **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. 3. **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. 4. **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. 5. **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. 6. **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: ```text best_skill: train_signal: selection_split: heldout_gate: accepted_patterns: rejected_patterns: patch_scope: reject_if: ``` 1. **Route card**: record the mode (`econ`, `epi`, or `ml-causal`), estimand, identification design, focal outcome/treatment, StatsPAI install extras, and required artifacts. 2. **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. 3. **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. 4. **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. 5. **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. 6. **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. 7. **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. ``` Paper section Step StatsPAI moves ─────────────────────────── ───── ──────────────────────────────────────────────── Pre-Analysis Plan −1 sp.power.* + freeze IdentificationPlan to disk §1. Data 0 data_contract + sample-construction log (footnote 4) §1.1 Descriptives (Table 1) 1 sp.sumstats · sp.balance_table · sp.describe §2. Empirical Strategy 2 write equation + identifying assumption + sp.causal_question (LLM-DAG addendum) 2.5 sp.llm_dag_propose · validate · constrained §3. Identification graphics 3 event-study · first-stage F · McCrary · love plot §4. Main Results (Table 2) 4 progressive controls + FE (sp.regtable / sp.causal) §5. Heterogeneity (Table 3) 5 sp.subgroup_analysis · sp.continuous_did · CATE §6. Mechanisms 6 sp.mediation · sp.decompose §7. Robustness gauntlet 7 placebo · Oster · honest_did · E-value · 2-way / Conley SE · spec_curve §8. Replication package 8 .to_latex() · .plot() · reproducibility stamp ``` > **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. ## Three domain modes (default = AER econ; alternates = epi & ML-causal) 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: | Mode | Reader convention | Identification stack | Reporting stack | Jump to | |---|---|---|---|---| | **Default — Applied Econ (AER / QJE / AEJ)** | "Show the equation + identifying assumption + design horse-race; controls visible; clustered SE" | DID / IV / RD / SCM / matching / `feols` HDFE | AER house-style multi-column `regtable` + 8-section paper layout | §−1 → §8 (entire playbook above) | | **Mode A — Epidemiology / Public Health** | "STROBE / TRIPOD-AI; target trial protocol; doubly-robust estimand; absolute & relative risk; KM survival" | Target-trial emulation · IPTW · g-formula · TMLE · Mendelian randomization · KM/AFT | Same `regtable` + `collect`, with risk-difference / hazard-ratio / E-value rows | §A. Epidemiology pipeline | | **Mode B — ML Causal Inference** | "DML / meta-learners / causal forest / DR-learner; CATE distribution; policy value" | DML · S/T/X/R/DR-Learner · GRF causal forest · Dragonnet/TARNet/CEVAE · BCF · matrix completion | `regtable` ML horse-race + `cate_plot` + policy-value table + `conformal_causal` PI | §B. ML causal pipeline | **How to invoke a non-default mode** (Claude / agent picks this up from the user's wording): | User says... | Mode the skill switches to | |---|---| | "Run a DID / IV / RD / event study", "AER table", "applied micro" | Default (AER econ) | | "Target trial emulation", "g-formula", "IPTW", "TMLE", "Mendelian randomization", "STROBE / TRIPOD", "公共健康 / 流行病学", "epi pipeline", "RWE study", "cohort study", "case-control" | Mode A (Epi) | | "DML", "double machine learning", "causal forest", "meta-learner", "CATE", "Dragonnet", "BCF", "policy learning", "conformal causal", "ML causal", "uplift modeling", "因果机器学习" | Mode B (ML causal) | | "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): | § | Artifact | StatsPAI primitive | Filenames (write all three) | |---|---|---|---| | §1 | **Figure 1**: raw trends / treatment rollout | `sp.parallel_trends_plot` · `sp.treatment_rollout_plot` | `figures/fig1_trends.png` | | §1 | **Table 1**: summary stats (full / treated / control + Δ) | `sp.sumstats` + `sp.mean_comparison(...).to_word()/.to_excel()` (or `sp.collect().add_summary().add_balance()`) | `tables/table1_summary.{tex,docx,xlsx}` | | §3 | **Figure 2**: identification graphic (event-study / first-stage / McCrary / RD scatter / SCM trajectory) | `sp.enhanced_event_study_plot` · `sp.binscatter` · `sp.rdplot` · `sp.rddensity().plot()` · `sp.synthdid_plot` | `figures/fig2_identification.png` | | §4 | **Table 2**: main results — progressive controls | `rt = sp.regtable(M1...M5, template="aer"); rt.to_word(...); rt.to_excel(...)` | `tables/table2_main.{tex,docx,xlsx}` | | §4 | **Table 2-bis**: design horse-race (OLS / IV / DID / DML) | `sp.regtable(ols, iv, did, dml, ...).to_word/.to_excel` | `tables/table2b_designs.{tex,docx,xlsx}` | | §4 | **Figure 3** (optional): coefficient plot across specs | `sp.coefplot(M1, M2, M3, M4)` | `figures/fig3_coef.png` | | §5 | **Table 3**: heterogeneity by subgroup | `sp.regtable(g_full, g_male, g_fem, g_q1...q4).to_word/.to_excel` | `tables/table3_heterogeneity.{tex,docx,xlsx}` | | §5 | **Figure 4**: dose-response / CATE | `sp.dose_response(...).plot()` · `sp.cate_plot` · `sp.cate_group_plot` | `figures/fig4_cate.png` | | §6 | **Table 4**: mechanisms (mediation / decomposition) | `sp.regtable(total, direct, indirect).to_word/.to_excel` | `tables/table4_mechanisms.{tex,docx,xlsx}` | | §7 | **Table A1**: robustness master (one row per check) | `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 | `rt = sp.regtable(M1, M2, ..., template="aer", title=...)` *(default: all coefs incl. intercept)*
`rt.to_word("table2.docx")`
`rt.to_excel("table2.xlsx")`
`rt.to_latex()` · `rt.to_markdown()` | `template`, `coef_labels`, `model_labels`, `panel_labels`, `dep_var_labels`, `stats`, `stars`, `add_rows`; opt-in filters: `drop=["Intercept"]` (suppress constant), `keep=[focal]` (focal-only) | | **2. Multi-panel paper format** (Tables 2 + 3 + A1 + A2 in one file) | Producing the *paper-tables block* — main + heterogeneity + robustness + placebo as a single document | `pt = sp.paper_tables(main=[M1...M5], heterogeneity=[H1,H2,H3], robustness=[R1...Rn], placebo=[P1,P2], template="aer")`
`pt.to_docx("paper_tables.docx")`
`pt.to_xlsx("paper_tables.xlsx")`
`pt.to_latex(...)` | `main`, `heterogeneity`, `robustness`, `placebo`, `template`, `coef_labels`, `model_labels_`, `keep` | | **3. Full session bundle** (Stata 15 `collect` equivalent) | Replication appendix that mixes summary stats + balance + multiple regression tables + headings + prose in **one** file | `c = sp.collect("Paper title", template="aer")`
`c.add_heading("§1. Descriptives")`
`c.add_summary(df, vars=...)`
`c.add_balance(df, treatment=, variables=...)`
`c.add_regression(M1, M2, ..., title="Table 2")`
`c.add_text("Notes ...")`
`c.save("paper.docx")` (auto-detect by extension; `.xlsx`/`.tex`/`.md`/`.html`/`.txt` all work) | `add_heading(level=)`, `add_summary(stats=, labels=)`, `add_balance(weights=, test=)`, `add_regression(**regtable_kwargs)`, `add_table(result)`, `add_text(...)` | **Journal templates** (apply the right SE label, star levels, and notes automatically): ```python sp.list_journal_templates() # → ('aer', 'qje', 'econometrica', 'restat', 'jf', 'aeja', 'jpe', 'restud') rt = sp.regtable(M1, M2, M3, template="qje") # QJE styling; default = full coef list (incl. intercept) rt.to_word("table2_qje.docx") # Opt-in filters: # • drop the constant only: sp.regtable(M1, M2, M3, template="qje", drop=["Intercept"]) # • focal-coefficient only: sp.regtable(M1, M2, M3, template="qje", keep=["x"]) sp.get_journal_template("aer") # inspect a preset # → {'label': 'American Economic Review', 'star_levels': (0.1, 0.05, 0.01), # 'se_label': 'Standard errors', 'stats': ('N', 'R-squared'), # 'notes_default': ('Standard errors in parentheses.', '*** p<0.01, ** p<0.05, * p<0.10.'), # 'font_name': 'Times New Roman'} # note: tuples, not lists ``` **Inline citations in prose** (drop a coefficient straight into a sentence): ```python sp.cite(M3, "training") # → "1.239*** (0.153)" sp.cite(M3, "training", output="latex") # → "1.239^{***}~(0.153)" (wrap in $...$ yourself) ``` > **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: 1. **CJK labels render as ▢▢▢ tofu** — the matplotlib default `DejaVu Sans` carries no Chinese / Japanese / Korean glyphs, so `ax.set_title("教育回报")` silently degrades into squares. 2. **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 ```python import matplotlib as mpl import matplotlib.pyplot as plt def setup_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 is not None: 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()`) ```python fig, ax = plt.subplots(figsize=(4, 2.5)) ax.plot([0, 1, 2], [-1, 0, 1]) ax.set_title("中文标题测试 — Card (1995) 教育回报") ax.set_xlabel("受教育年数 (years)") fig.tight_layout() fig.savefig("figures/_font_smoke_test.png", dpi=300) # delete after verifying ``` 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: > > ```python > fig, ax = sp.parallel_trends_plot(df, y="wage", time="year", treat="training", treat_time=2015) > fig.savefig("figures/fig1.png", dpi=300) > ``` > > Two exceptions to memorize: > - **`sp.binscatter(...)` returns a 3-tuple** `(fig, ax, binned_df)` — `fig, ax, _ = sp.binscatter(...)`. > - **`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 inline` *before* `setup_plot()`. The saved `.png` (driven by `savefig.dpi=300`) is sharp regardless. | | `sp.(...)` 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`. ```python # Always go through the dispatcher when you want auto-solve. The # `sp.power_` 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. ### 0.1 Sample-construction log (footnote 4) ```python sample_log = [] df0 = df_raw.copy(); sample_log.append(("0. raw", len(df0))) df1 = df0.dropna(subset=["wage"]); sample_log.append(("1. drop missing wage", len(df1))) df2 = df1[df1["age"].between(18, 65)]; sample_log.append(("2. drop age outside 18-65", len(df2))) df3 = df2[df2["industry"].isin(MANUF_CODES)]; sample_log.append(("3. keep manufacturing", len(df3))) df = df3 import json; json.dump(sample_log, open("artifacts/sample_construction.json", "w"), indent=2) ``` 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) ```python import pandas as pd, numpy as np, statspai as sp def data_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] if id and 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, } if id and 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" break return 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" assert all(v == 0 for 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: ```python import json; json.dump(contract, open("artifacts/data_contract.json", "w"), indent=2, default=str) ``` ## Step 1 — Descriptive statistics (Table 1) 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. ```python # 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 Excel open("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. ```python 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: ```python # (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). ## Step 2 — Empirical strategy (Section "Identification") This is the heart of an AER paper. Before any code, **write down the equation explicitly** and **state the identifying a …(truncated)