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: Wang & Rozelle (2026), Journal of Open Source Software 11(125), 10604, https://doi.org/10.21105/joss.10604; 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.

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:

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>

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

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:

How to invoke a non-default mode (Claude / agent picks this up from the user's wording):

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):

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:

Journal templates (apply the right SE label, star levels, and notes automatically):

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):

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

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())

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:

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

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.

0.1 Sample-construction log (footnote 4)

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)

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:

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.

# 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.

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).

Step 2 — Empirical strategy (Section "Identification")

This is the heart of an AER paper. Bef

…(truncated)