Full-empirical-analysis-skill

by brycewang-stanford4.6k

Communication

Install

The skill file stays in the publisher's repository.

Paste the copied command in a terminal. It installs this SKILL.md into that agent.

npx skills add brycewang-stanford/Auto-Empirical-Research-Skills --skill 00.1-Full-empirical-analysis-skill_Python -a cursor -y

Download skill skills/00.1-Full-empirical-analysis-skill_Python/SKILL.md on GitHub

Related skills

View all

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", "因果机器学习".

Questions

name: Full-empirical-analysis-skill description: 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:

import pandas as pd, numpy as np, json from scipy import stats 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 "n_missing": df[keys].isna().sum().to_dict(), # 3. missingness "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. 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") 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, f"dup (id, time): {contract['n_dupes_on_keys']}" assert all(v == 0 for v in contract["n_missing"].values()), \ f"NaNs on keys: {contract['n_missing']}" with open("artifacts/data_contract.json", "w") as f: json.dump(contract, f, indent=2, default=str)

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 in range(X.shape[1])] }) print(vif) # VIF > 10 is the classic red flag cond_number = np.linalg.cond(X.values) # > 30 = potential collinearity print(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}")

# 6a. Alternative specifications — progressive controls (M1 → M6) specs = [ "log_wage ~ training", # M1 raw "log_wage ~ training + age + edu", # M2 +covariates "log_wage ~ training + age + edu + tenure", # M3 +tenure "log_wage ~ training + age + edu + tenure | worker_id", # M4 +unit FE "log_wage ~ training + age + edu + tenure | worker_id + year", # M5 +time FE "log_wage ~ training + age + edu + tenure | worker_id + year + industry^year", # M6 +industry×year ] results = [pf.feols(f, data=df, vcov={"CRV1":"worker_id"}) for f in specs] pf.etable(results) # side-by-side publication table # 6b. Alternative cluster levels for cl in ["worker_id", "firm_id", "industry", "state"]: r = pf.feols("log_wage ~ training | worker_id+year", data=df, vcov={"CRV1": cl}) print(cl, r.coef()["training"], r.se()["training"]) # 6c. Subsample splits for col in ["male", "has_college", "high_tenure"]: for val in [0, 1]: sub = df[df[col]==val] r = pf.feols("log_wage ~ training | worker_id+year", data=sub) print(col, val, r.coef()["training"], r.se()["training"]) # 6d. Placebo — fake timing (treat 3 years before actual policy; should be ~0) df["fake_post"] = (df["year"] >= df["first_treat_year"] - 3).astype(int) placebo = pf.feols("log_wage ~ fake_post | worker_id + year", data=df) placebo.summary() # 6e. Placebo — permutation / randomization inference (500 draws) obs_coef = pf.feols("log_wage ~ training | worker_id+year", data=df).coef()["training"] draws = [] for s in range(500): df_s = df.copy() df_s["training_perm"] = df_s.groupby("worker_id")["training"].transform( lambda x: x.sample(frac=1, random_state=s).values) r = pf.feols("log_wage ~ training_perm | worker_id+year", data=df_s) draws.append(r.coef()["training_perm"]) p_perm = (np.abs(draws) >= abs(obs_coef)).mean() print(f"Permutation p = {p_perm:.3f}") # 6f. Oster (2019) δ* — selection on unobservables # See references/06-robustness.md §6.6 for the exact formula and its implementation. # ============================================================ # 6.g Pattern H — Robustness master table (Table A1, one column per check) # ============================================================ import pyfixest as pf base = pf.feols("log_wage ~ training + age + edu + tenure | industry + year", data=df, vcov={"CRV1":"firm_id"}) no99 = pf.feols("log_wage ~ training + age + edu + tenure | industry + year", data=df.query("wage < wage.quantile(0.99)"), vcov={"CRV1":"firm_id"}) balpan = pf.feols("log_wage ~ training + age + edu + tenure | industry + year", data=df.groupby("worker_id").filter(lambda g: len(g) == g["year"].nunique()), vcov={"CRV1":"firm_id"}) dropearly = pf.feols("log_wage ~ training + age + edu + tenure | industry + year", data=df.query("first_treat_year > 2008"), vcov={"CRV1":"firm_id"}) wfe = pf.feols("log_wage ~ training + age + edu + tenure | worker_id + year", data=df, vcov={"CRV1":"firm_id"}) cl2way = pf.feols("log_wage ~ training + age + edu + tenure | industry + year", data=df, vcov={"CRV1":"firm_id+year"}) logy = pf.feols("np.log1p(wage) ~ training + age + edu + tenure | industry + year", data=df, vcov={"CRV1":"firm_id"}) ihsy = pf.feols("np.arcsinh(wage) ~ training + age + edu + tenure | industry + year", data=df, vcov={"CRV1":"firm_id"}) pf.etable([base, no99, balpan, dropearly, wfe, cl2way, logy, ihsy], type="tex", file="tables/tableA1_robustness.tex", headers=["(1) Baseline", "(2) Drop top 1%", "(3) Balanced", "(4) Drop early", "(5) Worker FE", "(6) 2-way cluster", "(7) log Y", "(8) IHS Y"], keep="training", notes="Each column is one robustness check. β̂ on training is the focal coefficient.") # ============================================================ # 6.h Specification curve (Simonsohn-Simmons-Nelson 2020) # ============================================================ import itertools, numpy as np specs = [] for ctrls in [["age"], ["age","edu"], ["age","edu","tenure"], ["age","edu","tenure","firm_size"]]: for trans in ["log_wage", "ihs_wage"]: for mask_label, mask in [("all", slice(None)), ("manuf", df["industry"]=="manuf"), ("no99", df["wage"] < df["wage"].quantile(0.99))]: for cl in ["firm_id", "firm_id+year"]: d = df.loc[mask] if not isinstance(mask, slice) else df fit = pf.feols(f"{trans} ~ training + {' + '.join(ctrls)} | industry + year", data=d, vcov={"CRV1": cl}) b = fit.coef()["training"] se = fit.se()["training"] specs.append({ "label": f"c={'+'.join(ctrls)}/y={trans}/s={mask_label}/cl={cl}", "b": b, "lo": b - 1.96*se, "hi": b + 1.96*se}) specs_df = pd.DataFrame(specs).sort_values("b").reset_index(drop=True) fig, ax = plt.subplots(figsize=(10, 5)) ax.errorbar(specs_df.index, specs_df["b"], yerr=[specs_df["b"]-specs_df["lo"], specs_df["hi"]-specs_df["b"]], fmt="o", capsize=2, ms=3) ax.axhline(0, ls="--", c="gray") ax.set_xlabel("Specification (sorted by β̂)") ax.set_ylabel("Coefficient on training") ax.set_title("Figure 5. Specification curve") plt.tight_layout() plt.savefig("figures/fig5_spec_curve.pdf", dpi=300) plt.savefig("figures/fig5_spec_curve.png", dpi=300) # ============================================================ # 6.i Sensitivity dashboard — HonestDiD + Oster + E-value # ============================================================ # (a) HonestDiD — Rambachan-Roth (2023): Python via R callout (rpy2) or HonestDiD-py if installed. # (b) Oster δ — see references/06-robustness.md §6.6 for the formula implementation. # (c) E-value (VanderWeele-Ding 2017) — for risk-ratio outcomes; convert OLS to RR scale first. def evalue(rr_point, rr_lower): """Min strength of unmeasured confounding to nullify; (point, ci_lower).""" e_point = rr_point + np.sqrt(rr_point * (rr_point - 1)) e_ci = rr_lower + np.sqrt(rr_lower * (rr_lower - 1)) if rr_lower > 1 else 1.0 return e_point, e_ci print("E-value (point, CI):", evalue(1.45, 1.10))

