Classical end-to-end empirical analysis workflow in the traditional Python econometric stack — pandas + numpy + scipy + statsmodels + linearmodels + pyfixest + rdrobust + econml + causalml + matplotlib/seaborn. **Defaults to economics empirical-paper style** (AER / QJE / AEJ) — every run produces a publication-ready output set with a multi-column regression table (M1→M6 progressive controls/FE) as the centerpiece, plus Table 1 (descriptives), mechanism / heterogeneity / robustness tables, and event-study + coefficient + trend figures. Covers the full 8-step pipeline an applied economist or quantitative social scientist runs on every paper — (1) data cleaning, (2) variable construction & transformation, (3) descriptive statistics & Table 1, (4) statistical diagnostic tests, (5) baseline empirical modeling, (6) robustness battery, (7) further analysis (mechanism, heterogeneity, mediation, moderation), (8) publication-ready tables & figures. **Also covers two parallel domain modes that share the same 8-step scaf
Classical end-to-end empirical analysis workflow in the traditional Python econometric stack — pandas + numpy + scipy + statsmodels + linearmodels + pyfixest + rdrobust + econml + causalml + matplotlib/seaborn. **Defaults to economics empirical-paper style** (AER / QJE / AEJ) — every run produces a publication-ready output set with a multi-column regression table (M1→M6 progressive controls/FE) as the centerpiece, plus Table 1 (descriptives), mechanism / heterogeneity / robustness tables, and event-study + coefficient + trend figures. Covers the full 8-step pipeline an applied economist or quantitative social scientist runs on every paper — (1) data cleaning, (2) variable construction & transformation, (3) descriptive statistics & Table 1, (4) statistical diagnostic tests, (5) baseline empirical modeling, (6) robustness battery, (7) further analysis (mechanism, heterogeneity, mediation, moderation), (8) publication-ready tables & figures. **Also covers two parallel domain modes that share the same 8-step scaffolding** — **Mode A — Epidemiology / public health** (target-trial emulation via `zepid` / hand-rolled `pandas`, IPTW + g-formula + TMLE doubly-robust triplet via `zepid` / `econml` / `lifelines`, Mendelian randomization via `pymr` / `mrtool` (or `rpy2` → `MendelianRandomization`/`TwoSampleMR`), KM / AFT / Cox survival via `lifelines`, E-value sensitivity, principal stratification — STROBE / TRIPOD reporting), and **Mode B — ML causal inference** (DML via `econml.dml` / `doubleml`, S/T/X/R/DR meta-learners via `econml.metalearners` / `causalml`, causal forest via `econml.grf` / `causalml`, Dragonnet / TARNet / CEVAE neural causal via `causalml`, BCF via `pymc-bart` / `bcf-py`, matrix completion, CATE distribution + policy tree via `econml.policy` / `policytree-py`, off-policy evaluation, conformal causal via `mapie`, fairness audit via `fairlearn`, DAG learning via `causal-learn` / `cdt` / LLM-assisted). Prescribes which library to reach for at each step, shows the canonical code, and links to deeper `references/` files for variant-specific patterns. Use when the user asks for a **complete empirical analysis** in Python, wants to replicate an applied-economics paper from scratch, needs a reproducible workflow that is NOT opinionated on any single vertical package (contrast with StatsPAI), wants explicit control over every estimator and diagnostic, or asks "how do I write a full empirical pipeline in Python?". Also triggers when the user names a specific classical step in isolation — "winsorize at 1/99%", "run Breusch-Pagan", "build a Table 1 balance table", "do a placebo test", "event study plot", "mediation analysis" — and wants it wired into the broader pipeline. Mode A triggers on "target trial emulation", "IPTW", "TMLE", "Mendelian randomization", "STROBE", "公共健康", "流行病学". Mode B triggers on "DML", "double machine learning", "causal forest", "meta-learner", "Dragonnet", "BCF", "policy tree", "conformal causal", "fairness audit", "因果机器学习".
triggers
["full empirical analysis in Python","classical econometrics pipeline","traditional Python econometrics","end-to-end empirical workflow","pandas statsmodels linearmodels workflow","replicate an applied economics paper","data cleaning empirical","winsorize and standardize","variable construction","Table 1 summary statistics","balance table","correlation matrix","normality test","heteroskedasticity test","autocorrelation test","stationarity test","multicollinearity VIF","endogeneity test","baseline regression","panel fixed effects","DID workflow Python","event study","instrumental variables regression","regression discontinuity","propensity score matching","synthetic control python","double machine learning","robustness checks","placebo test","specification curve","alternative clustering","heterogeneity analysis","mechanism analysis","mediation analysis","moderation analysis","publication-ready regression table","coefplot","binscatter","event study plot","epidemiology pipeline python","public health causal inference python","target trial emulation python","g-formula python","IPTW marginal structural model python","TMLE doubly robust python","HAL-TMLE python","Mendelian randomization python","MR-Egger weighted median python","STROBE TRIPOD reporting python","E-value sensitivity python","Kaplan-Meier AFT survival python","lifelines survival python","zepid epidemiology","流行病学 python","公共健康 python","ML causal inference python","double machine learning DML python","econml DoubleML","meta-learner S T X R DR python","causal forest GRF python","causalml meta-learner","Dragonnet TARNet CEVAE python","Bayesian causal forest BCF python","CATE distribution python","policy tree python","off-policy evaluation python","conformal causal prediction python","mapie conformal","fairness audit python","fairlearn audit","causal discovery PC NOTEARS python","causal-learn cdt","因果机器学习 python"]
Full Empirical Analysis — Classical Python Workflow
This skill is the canonical 8-step pipeline an applied economist runs on every empirical paper, written in the traditional Python ecosystem — no opinionated one-stop wrapper. Every step calls libraries directly (pandas, numpy, scipy, statsmodels, linearmodels, pyfixest, rdrobust, econml, causalml, matplotlib, seaborn), so the agent — or the user reading the agent's code — has full visibility and can swap any component.
Companion skill: if the user prefers a single-import agent-native DSL (import statspai as sp), route to 00-StatsPAI_skill instead. This skill is the opposite philosophy: everything explicit, everything inspectable, every diagnostic run by hand, every plot shaped by the user.
Philosophy
Traditional stack, no magic. Agents should be able to read every line and know exactly which library / estimator / standard error family is at work.
Full pipeline, not just estimation. 80% of the time on a real paper is steps 1–4 and 6–8. This skill treats them as first-class, not an afterthought.
Rich outputs. Every step produces at least one table or figure — never a single point estimate in isolation.
Progressive disclosure. SKILL.md gives the canonical call at each step; references/ holds variant-specific depth (dozens of tests, estimator-specific diagnostics, plot recipes).
Reproducible. Every code block is runnable after pip install -r requirements.txt and df = pd.read_csv(...).
The default playbook above is AER-style applied econometrics — the AEA convention: written-out estimating equation, identifying assumption, 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 Steps 1–4 (cleaning / construction / descriptives / diagnostics) and Step 8 (tables & figures) — only Step 5 (estimator) and Step 6/7 (robustness / mechanism) swap libraries:
이 SKILL.md는 매우 커서 SkillsMP가 여기에는 첫 섹션만 미리 보여줍니다.GitHub에서 보기
"Mix" (e.g. "estimate DID + then ML CATE on the heterogeneity")
Default + Mode B in sequence — every estimator returns a coefficient + SE pair, drop them all into one pf.etable(...) for the horse-race column
The three modes share the same Step 1–4 cleaning / Table 1 / diagnostics scaffolding, the same Step 8 export stack, and the same DAG-first identification logic — switching modes only changes which Step-5 estimator family you reach for, not the surrounding paper structure. If you only want descriptive stats / Table 1 / a balance check, the AER tableone / gtsummary-style calls in Step 3 work identically across all three modes.
Default Output Spec — Economics Empirical Paper
This skill defaults to the applied-economics paper convention. Unless the user explicitly asks for a single point estimate, every run produces the full publication-ready output set below. Treat it as the contract of Step 8 — mandatory, not opt-in.
Required tables (always produced)
#
Table
Source / library
Saves to
T1
Summary statistics & balance (treated vs control, with SMD / p-values)
pandas.describe + custom table1() (Step 3)
tables/table1_balance.xlsx + .docx + .tex
T2 ★
Main results — multi-column regression M1→M6 (progressive controls + FE)
Mechanism / outcome ladder — same treatment, 3+ outcomes side-by-side
feols looped over y ∈ {Y1, Y2, Y3, Y_main} → pf.etable
tables/table3_mechanism.xlsx + .docx + .tex
T4
Heterogeneity — subgroup × main coef (gender, age, region, …)
subgroup feols × Wald → pf.etable (Step 7)
tables/table4_heterogeneity.xlsx + .docx + .tex
T5
Robustness battery — alt SE / alt cluster / alt sample / placebo, in one table
feols × variants → pf.etable (Step 6)
tables/table5_robustness.xlsx + .docx + .tex
★ Table 2 is the centerpiece of every economics paper. It is the multi-column regression table that walks the reader from raw correlation (M1) to the fully-specified design (M6: 2-way FE + interacted FE + cluster-robust SE). Do not collapse it into a single column. Do not report only the headline coefficient. The progression is the credibility argument: if M1→M6 is monotone and stable, the design is plausibly identifying; if it collapses on adding FE, that is the result.
The 8 steps mirror the canonical sections of an applied AER / QJE / AEJ paper. Each step is one paper section and emits a paper-ready artifact on disk:
Below is the canonical code at each step. All examples share one running narrative — a labor-economics panel where training (treatment) affects log_wage (outcome), with covariates age, edu, tenure, panel keys worker_id / firm_id / year. Column names and parameter values are illustrative — substitute the real ones from the user's DataFrame. Only library names and call shapes are normative.
When a step has many variants (e.g. staggered DID has five different estimators; heteroskedasticity has four classic tests), SKILL.md shows the one you reach for first and links to references/NN-<topic>.md for the rest. Read the reference file when the user's case doesn't fit the default.
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 leaves 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, the build system can use LaTeX, and editors can edit raw numbers in Excel). 所有图片必须同时保存 PNG（≥300 dpi）和 PDF 两种格式。
§
Artifact
Python primitive
Filenames
§1
Figure 1: raw trends / treatment rollout
df.groupby([time,treat])[y].mean().unstack().plot() · seaborn.heatmap for staggered rollout
Every Python estimator above (pf.feols / IV2SLS / att_gt via R-callout / CausalForestDML) returns a result object that can be passed straight into pf.etable(...) / pf.coefplot(...) / Stargazer(...). Don't hand-roll LaTeX from df.to_latex(), and don't render Word via python-docx directly — pf.etable / Stargazer apply book-tab borders, AER stars, and the right SE label automatically. For deeper export recipes (LaTeX / Word / Markdown variants, multi-panel .docx, full gtsummary-style flow), see references/08-tables-plots.md.
Export cookbook — LaTeX / Word / Excel in one block
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
Stack via repeat pf.etable(..., extralines=...) calls, or use gtsummary-style chained tables; for true single-file multi-panel, write to a .tex then concat
first panel: write; subsequent: append; surround with LaTeX \section{} headers
3. Full session bundle (the Stata collect / R gt equivalent)
Replication appendix that mixes summary stats + balance + multiple regression tables + headings + prose in one file
Compose programmatically with pylatex / python-docx / quarto — render once. Or use statsmodels.iolib.summary2.summary_col for a quick concat of tables.
Journal styling — pick the right signif_code and SE label. AEA convention is [0.1, 0.05, 0.01] and SE label "Cluster-robust standard errors in parentheses". Define a wrapper once at the top of master.py:
# top of master.py — journal house-style wrapper
AER_SIGNIF = [0.1, 0.05, 0.01]
AER_NOTES = ("Cluster-robust standard errors in parentheses. ""* p<0.10, ** p<0.05, *** p<0.01.")
defaer_table(models, *, file, headers=None, coef_map=None):
# 同时导出三种格式：.xlsx（用于编辑）、.docx（用于Word）、.tex（用于LaTeX）
base, ext = os.path.splitext(file)
for ext, type_ in [(".xlsx", "xlsx"), (".docx", "docx"), (".tex", "tex")]:
pf.etable(models, type=type_, file=base + ext,
headers=headers, coef_map=coef_map,
digits=3, signif_code=AER_SIGNIF, notes=AER_NOTES)
For the multi-panel .docx / .xlsx and Markdown / Quarto cookbook (single-file paper-tables bundle), see references/08-tables-plots.md.
Step −1 — Pre-Analysis Plan (pre-data; AEA RCT Registry style)
Before touching the data, write down (a) the population, (b) the design, (c) the minimum detectable effect (MDE) under the planned sample size and α=0.05, β=0.20. Persist the result as pap.json so a referee can verify the design was powered before, not after, the data were seen.
import json
from statsmodels.stats.power import TTestIndPower, NormalIndPower
from statsmodels.stats.proportion import samplesize_proportions_2indep_onetail
# Two-sample MDE for a continuous outcome (Cohen's d framing)
analysis = TTestIndPower()
n_required = analysis.solve_power(effect_size=0.20, power=0.80, alpha=0.05, ratio=1.0)
print(f"n per arm for d=0.20, 80% power: {n_required:.0f}")
# Solve for MDE given fixed n
mde = analysis.solve_power(nobs1=2000, power=0.80, alpha=0.05, ratio=1.0)
print(f"MDE (Cohen's d) at n=2000 per arm: {mde:.3f}")
# Cluster-randomized RCT — design effect = 1 + (m-1)·ICC
m, icc = 50, 0.05
deff = 1 + (m - 1) * icc
n_eff_required = analysis.solve_power(effect_size=0.20, power=0.80, alpha=0.05) * deff
print(f"n per arm under ICC={icc}, cluster size={m}: {n_eff_required:.0f}")
# DID: use Frison-Pocock / Bloom (1995) — see references/05-modeling.md §5.4# RD: power via Monte Carlo — see references/05-modeling.md §5.5# Persist the protocol — the referee will ask whether the design was powered ex ante
pap = {
"population": "manufacturing workers, 2010–2020",
"treatment": "training (binary, staggered adoption)",
"outcome": "log_wage",
"estimand": "ATT",
"design": "staggered DID, Callaway–Sant'Anna",
"alpha": 0.05,
"power_target": 0.80,
"mde_d": 0.20,
"n_planned": 12000,
"frozen_at": "2026-01-15T09:00:00Z",
"git_sha": "<paste>",
}
withopen("artifacts/pap.json", "w") as f:
json.dump(pap, f, indent=2)
Commit artifacts/pap.json in the repo before Step 1. AEA RCT Registry / OSF preregistration tools accept it as the analysis-plan exhibit.
Step 0 — Sample-construction log & 5-check data contract
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. The data contract is the cure for "mysterious sample-size shrinkage" bugs in the response letter.
0.1 Sample-construction log (footnote 4)
import pandas as pd, json
sample_log = []
df_raw = pd.read_csv("raw.csv")
sample_log.append(("0. raw", len(df_raw)))
df1 = df_raw.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","construction","transport"})]
sample_log.append(("3. keep target industries", len(df3)))
df = df3
for label, n in sample_log:
print(f" {label:<30s} N = {n:>10,d}")
withopen("artifacts/sample_construction.json", "w") as f:
json.dump(sample_log, f, indent=2)
Paste the printed lines verbatim as footnote 4 of the paper.
0.2 Five-check data contract (go / no-go gate)
import pandas as pd, numpy as np, json
from scipy import stats
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"n_missing": df[keys].isna().sum().to_dict(), # 3. missingness"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. dups
bal = df.groupby(id).size()
c["panel_balanced"] = bool((bal == bal.max()).all()) # 5. balance
c["n_dropped_by_balance"] = int((bal != bal.max()).sum())
if"first_treat_year"in df.columns:
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())
# MCAR sniff test (Rubin) — if missing(y) is associated with covariates,# listwise deletion biases the estimate. Use MI / IPW instead.
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}) "f"→ use MI / IPW, NOT listwise drop")
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, f"dup (id, time): {contract['n_dupes_on_keys']}"assertall(v == 0for v in contract["n_missing"].values()), \
f"NaNs on keys: {contract['n_missing']}"withopen("artifacts/data_contract.json", "w") as f:
json.dump(contract, f, indent=2, default=str)
If any assertion fires, stop and fix it in pandas. Estimators silently drop NaN rows downstream — this contract is the cheapest insurance against "why did N drop from 12,000 to 9,800 between Table 1 and Table 2?" referee questions.
import pandas as pd
import numpy as np
df = pd.read_csv("raw.csv")
# 1a. Inspect — always do this first
df.info() # dtypes + non-null counts
df.describe(include="all").T # numeric + categorical
df.isna().mean().sort_values(ascending=False) # missingness share per column# 1b. Fix dtypes (strings-that-should-be-numeric are the #1 silent bug)
df["year"] = pd.to_numeric(df["year"], errors="coerce")
df["wage"] = pd.to_numeric(df["wage"], errors="coerce")
df["gender"] = df["gender"].astype("category")
df["date"] = pd.to_datetime(df["date"], errors="coerce")
# 1c. Missing values — decide PER VARIABLE, never blanket-drop
key_vars = ["wage", "training", "worker_id", "year"]
df = df.dropna(subset=key_vars) # drop rows missing on keys
df["tenure"] = df["tenure"].fillna(df["tenure"].median()) # median-impute numeric covariate
df["union"] = df["union"].fillna("unknown") # explicit "unknown" for categorical# 1d. Outliers — flag first, winsorize in Step 2
df["wage_z"] = (df["wage"] - df["wage"].mean()) / df["wage"].std()
outlier_mask = df["wage_z"].abs() > 4print(f"{outlier_mask.sum()} rows flagged as |z|>4 on wage")
# 1e. Deduplicate on the panel key
dupes = df.duplicated(subset=["worker_id", "year"], keep=False)
assert dupes.sum() == 0, f"{dupes.sum()} duplicate (worker_id, year) pairs"# 1f. Merge auxiliary data — use validate= to catch silent m:m blowups
df = df.merge(firm_chars, on="firm_id", how="left", validate="many_to_one")
# 1g. Panel structure check — balanced vs. unbalanced
panel_summary = df.groupby("worker_id")["year"].agg(["count", "min", "max"])
print(panel_summary.describe())
is_balanced = (panel_summary["count"] == panel_summary["count"].max()).all()
print(f"Balanced: {is_balanced}")
Key principle: pandas + explicit decisions. Never silently drop rows inside an estimator — all row exclusions happen in Step 1 with a printed count.
Step 2 — Variable construction & transformation
Deeper patterns: references/02-data-transformation.md — log/ihs/Box–Cox, winsorizing vs. trimming, within-group standardization, one-hot vs. ordinal vs. target encoding, interaction terms, lag/lead operators, first differences, deflation with CPI.
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. Persist the strategy as strategy.md so it is a dated, version-controlled artifact — not a post-hoc rationalization written after seeing the coefficient.
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.5.1 Event-study figure + 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.
import pyfixest as pf
import matplotlib.pyplot as plt
# (a) Sun-Abraham via pyfixest::sunab — the modern primary for staggered DID
es = pf.feols("log_wage ~ sunab(first_treat_year, year) | worker_id + year",
data=df, vcov={"CRV1": "worker_id"})
# (b) Coefficient figure
fig = pf.iplot(es,
figsize=(7, 4),
title="Figure 2a. Event-study coefficients (95% CI; ref. e = -1)",
xlabel="Years relative to treatment",
ylabel="Coefficient (ATT)")
fig.savefig("figures/fig2a_event_study.pdf", dpi=300)
fig.savefig("figures/fig2a_event_study.png", dpi=300)
# (c) Numerical pre-trends F-test (Wald on the leads jointly = 0)import numpy as np
pre_idx = [k for k in es.coef().index if"rel_time::-"in k and"ref"notin k]
W = es.wald_test(pre_idx)
print(f"Pre-trends Wald χ² = {W.statistic:.2f}, p = {W.pvalue:.3f}")
# (d) Bacon decomposition (Goodman-Bacon 2021) — TWFE diagnostic# Pure-Python: callout to R::bacondecomp via rpy2 (Stata/R have first-class support)# OR use the `bacondecomp` Python port if installed.try:
from bacondecomp import bacon
bd = bacon(df, y="log_wage", D="training", time="year", id_var="worker_id")
bd.plot(); plt.savefig("figures/fig2a_bacon.pdf", dpi=300)
except ImportError:
print("bacondecomp not installed; use R callout or Stata's -bacondecomp-.")
Rule of thumb: first-stage F ≥ 10 for OLS-style inference; F ≥ 23 for AR-equivalent inference (Stock–Yogo / Lee 2022).
from linearmodels.iv import IV2SLS
iv = IV2SLS.from_formula(
"log_wage ~ 1 + age + edu + [training ~ Z1 + Z2]",
data=df).fit(cov_type="clustered", clusters=df["firm_id"])
print(iv.first_stage) # reports first-stage Fprint(iv.summary)
# Binscatter for the first-stage scatter (residualized on age + edu)from binsreg import binsreg # pip install binsreg
res = binsreg(y=df["training"], x=df["Z1"], w=df[["age","edu"]],
nbins=20, polyreg=2, ci=(3,3))
res.bins_plot.savefig("figures/fig2b_first_stage.pdf", dpi=300)
3.5.3 RD: McCrary density + canonical RD plot
The signature RD figure is rdplot (CCT-style binned scatter with local-polynomial fit on each side), paired with the McCrary manipulation test.
from rdrobust import rdplot, rdrobust
from rddensity import rddensity, rdplotdensity
# (a) Canonical RD plot — binned means + local poly on each side
rdp = rdplot(y=df["outcome"], x=df["running_var"], c=0,
p=4, kernel="triangular", binselect="esmv")
plt.savefig("figures/fig2c_rdplot.pdf", dpi=300)
# (b) McCrary density (Cattaneo-Jansson-Ma 2018)
dens = rddensity(X=df["running_var"], c=0)
print(dens)
rdplotdensity(dens, X=df["running_var"])
plt.savefig("figures/fig2c_mccrary.pdf", dpi=300)
3.5.4 Matching: love plot (standardized differences pre vs post)
import causalml.matching as cm
from causalml.matchimport NearestNeighborMatch
import seaborn as sns
# Pre-matching SMDs
pre_smd = ((df.loc[df.training==1, ["age","edu","tenure"]].mean()
- df.loc[df.training==0, ["age","edu","tenure"]].mean())
/ df[["age","edu","tenure"]].std())
# Match
psm = NearestNeighborMatch(replace=False, ratio=1, random_state=42)
matched = psm.match(data=df, treatment_col="training",
score_cols=["age","edu","tenure"])
# Post-matching SMDs
post_smd = ((matched.loc[matched.training==1, ["age","edu","tenure"]].mean()
- matched.loc[matched.training==0, ["age","edu","tenure"]].mean())
/ matched[["age","edu","tenure"]].std())
love = pd.DataFrame({"pre": pre_smd.abs(), "post": post_smd.abs()})
love.plot.barh(); plt.axvline(0.10, ls="--", c="r")
plt.title("Figure 2d. Love plot — |SMD| pre vs post matching (target < 0.10)")
plt.savefig("figures/fig2d_loveplot.pdf", dpi=300)
3.5.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.
# Pure Python: pysynth or sparseSC; for SDiD use rpy2 callout to R::synthdidfrom SyntheticControlMethods import Synth
sc = Synth(df, "outcome", "unit_id", "time", treatment_period=2015,
treated_unit=1, n_optim=10)
sc.plot(["original", "pointwise"], treated_label="Treated unit",
synth_label="Synthetic control")
plt.savefig("figures/fig2e_synth_trajectory.pdf", dpi=300)
Identification-specific checks (PT for DID, weak-IV F, density for RD, common support for matching) are also auto-run inside the Step-5 estimators — don't duplicate the numerics here, but DO produce the figures: a referee scans the figures first.
Step 4 — Diagnostic statistical tests
Deeper patterns: references/04-statistical-tests.md — every classical test with its null/alternative/decision rule (Shapiro–Wilk, Kolmogorov–Smirnov, Jarque–Bera, Breusch–Pagan, White, Goldfeld–Quandt, Durbin–Watson, Breusch–Godfrey, Ljung–Box, ADF, KPSS, Phillips–Perron, VIF, condition number, Hausman, Wu–Hausman, Sargan–Hansen).
Run diagnostics before taking estimates at face value. The 5 classes below cover 90% of applied work.
import statsmodels.api as sm
from statsmodels.stats.diagnostic import (
het_breuschpagan, het_white, acorr_breusch_godfrey, acorr_ljungbox,
)
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.stats.stattools import durbin_watson, jarque_bera
from scipy import stats as sps
# Fit a baseline OLS to get residuals for diagnostics
X = sm.add_constant(df[["training","age","edu","tenure"]])
y = df["log_wage"]
ols = sm.OLS(y, X, missing="drop").fit()
# 4a. Normality of residuals (informative but NOT required for OLS — CLT handles large N)
jb_stat, jb_p, skew, kurt = jarque_bera(ols.resid)
sw_stat, sw_p = sps.shapiro(ols.resid.sample(min(5000, len(ols.resid))))
print(f"Jarque-Bera p={jb_p:.3f} Shapiro p={sw_p:.3f} skew={skew:.2f} kurt={kurt:.2f}")
# 4b. Heteroskedasticity — Breusch-Pagan + White
bp = het_breuschpagan(ols.resid, ols.model.exog)
wh = het_white (ols.resid, ols.model.exog)
print(f"Breusch-Pagan p={bp[1]:.3f} White p={wh[1]:.3f}")
# → if p<0.05, use robust / cluster-robust SEs (you already should)# 4c. Autocorrelation (time-series / panel) — Durbin-Watson + Breusch-Godfrey + Ljung-Box
dw = durbin_watson(ols.resid) # ~2 = no AR(1)
bg = acorr_breusch_godfrey(ols, nlags=4) # general AR(p)
lb = acorr_ljungbox(ols.resid, lags=[4,8], return_df=True)
print(f"Durbin-Watson={dw:.2f} Breusch-Godfrey p={bg[1]:.3f}")
print(lb)
# 4d. Multicollinearity — VIF + condition number
vif = pd.DataFrame({
"var": X.columns,
"VIF": [variance_inflation_factor(X.values, i) for i inrange(X.shape[1])]
})
print(vif) # VIF > 10 is the classic red flag
cond_number = np.linalg.cond(X.values) # > 30 = potential collinearityprint(f"Condition number = {cond_number:.1f}")
# 4e. Stationarity (time-series) — ADF + KPSS (dual test: ADF rejects unit root, KPSS accepts stationarity)from statsmodels.tsa.stattools import adfuller, kpss
adf_stat, adf_p, *_ = adfuller(df["log_wage"].dropna(), autolag="AIC")
kpss_stat, kpss_p, *_ = kpss (df["log_wage"].dropna(), regression="c", nlags="auto")
print(f"ADF p={adf_p:.3f} KPSS p={kpss_p:.3f}")
Decision table (classic rules of thumb):
Test
Null
Action if rejected
Jarque–Bera / Shapiro
residuals ~ Normal
usually ignore when N large; bootstrap CIs if small-N inference matters
Breusch–Pagan / White
homoskedastic errors
use cov_type="HC3" or cluster SEs
Durbin–Watson / Breusch–Godfrey
no autocorrelation
use HAC (Newey–West) or cluster by unit
VIF > 10 / cond# > 30
—
drop / combine collinear regressors
ADF rejects + KPSS fails to reject
series is stationary
fit levels
ADF fails to reject
unit root
first-difference or cointegration test
Step 5 — Baseline empirical modeling (Section 4: Main Results)
Deeper patterns: references/05-modeling.md — every classical estimator with API: OLS, WLS, GLS, logit/probit, panel FE / RE / PO, clustered SEs, 2SLS / LIML / GMM, DID (2×2, TWFE, event study, CS, SA, BJS, SDiD), RD (sharp, fuzzy, kink, multi-cutoff), Synthetic Control, PSM / IPW / EB, DML / causal forest / DR-Learner.
This is the densest section of an applied paper. A modern AER §4 typically contains 2–3 multi-regression tables and one coefficient plot:
2SLS / IV → IV2SLS.from_formula("y ~ X + [D ~ Z]", df).fit(...) or pf.feols("y ~ X | D ~ Z", df)
DID / event-study → pf.feols("y ~ sunab(G, t) | i + t", df) for SA; R callout to did::att_gt for CS
Pick the estimator by identification strategy (not by "what's trendy"):
Observational cross-section, selection on observables → OLS + controls | PSM / IPW / DML
Observational panel, policy shock, parallel trends → DID (TWFE / CS / SA / BJS / SDiD)
Exogenous instrument for endogenous X → 2SLS / LIML / GMM (linearmodels / pyfixest)
Discontinuity in assignment rule → Sharp / Fuzzy / Kink RD (rdrobust)
N=1 treated unit, long panel → Synthetic Control (pysynth / SDiD)
High-dim controls or heterogeneous effects → DML / Causal Forest (econml)
Binary outcome → Logit / Probit (statsmodels)
Count outcome → Poisson / NegBin (pyfixest / statsmodels)
Canonical calls (details in references/05-modeling.md). The eight regression-table patterns A–H below are the AER table cookbook — pf.etable(*models, ...) is the workhorse, equivalent to Stata's outreg2/esttab and R's modelsummary.
5.A 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 6.f).
import pyfixest as pf
m1 = pf.feols("log_wage ~ training", data=df, vcov={"CRV1":"firm_id"})
m2 = pf.feols("log_wage ~ training + age + edu", data=df, vcov={"CRV1":"firm_id"})
m3 = pf.feols("log_wage ~ training + age + edu + tenure + firm_size", data=df, vcov={"CRV1":"firm_id"})
m4 = pf.feols("log_wage ~ training + age + edu + tenure + firm_size | industry + year",
data=df, vcov={"CRV1":"firm_id"})
m5 = pf.feols("log_wage ~ training + age + edu + tenure + firm_size | worker_id + year",
data=df, vcov={"CRV1":"firm_id"})
m6 = pf.feols("log_wage ~ training + age + edu + tenure + firm_size | worker_id + year + industry^year",
data=df, vcov={"CRV1":"firm_id"})
pf.etable([m1, m2, m3, m4, m5, m6],
type="tex", file="tables/table2_main.tex",
headers=["(1) Baseline", "(2) +Demog", "(3) +Labor-mkt",
"(4) Ind×Yr FE", "(5) Worker FE", "(6) Worker FE+Ind×Yr"],
digits=3, signif_code=[0.1, 0.05, 0.01],
notes="Cluster-robust SE in parentheses, clustered at firm_id.")
pf.etable([m1, m2, m3, m4, m5, m6], type="docx", file="tables/table2_main.docx")
AER convention: show ALL controls (and the intercept). Pass NEITHER keep= NOR drop= so every parameter is visible. Use keep=["training"] only when a focal-coefficient-only table is intentional (interaction-form heterogeneity, IV first-stage triplet); use drop=["Intercept"] only when you want to suppress the constant for paper aesthetics.
5.B 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.
from linearmodels.iv import IV2SLS
from econml.dml import LinearDML
from causalml.matchimport NearestNeighborMatch
ols = pf.feols("log_wage ~ training + age + edu + tenure | industry + year",
data=df, vcov={"CRV1":"firm_id"})
iv = IV2SLS.from_formula("log_wage ~ 1 + age + edu + tenure + [training ~ Z1 + Z2]",
data=df).fit(cov_type="clustered", clusters=df["firm_id"])
did = pf.feols("log_wage ~ sunab(first_treat_year, year) | worker_id + year",
data=df, vcov={"CRV1":"firm_id"})
dml = LinearDML().fit(df["log_wage"], df["training"],
X=df[["age","edu","tenure","firm_size"]])
psm = NearestNeighborMatch(replace=False, ratio=1).match(
df, treatment_col="training", score_cols=["age","edu","tenure"])
# Wrap non-pyfixest models or use Stargazer for a multi-source table:from stargazer.stargazer import Stargazer
table = Stargazer([ols.fit, iv, dml._final_estimator])
table.title("Table 2-bis. Convergent evidence across designs")
table.custom_columns(["(1) OLS+FE", "(2) 2SLS", "(3) DML"], [1,1,1])
open("tables/table2b_designs.tex", "w").write(table.render_latex())
5.C Pattern C — Multi-outcome table (same X, several Y's)
ys = ["log_wage", "weeks_employed", "left_firm", "promoted"]
multi_y = [pf.feols(f"{y} ~ training + age + edu + tenure | industry + year",
data=df, vcov={"CRV1":"firm_id"}) for y in ys]
pf.etable(multi_y, type="tex", file="tables/table2c_multi_outcome.tex",
headers=ys, keep="training",
notes="Each column is a separate regression on the labelled outcome.")
5.D Pattern D — Stacked Panel A / Panel B table
Same model family, two horizons (short-run / long-run) or two samples. Stack vertically with two pf.etable calls + LaTeX glue.
panelA = [pf.feols("wage_t1 ~ training + X | industry + year", data=df, vcov={"CRV1":"firm_id"}),
pf.feols("wage_t1 ~ training + X | worker_id + year", data=df, vcov={"CRV1":"firm_id"})]
panelB = [pf.feols("wage_t5 ~ training + X | industry + year", data=df, vcov={"CRV1":"firm_id"}),
pf.feols("wage_t5 ~ training + X | worker_id + year", data=df, vcov={"CRV1":"firm_id"})]
# Write Panel A
pf.etable(panelA, type="tex", file="tables/table2d_horizons.tex",
headers=["(1) Industry FE", "(2) Worker FE"], keep="training",
custom_row=[("Horizon", "1 year", "1 year")])
# Append Panel B (manually concat via texdoc-style glue, or use a wrapper)
5.E 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.
fs = pf.feols("training ~ Z + age + edu | industry + year", data=df, vcov={"CRV1":"firm_id"})
rf = pf.feols("log_wage ~ Z + age + edu | industry + year", data=df, vcov={"CRV1":"firm_id"})
iv2 = IV2SLS.from_formula("log_wage ~ 1 + age + edu + [training ~ Z]",
data=df).fit(cov_type="clustered", clusters=df["firm_id"])
pf.etable([fs, rf], type="tex", file="tables/table2e_iv_triplet.tex",
headers=["(1) First stage", "(2) Reduced form"], keep=["Z"],
notes=("First-stage F = " + f"{fs.fitstat()['F']:.2f}"))
# Append the 2SLS column (linearmodels IV2SLS output) via Stargazer or by hand.
IV triplet is intentionally focal: show only Z + endogenous regressor so the reader can eyeball the Wald ratio. Drop keep=["Z"] only if a referee asks for the full coefficient list.
5.F Pattern F — Causal-orchestrator main via did::att_gt / synth / econml.DML
For DID / SCM / matching mains, the modern Python estimator returns a self-contained estimate + automatic placebos / pre-trends / overlap diagnostics. Pipe into pf.etable via a thin adapter, or use Stargazer directly.
# DML with full diagnosticsfrom econml.dml import CausalForestDML
dml = CausalForestDML(n_estimators=2000, min_samples_leaf=5)
dml.fit(df["log_wage"], df["training"], X=df[["age","edu","tenure","firm_size"]])
ate = dml.ate(df[["age","edu","tenure","firm_size"]])
ci = dml.ate_interval(df[["age","edu","tenure","firm_size"]])
print(f"DML ATE = {ate:.3f} (95% CI [{ci[0]:.3f}, {ci[1]:.3f}])")
5.G Pattern G — Subgroup pf.etable (Table 3, see Step 7)
One column per subgroup. Detailed code in §Step 7 — Heterogeneity.
5.H Pattern H — Robustness master (Table A1, see Step 6)
Stack every robustness specification next to the baseline. Detailed code in §Step 6.f.
Canonical estimator commands (the underlying primitives)
import statsmodels.formula.api as smf
import pyfixest as pf
from linearmodels.iv import IV2SLS
from rdrobust import rdrobust
# 5a. OLS with cluster-robust SEs (use when no FE or just one FE)
ols = smf.ols("log_wage ~ training + age + edu + tenure", data=df).fit(
cov_type="cluster", cov_kwds={"groups": df["firm_id"]})
print(ols.summary())
# 5b. Panel FE (unit + time) — reach for pyfixest first; fastest & mirrors R's fixest
fe = pf.feols("log_wage ~ training + age + edu + tenure | worker_id + year",
data=df, vcov={"CRV1": "worker_id"})
fe.summary()
# 5c. 2×2 DID
did = smf.ols("log_wage ~ treated * post + age + edu", data=df).fit(
cov_type="cluster", cov_kwds={"groups": df["worker_id"]})
# 5d. Event study (dynamic DID, base period = −1)
es = pf.feols("log_wage ~ i(rel_time, ref=-1) | worker_id + year",
data=df, vcov={"CRV1": "worker_id"})
pf.iplot(es)
# 5e. Staggered DID — Callaway–Sant'Anna (via diff-diff) OR Sun–Abraham (via pyfixest interactions)# See references/05-modeling.md §5.4 for the full staggered DID playbook.# 5f. IV / 2SLS
iv = IV2SLS.from_formula(
"log_wage ~ 1 + age + edu + [training ~ draft_lottery]", data=df
).fit(cov_type="clustered", clusters=df["firm_id"])
print(iv.first_stage) # first-stage F > 10 (ideally > 104)print(iv.summary)
# 5g. Sharp RD
rd = rdrobust(y=df["outcome"], x=df["running_var"], c=0,
kernel="triangular", bwselect="mserd")
print(rd)
# 5h. Binary outcome
logit = smf.logit("employed ~ training + age + edu", data=df).fit()
print(logit.summary())
print(logit.get_margeff().summary()) # marginal effects — the interpretable quantity
# 7a. Heterogeneity via full interaction (cleanest: lets you test the interaction coefficient)
het = pf.feols("log_wage ~ training + training:female + age + edu | worker_id + year",
data=df, vcov={"CRV1": "worker_id"})
het.summary() # the interaction coefficient IS the heterogeneity test# 7b. Subgroup estimation with Wald test of equalityfrom scipy.stats import chi2
male_r = pf.feols("log_wage ~ training | worker_id+year", data=df[df.female==0])
female_r = pf.feols("log_wage ~ training | worker_id+year", data=df[df.female==1])
diff = male_r.coef()["training"] - female_r.coef()["training"]
se = np.sqrt(male_r.se()["training"]**2 + female_r.se()["training"]**2)
wald = (diff/se)**2print(f"Wald = {wald:.2f}, p = {1-chi2.cdf(wald,1):.3f}")
# 7c. Triple-difference (DDD) — heterogeneity by a THIRD dimension
ddd = pf.feols(
"log_wage ~ treated*post*high_exposure | worker_id + year",
data=df, vcov={"CRV1":"firm_id"})
# 7d. Mechanism — "outcome ladder" (same treatment, three sequential outcomes)for out in ["hours_worked", "productivity", "log_wage"]:
r = pf.feols(f"{out} ~ training | worker_id + year", data=df)
print(out, r.coef()["training"], r.se()["training"])
# 7e. Mediation — Baron-Kenny (ok for simple linear setting; use Imai for rigor)# Total: Y = a + c·T + ε# Step 1: M = a1 + b·T + ε (does T affect M?)# Step 2: Y = a2 + c'·T + d·M + ε (direct effect of T holding M fixed)# Mediated effect = b · d
b_coef = smf.ols("hours_worked ~ training + age+edu", data=df).fit().params["training"]
d_coef = smf.ols("log_wage ~ training + hours_worked + age+edu", data=df) \
.fit().params["hours_worked"]
print(f"Indirect effect via hours = {b_coef*d_coef:.3f}")
# 7f. Moderation — add interaction + marginal-effect plot; see references/07 for the full recipe.# 7g. Heterogeneous treatment effects via causal forest (high-dim moderators)from econml.dml import CausalForestDML
cf = CausalForestDML(n_estimators=1000, min_samples_leaf=5)
cf.fit(df["log_wage"], df["training"], X=df[["age","edu","tenure","firm_size"]])
tau = cf.effect(df[["age","edu","tenure","firm_size"]]) # per-unit CATE
cf.feature_importances_ # which X drives heterogeneity
Step 8 — Publication tables & figures
This step is mandatory — every analysis run produces all 5 required tables (T1–T5) and all 4 required figures (F1–F4) defined in the Default Output Spec at the top of this skill. Do not skip Step 8 because "the regression already ran". A coefficient without a table and a figure is not how applied economics communicates a result.
Deeper patterns: references/08-tables-plots.md — stargazer and pf.etable() for regression tables; coefficient plots with CIs; event-study plots (pre/post coefficients with reference line); binscatter; forest plots for subgroup analysis; RD plots; LaTeX / Word / Excel export.
# ============================================================# 8a. ★ TABLE 2 — Main results, multi-column regression M1→M6# (the centerpiece of every economics paper)# ============================================================from stargazer.stargazer import Stargazer
table = Stargazer([r.fit for r in [ols_m1, ols_m2, ols_m3, ols_m4, ols_m5, ols_m6]])
table.title("Effect of training on log wage")
table.custom_columns(["(1)","(2)","(3)","(4)","(5)","(6)"], [1]*6)
open("tables/table2_main.tex","w").write(table.render_latex())
# Or pyfixest's etable — handles FE indicators automatically (preferred):
pf.etable([m1, m2, m3, m4, m5, m6],
type="tex", file="tables/table2_main.tex",
headers=["(1) Raw","(2) +Demog","(3) +Tenure",
"(4) +Unit FE","(5) +2-way FE","(6) +Ind×Year FE"],
digits=3, signif_code=[0.1, 0.05, 0.01],
notes="Cluster-robust SE in parentheses, clustered at worker_id.")
pf.etable([m1, m2, m3, m4, m5, m6], type="docx", file="tables/table2_main.docx")
# ============================================================# 8b. TABLE 1 — Summary statistics & balance# ============================================================# Built in Step 3 as `t1` (DataFrame). Export both formats:
t1.to_latex("tables/table1_balance.tex", float_format="%.3f", index=False)
t1.to_excel("tables/table1_balance.xlsx", index=False)
# .docx version via python-docx or pandas → docx through tabulate# ============================================================# 8c. TABLE 3 — Mechanism / outcome ladder (3+ outcomes)# ============================================================
ladder = [pf.feols(f"{y} ~ training + age + edu + tenure | worker_id + year",
data=df, vcov={"CRV1":"worker_id"})
for y in ["hours_worked", "productivity", "log_wage"]]
pf.etable(ladder, type="tex", file="tables/table3_mechanism.tex",
headers=["Hours worked", "Productivity", "Log wage"],
notes="Each column is a separate regression on the labelled outcome.")
# ============================================================# 8d. TABLE 4 — Heterogeneity (subgroup × main coef)# ============================================================
het_specs = {
"All": df,
"Female=0": df[df.female==0],
"Female=1": df[df.female==1],
"Age<40": df[df.age<40],
"Age>=40": df[df.age>=40],
"Manufacturing": df[df.industry.eq("manufacturing")],
}
het_models = [pf.feols("log_wage ~ training + age + edu + tenure | worker_id + year",
data=d, vcov={"CRV1":"worker_id"})
for d in het_specs.values()]
pf.etable(het_models, type="tex", file="tables/table4_heterogeneity.tex",
headers=list(het_specs.keys()),
notes="Cluster-robust SE at worker_id. Wald p-values for cross-subgroup equality ""should accompany this table — see references/07.")
# ============================================================# 8e. TABLE 5 — Robustness battery (alt SE / cluster / sample / placebo)# ============================================================
rob = {
"Baseline": pf.feols("log_wage ~ training | worker_id + year", df,
vcov={"CRV1":"worker_id"}),
"Cluster=firm": pf.feols("log_wage ~ training | worker_id + year", df,
vcov={"CRV1":"firm_id"}),
"Two-way clust": pf.feols("log_wage ~ training | worker_id + year", df,
vcov={"CRV3x1":["worker_id","firm_id"]}),
"Winsor 1/99": pf.feols("log_wage ~ training | worker_id + year",
df.assign(log_wage=df.log_wage.clip(*df.log_wage.quantile([.01,.99]))),
vcov={"CRV1":"worker_id"}),
"Drop manuf": pf.feols("log_wage ~ training | worker_id + year",
df[df.industry!="manufacturing"], vcov={"CRV1":"worker_id"}),
"Placebo (-3)": pf.feols("log_wage ~ fake_post | worker_id + year",
df, vcov={"CRV1":"worker_id"}),
}
pf.etable(list(rob.values()), type="tex", file="tables/table5_robustness.tex",
headers=list(rob.keys()))
# ============================================================# 8f. ★ FIGURE 3 — Coefficient plot across M1→M6# ============================================================
fig, ax = plt.subplots(figsize=(6, 3.5))
labels, betas, lows, highs = [], [], [], []
for name, r in [("(1)",m1),("(2)",m2),("(3)",m3),("(4)",m4),("(5)",m5),("(6)",m6)]:
b = r.coef()["training"]; se = r.se()["training"]
labels.append(name); betas.append(b)
lows.append(b-1.96*se); highs.append(b+1.96*se)
betas = np.array(betas); lows = np.array(lows); highs = np.array(highs)
ax.errorbar(labels, betas, yerr=[betas-lows, highs-betas], fmt="o", capsize=3, color="navy")
ax.axhline(0, ls="--", color="gray", alpha=.6)
ax.set_ylabel("ATT on log wage"); ax.set_xlabel("Specification")
plt.tight_layout()
plt.savefig("figures/fig3_coefplot.pdf"); plt.savefig("figures/fig3_coefplot.png", dpi=300)
# ============================================================# 8g. FIGURE 2 — Event-study plot (dynamic DID, base period = -1)# ============================================================
fig, ax = plt.subplots(figsize=(7, 4))
pf.iplot(es, ax=ax)
ax.axhline(0, ls="--", color="gray"); ax.axvline(-0.5, ls=":", color="gray")
ax.set_xlabel("Years relative to treatment"); ax.set_ylabel("Coefficient (ATT)")
plt.tight_layout()
plt.savefig("figures/fig2_event_study.pdf"); plt.savefig("figures/fig2_event_study.png", dpi=300)
# ============================================================# 8h. FIGURE 4 — Sensitivity / robustness curve (spec curve)# ============================================================# Loop over 32 spec combinations and rank by point estimate
specs_curve = []
for fe in ["", "| worker_id", "| worker_id + year", "| worker_id + year + industry^year"]:
for ctrl in [[], ["age"], ["age","edu"], ["age","edu","tenure"]]:
f = "log_wage ~ training" + ("+" + "+".join(ctrl) if ctrl else"") + " " + fe
try:
r = pf.feols(f, df, vcov={"CRV1":"worker_id"})
specs_curve.append({"spec": f.strip(), "b": r.coef()["training"],
"se": r.se()["training"]})
except Exception:
pass
sc = pd.DataFrame(specs_curve).sort_values("b").reset_index(drop=True)
fig, ax = plt.subplots(figsize=(7, 4))
ax.errorbar(range(len(sc)), sc["b"], yerr=1.96*sc["se"], fmt="o", ms=3, color="navy", alpha=.7)
ax.axhline(0, ls="--", color="gray")
ax.set_xlabel("Specification rank (sorted by point estimate)")
ax.set_ylabel("Coefficient on training")
plt.tight_layout()
plt.savefig("figures/fig4_sensitivity.pdf"); plt.savefig("figures/fig4_sensitivity.png", dpi=300)
# ============================================================# 8i. FIGURE 1 — Trend / motivation (treated vs control over time)# ============================================================# Already built in Step 3; re-export with paper-grade styling.
fig, ax = plt.subplots(figsize=(7, 4))
trend = df.groupby(["year","training"])["log_wage"].mean().unstack()
trend.plot(ax=ax, marker="o", color={0:"darkred", 1:"navy"})
ax.axvline(policy_year, ls="--", color="gray", label="Policy")
ax.set_ylabel("Mean log wage"); ax.set_xlabel("Year")
ax.legend(["Control","Treated","Policy"])
plt.tight_layout()
plt.savefig("figures/fig1_trend.pdf"); plt.savefig("figures/fig1_trend.png", dpi=300)
# ============================================================# 8j. Auxiliary plots (optional — produce when relevant)# ============================================================from binsreg import binsreg # binscatter
binsreg(y=df["log_wage"], x=df["tenure"], w=df[["age","edu"]], nbins=20)
plt.savefig("figures/figA_binscatter.pdf")
from rdrobust import rdplot # RD plot (only when running_var exists)# rdplot(y=df["outcome"], x=df["running_var"], c=0); plt.savefig("figures/figA_rdplot.pdf")# Forest plot for subgroups → see references/08-tables-plots.md §8.6 for the full recipe.
Deliverables checklist (verify before declaring the run complete):
The single artifact a journal's replication office (or a future co-author) needs to reproduce the headline number. Persist Python version, seed, dataset hash, baseline coefficient + CI, and pointers to the protocol/contract:
import json, sys, hashlib, pyfixest
# Get baseline result (assumes `base` is the headline pf.feols/sm.OLS object)
b_hat = float(base.coef()["training"])
se = float(base.se()["training"])
lo, hi = b_hat - 1.96*se, b_hat + 1.96*se
dataset_sha = hashlib.sha256(
pd.util.hash_pandas_object(df, index=True).values.tobytes()
).hexdigest()[:16]
stamp = {
"python_version": sys.version,
"pyfixest_version": pyfixest.__version__,
"seed": 42,
"dataset_sha256_16": dataset_sha,
"n_obs": int(base._N),
"estimand": "ATT",
"estimator": "pf.feols",
"estimate": b_hat,
"se_cluster": se,
"ci95": [lo, hi],
"pre_registration": "artifacts/strategy.md",
"data_contract": "artifacts/data_contract.json",
"sample_log": "artifacts/sample_construction.json",
"paper_bundle": "tables/table2_main.tex",
}
withopen("artifacts/result.json", "w") as f:
json.dump(stamp, f, indent=2)
Commit artifacts/result.json alongside the paper PDF. A referee should be able to run python master.py and bit-identically reproduce this JSON.
§A — Epidemiology / Public Health Mode
When the user's wording flags Mode A (target-trial emulation / IPTW / TMLE / MR / STROBE / 流行病学 / 公共健康 / RWE / cohort), the 8 steps still apply — but Step 5 swaps the OLS-and-FE stack for the doubly-robust + survival + MR triplet, and the deliverables follow STROBE / TRIPOD-AI conventions. Steps 1–4 (cleaning, construction, Table 1, diagnostics) and Step 8 (tables/figures export) are identical to the Default mode.
Library footprint (install on top of the Default stack):
pip install zepid # IPTW, g-formula, TMLE, AIPW, E-value
pip install lifelines # KM, Cox, AFT, RMST
pip install scikit-survival # alternative survival stack (scaling-friendly)# Mendelian randomization — Python coverage is thin; for IVW/Egger/weighted-median:
pip install pymr # if available# Or call R from Python:
pip install rpy2 # then import TwoSampleMR / MendelianRandomization via rpy2
A.0 Cohort construction + target-trial protocol
Write the protocol before touching the data. Save it as protocol.yml and quote it in the paper.
# protocol.yml — target-trial emulation skeleton
target_trial = {
"eligibility": {"age": "40-75", "no_prior_event": True, "ascertained_at": "t0"},
"treatment": {"A=1": "statin initiation", "A=0": "no initiation"},
"assignment": "random at t0 (emulated by IPTW on baseline covariates)",
"follow_up_start":"t0 (treatment initiation date)",
"outcome": "incident MI within 5 years",
"estimand": "intention-to-treat ATE on risk difference + hazard ratio",
"censoring": {"loss_to_FU": True, "competing_risk": "death from non-MI causes"},
}
# Cohort construction in pandas — eligibility + index date + censoring date
cohort = (df
.query("age >= 40 & age <= 75 & prior_MI == 0") # eligibility
.assign(t0 = lambda d: d["statin_initiation_date"].fillna(d["enrollment_date"]),
event_5y = lambda d: ((d["MI_date"] - d["t0"]).dt.days <= 365*5).astype(int),
time_at_risk = lambda d: ((d["censor_date"] - d["t0"]).dt.days.clip(0, 365*5))))
A.1 Table 1 by exposure (identical to Default Step 3)
Use the same tableone / gtsummary-style call from Step 3, just groupby="A" (treatment indicator). E-values for unmeasured confounding go in the footer.
A.2 DAG + propensity-score overlap (positivity check)
# Estimate PS, plot overlap (positivity), love-plot SMDsfrom sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
X = df[["age", "edu", "smoke", "bmi", "ldl", "sbp"]]
y = df["A"]
ps = LogisticRegression(max_iter=1000).fit(StandardScaler().fit_transform(X), y).predict_proba(StandardScaler().fit_transform(X))[:, 1]
df["ps"] = ps
# Overlap plotimport matplotlib.pyplot as plt, seaborn as sns
sns.kdeplot(data=df, x="ps", hue="A", common_norm=False)
plt.savefig("figures/figA2_ps_overlap.pdf")
# Love plot (SMDs before vs after IPTW) — see zepid.causal.ipw.diagnostics