# 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 equality from 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)**2 print(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

# ============================================================ # 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"])["l

Mode	Reader convention	Step-5 estimator 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 / `pyfixest.feols` HDFE	AER house-style multi-column `pf.etable` / `Stargazer` + 8-section paper layout	Steps 1 → 8 (entire playbook below)
Mode A — Epidemiology / Public Health	"STROBE / TRIPOD-AI; target trial protocol; doubly-robust estimand; absolute & relative risk; KM survival"	Target-trial emulation · IPTW (`zepid`) · g-formula (`zepid` / hand-rolled) · TMLE (`zepid.causal.gformula` / `econml`) · Mendelian randomization (`pymr` / `rpy2`+`TwoSampleMR`) · KM/AFT (`lifelines`)	Same `pf.etable` + 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 (`econml.dml` / `doubleml`) · S/T/X/R/DR-Learner (`econml.metalearners` / `causalml`) · GRF causal forest (`econml.grf`) · Dragonnet/TARNet/CEVAE (`causalml`) · BCF (`pymc-bart`/`bcf-py`) · matrix completion	`pf.etable` ML horse-race + CATE plot + policy-value table + conformal PI (`mapie`)	§B. ML causal pipeline

User says...	Mode the skill switches to
"Run a DID / IV / RD / event study", "AER table", "applied micro"	Default (AER econ) — Steps 1 → 8
"Target trial emulation", "g-formula", "IPTW", "TMLE", "Mendelian randomization", "STROBE / TRIPOD", "公共健康 / 流行病学", "epi pipeline", "RWE study", "cohort study", "case-control"	Mode A (Epi) — §A
"DML", "double machine learning", "causal forest", "meta-learner", "CATE", "Dragonnet", "BCF", "policy learning", "conformal causal", "fairness audit", "ML causal", "uplift modeling", "因果机器学习"	Mode B (ML causal) — §B
"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

#	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)	`pyfixest.feols` × 6 specs → `pf.etable()` / `Stargazer` (Step 5–6)	`tables/table2_main.xlsx` + `.docx` + `.tex`
T3	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`

#	Figure	Source / library	Saves to
F1	Trend / motivation — treated vs control over time, with policy line	`df.groupby([year, treat])[y].mean().unstack().plot()` (Step 3)	`figures/fig1_trend.png` (300 dpi, 必须导出 PNG) + `.pdf`
F2	Event-study coefficients with 95% CI, base period at –1	`pf.feols(... ~ i(rel_time, ref=-1) ...)` → `pf.iplot()` (Step 5)	`figures/fig2_event_study.png` (300 dpi, 必须导出 PNG) + `.pdf`
F3	Coefficient plot across specs M1→M6	`matplotlib.errorbar` over the 6 fitted models (Step 8)	`figures/fig3_coefplot.png` (300 dpi, 必须导出 PNG) + `.pdf`
F4	Robustness / sensitivity curve — spec curve, HonestDiD, or cluster-comparison	spec_curve loop or `honest_did` plot (Step 6)	`figures/fig4_sensitivity.png` (300 dpi, 必须导出 PNG) + `.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	`figures/fig1_trend.png`(300dpi)+`.pdf`
§1	Table 1: summary stats (full / treated / control + Δ + SMD)	`table1(df, by=, cols=)` (Step 3.b) → write LaTeX/Word/Excel	`tables/table1_balance.xlsx/.docx/.tex`
§3	Figure 2: identification graphic (event-study / first-stage / McCrary / RD scatter / SCM trajectory)	`pf.iplot(es)` · `binscatter` · `rdplot` · `rddensity` · `synthdid`	`figures/fig2_event_study.png`(300dpi)+`.pdf`
§4	Table 2: main results — progressive controls M1→M6	`pf.etable([m1...m6])` · `Stargazer([m1.fit ... m6.fit])`	`tables/table2_main.xlsx/.docx/.tex`
§4	Table 2-bis: design horse-race (OLS / IV / DID / DML)	`pf.etable([ols, iv, did, dml])`	`tables/table2b_designs.xlsx/.docx/.tex`
§4	Figure 3: coefficient plot across specs	`pf.coefplot([m1...m6], coefs=["training"])`	`figures/fig3_coefplot.png`(300dpi)+`.pdf`
§5	Table 3: heterogeneity by subgroup	`pf.etable(g_full, g_male, g_fem, g_q1...q4)`	`tables/table3_heterogeneity.xlsx/.docx/.tex`
§5	Figure 4: dose-response / CATE	`econml.CausalForestDML(...).effect()` + matplotlib hist	`figures/fig4_cate.png`(300dpi)+`.pdf`
§6	Table 4: mechanism / outcome ladder	loop `pf.feols` over outcomes → `pf.etable`	`tables/table4_mechanism.xlsx/.docx/.tex`
§7	Table A1: robustness master (one column per check)	`pf.etable([base, no99, balpan, dropearly, wfe, cl2way, logy, ihsy, psm, ebal])`	`tables/tableA1_robustness.xlsx/.docx/.tex`
§7	Figure 5: spec curve	hand-rolled `itertools.product` + matplotlib `errorbar`	`figures/fig5_spec_curve.png`(300dpi)+`.pdf`
§7	Figure 6: sensitivity (HonestDiD / Oster / E-value)	`HonestDiD` · `oster_bound` · `evalue`	`figures/fig6_sensitivity.png`(300dpi)+`.pdf`
§8	Replication bundle: all tables in one document	`pf.etable([...] + extra=[...], type="docx")` · pylatex / `texdoc`-style multi-panel	`replication/paper_tables.xlsx/.docx/.tex`

Tier	Use when	API	Hot kwargs
1. Single multi-column table	Exporting one Table 2 / Table 3 / Table A1 with progressive columns	`pf.etable([m1,...,mN], type="tex"/"docx"/"xlsx", file="...", headers=[...], digits=3, signif_code=[0.1,0.05,0.01])` — or `Stargazer([m1.fit,...,mN.fit]).render_latex()` for statsmodels-only	`keep=`, `drop=`, `coef_map=`, `headers=`, `digits=`, `signif_code=`, `fixef_rm=`, `notes=`
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-style template + per-section heading + footnote macros

Design	Estimating equation	Identifying assumption	Python estimator
2×2 DID	`Y_it = α_i + λ_t + β·D_it + X'γ + ε_it`	parallel trends conditional on X	`pf.feols("y ~ i(treated, post, ref=0)
Event-study (CS / SA)	`Y_it = α_i + λ_t + Σ_{e≠-1} β_e · 1{t-G_i = e} + ε_it`	no anticipation + group-time PT	`pf.feols("y ~ i(rel, ref=-1)
2SLS	`Y_i = α + β·D_i + X'γ + ε_i; D_i = π·Z_i + X'δ + u_i`	exclusion + relevance + monotonicity	`linearmodels.iv.IV2SLS.from_formula("y ~ X + [D ~ Z]", df)`
Sharp RD	`Y_i = α + β·1{X_i ≥ c} + f(X_i) + ε_i` (local poly)	continuity of E[Y(0)\|X] at c, no manipulation	`rdrobust(y, x, c=0)` (+ `rddensity`)
SCM	`Ŷ_1t(0) = Σ_j ŵ_j Y_jt`, τ_t = `Y_1t − Ŷ_1t(0)` for t≥T_0	pre-period fit + interpolation validity	`pysynth` · `synthdid` (R callout)
Selection-on-observables (matching/IPW/DML)	`Y_i = m(X_i) + β·D_i + ε_i` (Robinson partialling-out)	unconfoundedness + overlap	`econml.dml.DML` · `econml.dr.DRLearner` · `causalml.matching.PSM`

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

Install

Related skills

Questions

Mode A — Epidemiology / public health

Mode B — ML causal inference

Full Empirical Analysis — Classical Python Workflow

Philosophy

SkillOpt-style execution gate

Three domain modes (default = AER econ; alternates = epi & ML-causal)

Default Output Spec — Economics Empirical Paper

Required tables (always produced)

Required figures (always produced)

Output file layout (default)

When to deviate

Required Libraries

The 8 Steps — Canonical Pipeline (mapped to AER paper sections)

Paper-ready figure & table inventory (what to produce by section)

Export cookbook — LaTeX / Word / Excel in one block

Step −1 — Pre-Analysis Plan (pre-data; AEA RCT Registry style)

Step 0 — Sample-construction log & 5-check data contract

0.1 Sample-construction log (footnote 4)

0.2 Five-check data contract (go / no-go gate)

Step 1 — Data cleaning

Step 2 — Variable construction & transformation

Step 2.5 — Empirical strategy (write the equation + identifying assumption)

Equation × identifying assumption × Python estimator (decision table)

Design picker (when the user is unsure)

Pre-registration strategy.md template

Step 3 — Descriptive statistics & Table 1

Step 3.5 — Identification graphics (Section "Identification, graphical evidence")

3.5.1 Event-study figure + numerical pre-trends test (DID identification)

3.5.2 First-stage F-statistic + scatter (IV identification)

3.5.3 RD: McCrary density + canonical RD plot

3.5.4 Matching: love plot (standardized differences pre vs post)

3.5.5 SCM: synthetic-control trajectory + gap plot

Step 4 — Diagnostic statistical tests

Step 5 — Baseline empirical modeling (Section 4: Main Results)

5.A Pattern A — Progressive controls (the canonical Table 2)

5.B Pattern B — Design horse race (Table 2-bis)

5.C Pattern C — Multi-outcome table (same X, several Y's)

5.D Pattern D — Stacked Panel A / Panel B table

5.E Pattern E — IV reporting triplet (first-stage / reduced-form / 2SLS)

5.F Pattern F — Causal-orchestrator main via did::att_gt / synth / econml.DML

5.G Pattern G — Subgroup pf.etable (Table 3, see Step 7)

5.H Pattern H — Robustness master (Table A1, see Step 6)

Canonical estimator commands (the underlying primitives)

Step 6 — Robustness battery

Step 7 — Further analysis (mechanism / heterogeneity / mediation / moderation)

Step 8 — Publication tables & figures

Pre-registration `strategy.md` template

5.F Pattern F — Causal-orchestrator main via `did::att_gt` / `synth` / `econml.DML`

5.G Pattern G — Subgroup `pf.etable` (Table 3, see Step 7)