关于 Full-empirical-analysis-skill
面向 Python 的实证研究技能，提供应用经济学、流行病学和因果机器学习的八步工作流。它强调直接调用并检查各类库、执行诊断与稳健性分析，并产出论文用表格和图形。
主要功能
- 01八步实证分析工作流
- 02经济学、流行病学与因果机器学习模式
- 03显式使用传统 Python 库
- 04数据检查与诊断
- 05稳健性与敏感性分析
- 06生成论文用表格和图形
适用场景
- 01构建应用经济学论文分析
- 02运行流行病学分析流程
- 03比较因果机器学习估计器
- 04生成可复现研究成果
安装
复制下面的命令，在终端中运行。
$
npx skills add https:/安装后在 Claude Code 中输入/full-empirical-analysis-skill确认已加载。
包含的文件
除 SKILL.md 外，安装时一起下载的参考资料和脚本。
README.mdreferences/03-descriptive-stats.md02-data-transformation.md08-tables-plots.md05-modeling.md07-further-analysis.md06-robustness.md04-statistical-tests.md01-data-cleaning.md详细介绍
完整实证分析——传统 Python 工作流
此技能是应用经济学家处理每篇实证论文时采用的标准 8 步流程，使用传统 Python 生态系统编写——不依赖固执己见的一站式封装。每一步都直接调用库（pandas、numpy、scipy、statsmodels、linearmodels、pyfixest、rdrobust、econml、causalml、matplotlib、seaborn），因此 Agent——或阅读 Agent 代码的用户——都能全面了解具体实现，并可替换任意组件。
配套技能：如果用户更喜欢单次导入、原生支持 Agent 的 DSL（import statspai as sp），请改用 00-StatsPAI_skill。本技能的理念正好相反：一切都明确呈现，一切都可检查，每项诊断都手动运行，每张图都由用户定制。
理念
- 传统技术栈，不搞黑箱。 Agent 应当能够读懂每一行，并准确了解当前使用的库、估计器和标准误类型。
- 完整流程，而不只是估计。 在实际论文中，80% 的时间都花在步骤 1–4 和 6–8 上。本技能将这些步骤视为核心内容，而非事后补充。
- 丰富的输出。 每一步至少生成一张表格或一幅图——绝不孤立地只给出一个点估计。
- 渐进式披露。 SKILL.md 提供每一步的标准调用方式；
references/则包含针对不同变体的详细内容（数十种检验、针对特定估计器的诊断、绘图方案）。 - 可复现。 每个代码块都可以在运行
pip install -r requirements.txt和df = pd.read_csv(...)后执行。
SkillOpt 风格的执行关卡
将这份详细指南作为种子技能，而不是必须逐项执行的脚本。遵循 SkillOpt 准则：将每项局部分析代码更改视为候选补丁；它必须通过选择检查和留出检查，才能成为可复用的样板代码。在编写或修改分析代码之前，将用户的请求压缩成一张任务专属的 best_skill 卡片：
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>
- 路由卡片：记录模式（
econ、epi或ml-causal）、估计对象、识别设计、核心结果变量／处理变量、软件包栈以及所需产物。 - 限制改动范围：每次只更改一项决策（样本规则、估计器、聚类方式、导出格式或稳健性检验）。优先采用能够通过验证的最小补丁。
- 遵循选择样本拆分纪律：将用户当前遇到的问题或所需产物视为选择样本。至少预留一个替代结果变量、样本时间窗口、估计器系列或导出目标作为留出关卡。
- 留出关卡：运行代码前定义检查项：行数、键的唯一性、处理变量支持度、缺失值阈值、预期的表格／图片文件，以及一项不得被此次更改破坏的非核心稳健性检验／模型设定。
- 拒绝缓冲区：如果候选模型设定未通过关卡，就在
analysis_log.md中记录失败情况、代码差异和关卡输出；恢复到上一个通过检查的设定，并且不要重试同一未经检查的模式。 - 慢速／元更新：任务结束时，根据整个过程写下
accepted_patterns和rejected_patterns。不要仅凭一次通过的运行就扩展规范项目模板。 - 仅在验证后推广：只有在当前数据和至少一个替代结果变量／样本／模型设定上通过验证后，才能将一次性修复转为可复用的项目样板代码。
三种领域模式（默认 = AER 经济学；可选 = 流行病学与 ML 因果推断）
上文的默认流程是 AER 风格的应用计量经济学——遵循 AEA 惯例：写出估计方程、识别假设、设计对比检验，以及完整的稳健性检验。这项技能还附带了另外两种主流因果推断传统的并行子流程；每种流程都复用相同的步骤 1–4（清理／构建／描述性分析／诊断）和步骤 8（表格与图表）——只有步骤 5（估计器）和步骤 6/7（稳健性／机制）会切换工具库：
| 模式 | 读者惯例 | 步骤 5 估计器栈 | 报告栈 | 跳转到 |
|---|---|---|---|---|
| 默认 — 应用经济学（AER / QJE / AEJ） | “展示方程 + 识别假设 + 设计对比检验；控制变量清晰可见；聚类标准误” | DID / IV / RD / SCM / 匹配 / pyfixest.feols HDFE | AER 风格的多列 pf.etable / Stargazer + 8 节论文布局 | 步骤 1 → 8（下方完整流程） |
| 模式 A — 流行病学 / 公共卫生 | “STROBE / TRIPOD-AI；目标试验方案；双重稳健估计目标；绝对与相对风险；KM 生存分析” | 目标试验模拟 · IPTW（zepid）· g-formula（zepid / 手动实现）· TMLE（zepid.causal.gformula / econml）· 孟德尔随机化（pymr / rpy2+TwoSampleMR）· KM/AFT（lifelines） | 同样使用 pf.etable + 风险差异 / 风险比 / E-value 行 | §A. 流行病学流程 |
| 模式 B — ML 因果推断 | “DML / 元学习器 / 因果森林 / DR-learner；CATE 分布；策略价值” | DML（econml.dml / doubleml）· S/T/X/R/DR-Learner（econml.metalearners / causalml）· GRF 因果森林（econml.grf）· Dragonnet/TARNet/CEVAE（causalml）· BCF（pymc-bart/bcf-py）· 矩阵补全 | pf.etable ML 对比检验 + CATE 图 + 策略价值表 + 保形 PI（mapie） | §B. ML 因果流程 |
如何调用非默认模式（Claude / Agent 会根据用户的措辞选择）：
| 用户说…… | 技能切换到的模式 |
|---|---|
| “运行 DID / IV / RD / 事件研究”，“AER 表格”，“应用微观经济学” | 默认（AER 经济学）— 步骤 1 → 8 |
| “目标试验模拟”，“g-formula”，“IPTW”，“TMLE”，“孟德尔随机化”，“STROBE / TRIPOD”，“公共健康 / 流行病学”，“epi pipeline”，“RWE 研究”，“队列研究”，“病例对照” | 模式 A（流行病学）— §A |
| “DML”，“双重机器学习”，“因果森林”，“元学习器”，“CATE”，“Dragonnet”，“BCF”，“策略学习”，“保形因果推断”，“公平性审计”，“ML 因果推断”，“uplift 建模”，“因果机器学习” | 模式 B（ML 因果推断）— §B |
| “混合”（例如“先估计 DID，再对异质性进行 ML CATE 分析”） | 按顺序使用默认模式 + 模式 B — 每个估计器都会返回一个系数 + SE 配对值，将它们全部放入一个 pf.etable(...)，作为对比检验列 |
三种模式共用相同的步骤 1–4 清理／表 1／诊断框架、相同的步骤 8 导出栈，以及相同的 DAG 优先识别逻辑——切换模式只会改变步骤 5 使用的估计器类别，不会改变论文的整体结构。如果你只需要描述性统计／表 1／平衡性检验，那么步骤 3 中 AER 风格的 tableone / gtsummary 调用在三种模式下都同样适用。
默认输出规范 — 经济学实证论文
这项技能默认采用应用经济学论文惯例。除非用户明确要求单个点估计，否则每次运行都会生成下方完整的、可直接用于发表的输出集。请将其视为步骤 8 的约定——必须执行，不可选择跳过。
必须生成的表格（始终生成）
| # | 表格 | 来源 / 库 | 保存至 |
|---|---|---|---|
| T1 | 汇总统计与平衡性（处理组与对照组，含 SMD / p 值） | pandas.describe + 自定义 table1()（步骤 3） | tables/table1_balance.xlsx + .docx + .tex |
| T2 ★ | 主要结果 — 多列回归 M1→M6（逐步添加控制变量 + FE） | pyfixest.feols × 6 个设定 → pf.etable() / Stargazer（步骤 5–6） | tables/table2_main.xlsx + .docx + .tex |
| T3 | 机制 / 结果变量阶梯 — 相同处理变量，3 个或更多结果变量并排展示 | 对 y ∈ {Y1, Y2, Y3, Y_main} 循环运行 feols → pf.etable | tables/table3_mechanism.xlsx + .docx + .tex |
| T4 | 异质性 — 子组 × 主要系数（性别、年龄、地区等） | 子组 feols × Wald → pf.etable（步骤 7） | tables/table4_heterogeneity.xlsx + .docx + .tex |
| T5 | 稳健性检验汇总 — 替代 SE / 替代聚类 / 替代样本 / 安慰剂，全部放在同一张表中 | feols × 各种设定 → pf.etable（步骤 6） | tables/table5_robustness.xlsx + .docx + .tex |
★ 表 2 是每篇经济学论文的核心。 这是一张多列回归表，带领读者从原始相关性（M1）逐步了解完整设定（M6：双向 FE + 交互 FE + 聚类稳健 SE）。不要将其合并为单列。不要只报告最受关注的系数。逐步变化本身就是可信度论证：如果 M1→M6 单调且稳定，则该设计很可能具有识别能力；如果加入 FE 后结果崩塌，那本身就是结果。
标准 6 列，顺序如下：
- M1 原始双变量回归（
y ~ treat）- M2 + 人口统计特征（
+ age + edu）- M3 + 行业控制变量（
+ tenure / firm_size / industry）- M4 + 个体 FE（
| worker_id）- M5 + 双向 FE（
| worker_id + year）- M6 + 交互 FE（
| worker_id + year + industry^year），并使用聚类稳健 SE
必须生成的图（始终生成）
| # | 图 | 来源 / 库 | 保存至 |
|---|---|---|---|
| F1 | 趋势 / 研究动机 — 处理组与对照组随时间的变化，并标出政策时点 | df.groupby([year, treat])[y].mean().unstack().plot()（步骤 3） | figures/fig1_trend.png（300 dpi，必须导出 PNG）+ .pdf |
| F2 | 事件研究系数及 95% 置信区间，基准期为 –1 | pf.feols(... ~ i(rel_time, ref=-1) ...) → pf.iplot()（步骤 5） | figures/fig2_event_study.png（300 dpi，必须导出 PNG）+ .pdf |
| F3 | M1→M6 各设定的系数图 | 对 6 个拟合模型使用 matplotlib.errorbar（步骤 8） | figures/fig3_coefplot.png（300 dpi，必须导出 PNG）+ .pdf |
| F4 | 稳健性 / 敏感性曲线 — 设定曲线、HonestDiD 或聚类对比 | 设定曲线循环或 honest_did 图 | figures/fig4_sensitivity.png（300 dpi，必须导出 PNG）+ .pdf |
输出文件布局（默认）
project/
├── tables/ table1_balance.xlsx/.docx/.tex table2_main.xlsx/.docx/.tex
│ table3_mechanism.xlsx/.docx/.tex table4_heterogeneity.xlsx/.docx/.tex
│ table5_robustness.xlsx/.docx/.tex
└── figures/ fig1_trend.png(300dpi)+.pdf fig2_event_study.png(300dpi)+.pdf
fig3_coefplot.png(300dpi)+.pdf fig4_sensitivity.png(300dpi)+.pdf
关键输出规则（必须遵守）：
- 图片格式：所有图片必须同时导出 PNG 格式（≥300 dpi） 和 PDF 格式（用于 LaTeX 排版）
- 表格格式：所有回归表格必须同时导出 Excel（.xlsx）、Word（.docx） 和 LaTeX（.tex） 三种格式
- PNG 用于幻灯片、Markdown 文档、邮件等场景；PDF 用于学术论文排版
何时偏离标准要求
- 单次快速估计 — 仅生成相关单元格，但需说明标准交付内容是上述完整集合，并提出可以执行完整生成。
- 研究设计不支持某张图（横截面数据 → 无事件研究）— 跳过该图，并打印说明原因；不得悄然省略。
- 处理组单位数 N=1（合成控制） — 用 SCM 轨迹 + 安慰剂分布替代 F1/F2；T1–T5 仍适用。
必需的库
pip install pandas numpy scipy matplotlib seaborn \
statsmodels linearmodels pyfixest \
rdrobust rddensity \
econml causalml \
stargazer # publication-ready regression tables
# Optional but commonly needed:
pip install missingno # missing-data visualization
pip install pyreadstat # Stata .dta / SPSS .sav import
pip install arch # GARCH, unit-root tests, HAC
pip install pingouin # clean stats tests wrapper
pip install pysynth # synthetic control (N=1 treated)
8 个步骤——规范流程（对应 AER 论文各章节）
┌──────────────────────────────────────────────────────────────────────┐
│ Step −1 Pre-Analysis Plan (PAP) statsmodels.stats.power / mde │
│ Step 0 Sample log + data contract sample_log/asserts/JSON dump │
│ Step 1 Data cleaning missing / outliers / dtype / join │
│ Step 2 Variable construction log / winsorize / std / encode │
│ Step 2.5 Empirical strategy equation × ID assumption + pre-reg│
│ Step 3 Descriptive statistics Table 1 / corr / distribution │
│ Step 3.5 Identification graphics event-study/1st-stage/McCrary/love│
│ Step 4 Diagnostic tests normality / hetero / autocorr / VIF│
│ Step 5 Baseline modeling OLS / panel / IV / DID / RD / SC │
│ Step 6 Robustness battery placebo / subsample / spec curve │
│ Step 7 Further analysis mechanism / heterogeneity / mediation│
│ Step 8 Tables & figures stargazer / coefplot / event study│
└──────────────────────────────────────────────────────────────────────┘
这 8 个步骤对应应用微观经济学领域 AER / QJE / AEJ 论文的规范章节。每个步骤对应论文中的一个章节，并会在磁盘上生成可直接用于论文的产物：
Paper section Step Python moves
─────────────────────────── ───── ────────────────────────────────────────────────
Pre-Analysis Plan −1 statsmodels.stats.power + freeze pap.json
§1. Data 0 sample_log + 5-check data contract → JSON
§1. Data 1 pandas read_*/dropna/dtype/merge(validate=)
§1. Data 2 np.log/np.arcsinh/winsorize/groupby.shift/diff
§1.1 Descriptives (Table 1) 3 df.describe() · table1_with_smd · seaborn
§2. Empirical Strategy 2.5 equation × ID assumption table + pre-reg
§3. Identification graphics 3.5 pf.iplot · binscatter · rdplot · mccrary · love
§3.5 Diagnostics 4 statsmodels.diagnostic + scipy.stats
§4. Main Results (Table 2) 5 pf.feols/IV2SLS/CausalForest · pf.etable / Stargazer
§5. Heterogeneity (Table 3) 7 pf.feols(... + i(.):X) · marginaleffects-py
§6. Mechanisms / Channels 7 Baron-Kenny · econml.dml · outcome ladder
§7. Robustness gauntlet 6 placebo · oster · honestdid · spec_curve · 2-way
§8. Replication package 8 Stargazer.render_latex · pf.etable("docx") · result.json
以下是每个步骤的规范代码。所有示例都围绕同一个案例展开——一个劳动经济学面板数据案例，其中 training（处理变量）影响 log_wage（结果变量），协变量为 age、edu、tenure，面板键为 worker_id / firm_id / year。列名和参数值仅供说明——请替换为用户 DataFrame 中的实际值。只有库名称和调用形式是规范要求。
当某个步骤有许多变体时（例如，交错 DID 有五种不同的估计量；异方差有四种经典检验），SKILL.md 会介绍首选方案，并链接到
references/NN-<topic>.md，供查阅其他方案。如果用户的情况不符合默认方案，请阅读参考文件。
按章节整理的论文级图表清单（需要制作的内容）
现代 AER 论文通常包含 5–7 张图和 3–5 张正文表格，另加一张附录稳健性表格。以下每个步骤都会在磁盘上留下至少一个编号产物。默认文件名假设同时导出 .tex / .docx / .xlsx（Agent 应生成这三种格式，以便合作者在 Word 中编辑、构建系统使用 LaTeX，编辑者在 Excel 中编辑原始数据）。所有图片必须同时保存为 PNG（≥300 dpi）和 PDF 两种格式。
| § | 产物 | Python 基本操作 | 文件名 |
|---|---|---|---|
| §1 | 图 1：原始趋势 / 处理推出情况 | df.groupby([time,treat])[y].mean().unstack().plot() · seaborn.heatmap 用于分期推出情况 | figures/fig1_trend.png(300dpi)+.pdf |
| §1 | 表 1：汇总统计（全样本 / 处理组 / 对照组 + Δ + SMD） | table1(df, by=, cols=)（步骤 3.b）→ 写入 LaTeX/Word/Excel | tables/table1_balance.xlsx/.docx/.tex |
| §3 | 图 2：识别图示（事件研究 / 第一阶段 / McCrary / RD 散点图 / SCM 轨迹） | pf.iplot(es) · binscatter · rdplot · rddensity · synthdid | figures/fig2_event_study.png(300dpi)+.pdf |
| §4 | 表 2：主要结果——逐步加入控制变量 M1→M6 | pf.etable([m1...m6]) · Stargazer([m1.fit ... m6.fit]) | tables/table2_main.xlsx/.docx/.tex |
| §4 | 表 2-bis：设计对比（OLS / IV / DID / DML） | pf.etable([ols, iv, did, dml]) | tables/table2b_designs.xlsx/.docx/.tex |
| §4 | 图 3：各模型设定的系数图 | pf.coefplot([m1...m6], coefs=["training"]) | figures/fig3_coefplot.png(300dpi)+.pdf |
| §5 | 表 3：按子组划分的异质性 | pf.etable(g_full, g_male, g_fem, g_q1...q4) | tables/table3_heterogeneity.xlsx/.docx/.tex |
| §5 | 图 4：剂量—反应 / CATE | econml.CausalForestDML(...).effect() + matplotlib 直方图 | figures/fig4_cate.png(300dpi)+.pdf |
| §6 | 表 4：机制 / 结果阶梯 | 循环对各结果运行 pf.feols → pf.etable | tables/table4_mechanism.xlsx/.docx/.tex |
| §7 | 表 A1：稳健性总表（每项检验一列） | pf.etable([base, no99, balpan, dropearly, wfe, cl2way, logy, ihsy, psm, ebal]) | tables/tableA1_robustness.xlsx/.docx/.tex |
| §7 | 图 5：模型设定曲线 | 手动编写 itertools.product + matplotlib errorbar | figures/fig5_spec_curve.png(300dpi)+.pdf |
| §7 | 图 6：敏感性分析（HonestDiD / Oster / E-value） | HonestDiD · oster_bound · evalue | figures/fig6_sensitivity.png(300dpi)+.pdf |
| §8 | 复现材料包：将所有表格汇总至一份文档 | pf.etable([...] + extra=[...], type="docx") · pylatex / texdoc 风格的多面板文档 | replication/paper_tables.xlsx/.docx/.tex |
上述每个 Python 估计器（通过 R-callout 调用的
pf.feols/IV2SLS/att_gt/CausalForestDML）都会返回结果对象，可直接传入pf.etable(...)/pf.coefplot(...)/Stargazer(...)。不要用df.to_latex()手动编写 LaTeX，也不要直接通过python-docx渲染 Word——pf.etable/Stargazer会自动应用 book-tab 边框、AER 显著性星号和正确的标准误标签。有关更深入的导出方案（LaTeX / Word / Markdown 变体、多面板.docx、完整的gtsummary风格工作流），请参阅references/08-tables-plots.md。
导出指南——在一个区块中导出 LaTeX / Word / Excel
关键规则（必须遵守）：每个表格必须同时导出三种格式——Excel(.xlsx)、Word(.docx)、LaTeX(.tex)。每个图片必须同时保存 PNG(≥300dpi) 和 PDF 两种格式。
分为三个层级，根据范围选择：
| 层级 | 适用场景 | API | 常用 kwargs |
|---|---|---|---|
| 1. 单个多列布局表格 | 导出一个包含递进列的 Table 2 / Table 3 / Table A1 | pf.etable([m1,...,mN], type="tex"/"docx"/"xlsx", file="...", headers=[...], digits=3, signif_code=[0.1,0.05,0.01]) — 或者，对于仅使用 statsmodels 的情况，使用 Stargazer([m1.fit,...,mN.fit]).render_latex() | keep=, drop=, coef_map=, headers=, digits=, signif_code=, fixef_rm=, notes= |
| 2. 多面板论文格式（将 Tables 2 + 3 + A1 + A2 放在一个文件中） | 生成论文表格区块——将主要结果 + 异质性 + 稳健性 + 安慰剂检验整合为一份文档 | 通过重复调用 pf.etable(..., extralines=...) 进行堆叠，或使用 gtsummary 风格的链式表格；若要生成真正的单文件多面板文档，则写入 .tex 后再合并 | 第一个面板：写入；后续面板：追加；使用 LaTeX \section{} 标题包围各面板 |
3. 完整会话包（相当于 Stata collect / R gt） | 生成复现附录，将汇总统计 + 平衡性检验 + 多个回归表 + 标题 + 正文整合到一个文件中 | 使用 pylatex / python-docx / quarto 以编程方式组合，并一次性渲染。或者使用 statsmodels.iolib.summary2.summary_col 快速合并表格。 | 期刊风格模板 + 各节标题 + 脚注宏 |
期刊样式——选择合适的 signif_code 和标准误标签。 AEA 惯例采用 [0.1, 0.05, 0.01]，标准误标签为 "Cluster-robust standard errors in parentheses"。在 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.")
def aer_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)
多面板 .docx / .xlsx 以及 Markdown / Quarto 指南（单文件论文表格包）请参阅 references/08-tables-plots.md。
步骤 −1 —— 分析前计划（数据分析前；AEA RCT Registry 风格）
在接触数据之前，写明 (a) 研究总体，(b) 研究设计，(c) 在计划样本量和 α=0.05、β=0.20 条件下的最小可检测效应（MDE）。将结果保存为 pap.json，以便审稿人确认研究设计在查看数据之前就已具备足够的统计功效。
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>",
}
with open("artifacts/pap.json", "w") as f:
json.dump(pap, f, indent=2)
在步骤 1 之前，将 artifacts/pap.json 提交到代码库。AEA RCT Registry / OSF 预注册工具可将其作为分析计划附件。
步骤 0 —— 样本构建日志与五项检查数据契约
AER §1 Data 部分有三项任务：(a) 描述数据来源，(b) 记录所有样本限制（即“脚注 4”样本日志），(c) 固定面板结构。数据契约可以解决回复信中“样本量神秘缩减”的问题。
0.1 样本构建日志（脚注 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}")
with open("artifacts/sample_construction.json", "w") as f:
json.dump(sample_log, f, indent=2)
将打印出的各行原样粘贴为论文脚注 4。
0.2 五项检查数据契约（继续 / 停止关卡）
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)
如果任何断言触发，停止并在 pandas 中修复。估计器会在后续处理中静默丢弃 NaN 行——这份契约是防范审稿人询问“为什么 N 从表 1 的 12,000 降到表 2 的 9,800？”的最简便保障。
步骤 1 —— 数据清理
更多模式： references/01-data-cleaning.md —— 缺失值处理策略（MCAR/MAR/MNAR）、异常值检测（IQR / z-score / Mahalanobis）、dtype 强制转换、合并注意事项、面板结构验证、去重。
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() > 4
print(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}")
关键原则：使用 pandas 并明确作出决策。绝不要在估计器内部静默删除行——所有行排除都必须在步骤 1 中完成，并打印数量。
步骤 2 —— 变量构建与转换
更多模式： references/02-data-transformation.md —— log/ihs/Box–Cox、缩尾与截尾、组内标准化、独热编码与序数编码及目标编码、交互项、滞后 / 超前算子、一阶差分、使用 CPI 进行平减。
# 2a. Log / IHS transform (skewed positive variables)
df["log_wage"] = np.log(df["wage"].clip(lower=1)) # floor at 1 to avoid -inf
df["ihs_assets"] = np.arcsinh(df["assets"]) # handles zero / negative
# 2b. Winsorize (top/bottom 1%) — reduces outlier influence without dropping rows
from scipy.stats.mstats import winsorize
df["wage_w"] = winsorize(df["wage"], limits=[0.01, 0.01]).data
# 2c. Standardize (z-score) — for interpretability or ML
df["age_std"] = (df["age"] - df["age"].mean()) / df["age"].std()
# 2d. Categorical encoding
df = pd.get_dummies(df, columns=["industry"], prefix="ind", drop_first=True)
# 2e. Interaction & polynomial
df["age_sq"] = df["age"] ** 2
df["training_x_edu"] = df["training"] * df["edu"]
# 2f. Panel operators — first difference, lag, lead, within-group mean
df = df.sort_values(["worker_id", "year"])
df["wage_l1"] = df.groupby("worker_id")["log_wage"].shift(1) # lag
df["wage_f1"] = df.groupby("worker_id")["log_wage"].shift(-1) # lead
df["d_wage"] = df.groupby("worker_id")["log_wage"].diff() # Δy_it
df["wage_mean"] = df.groupby("worker_id")["log_wage"].transform("mean") # within-unit mean
# 2g. Treatment timing variables for staggered DID
df["first_treat_year"] = df.groupby("worker_id")["training"] \
.transform(lambda s: s.idxmax() if s.any() else np.nan)
df["rel_time"] = df["year"] - df["first_treat_year"]
步骤 2.5 — 实证策略（写出方程 + 识别假设）
这是 AER 论文的核心。在编写任何代码之前，明确写出方程并说明识别假设。含糊的识别表述是审稿人拒绝应用论文最常见的原因。将策略保存为 strategy.md，使其成为有日期、受版本控制的文档——而不是看到系数后才事后编写的合理化解释。
方程 × 识别假设 × Python 估计器（决策表）
| 设计 | 估计方程 | 识别假设 | Python 估计器 |
|---|---|---|---|
| 2×2 DID | Y_it = α_i + λ_t + β·D_it + X'γ + ε_it | 给定 X 后满足平行趋势 | `pf.feols("y ~ i(treated, post, ref=0) |
| 事件研究（CS / SA） | Y_it = α_i + λ_t + Σ_{e≠-1} β_e · 1{t-G_i = e} + ε_it | 无预期效应 + 组别-时间平行趋势 | `pf.feols("y ~ i(rel, ref=-1) |
| 2SLS | Y_i = α + β·D_i + X'γ + ε_i; D_i = π·Z_i + X'δ + u_i | 排除限制 + 相关性 + 单调性 | linearmodels.iv.IV2SLS.from_formula("y ~ X + [D ~ Z]", df) |
| Sharp RD | Y_i = α + β·1{X_i ≥ c} + f(X_i) + ε_i（局部多项式） | E[Y(0)|X] 在 c 处连续，且不存在操纵 | rdrobust(y, x, c=0)（+ rddensity） |
| SCM | Ŷ_1t(0) = Σ_j ŵ_j Y_jt，τ_t = Y_1t − Ŷ_1t(0)，适用于 t≥T_0 | 处理前拟合 + 插值有效性 | pysynth · synthdid（R 调用） |
| 基于可观测变量的选择（匹配/IPW/DML） | Y_i = m(X_i) + β·D_i + ε_i（Robinson 偏差剔除） | 无混杂 + 重叠性 | econml.dml.DML · econml.dr.DRLearner · causalml.matching.PSM |
设计选择器（用户不确定时）
┌─ running var + cutoff ───────────────── RDD (rdrobust)
│
├─ exogenous instrument Z ─────────────── IV/2SLS (linearmodels.IV2SLS)
data + question ─┤
├─ pre/post × treat/control ─┬ 2 periods ── 2×2 DID (pf.feols + i())
│ └ staggered ── CS / SA / BJS (att_gt / sunab / did_imputation)
│
├─ 1 treated unit + donor pool + long pre ── SCM (pysynth / synthdid)
│
├─ high-dim X, selection-on-observables ── ML causal (econml.dml / causalml — see §B)
│
└─ none of the above ──────────────────── matching + sensitivity (causalml + evalue)
预注册 strategy.md 模板
from pathlib import Path
strategy = """\
# Empirical Strategy (pre-registration)
**Frozen**: 2026-01-15 (Git SHA: <paste>)
**Population**: manufacturing workers, 2010–2020, balanced panel
**Treatment**: training (binary, staggered adoption)
**Outcome**: log_wage (CPI-deflated 2010 USD)
**Estimand**: ATT on the treated, dynamic horizon -4..+4
## Estimating equation (paste from §2.5 row that matches the design)
log_wage_it = α_i + λ_t + Σ_{e≠-1} β_e · 1{t - G_i = e} + ε_it
## Identifying assumption
1. No anticipation: E[Y_it(0) | t < G_i] = E[Y_it(0) | never-treated]
2. Group-time PT: Δ E[Y_it(0)] is the same across treatment cohorts
## Auto-flagged threats (must defend in §2)
- Selection of G_i on Y_i(0) → bacondecomp + honestdid sensitivity
- Spillover within firm → cluster at firm_id, also try firm_id × year
- Anticipation in pre-period → include lead in event study
## Fallback estimators (Step 6 robustness)
- Sun–Abraham via `pf.feols("y ~ sunab(G, t) | i + t", df)` (or R callout to `fixest::sunab`)
- Borusyak-Jaravel-Spiess via `did_imputation` (R callout)
- Synthetic DID via `synthdid` (R callout)
"""
Path("artifacts/strategy.md").write_text(strategy)
在运行步骤 5 / 步骤 6 之前，将 artifacts/strategy.md 提交到仓库。该文件的 git log 就是分析计划。
步骤 3 — 描述性统计与表 1
更深入的模式：references/03-descriptive-stats.md — 含标准化均值差（SMD）的分层表 1、加权描述性统计、分布比较（ECDF、QQ）、相关性热图 + 聚类、面板平衡检查、时间序列季节性图。
import seaborn as sns
import matplotlib.pyplot as plt
# 3a. Full-sample summary — classic "Table 1, column 1"
stats = df[["log_wage","age","edu","tenure","training"]].describe().T
stats["N"] = df[["log_wage","age","edu","tenure","training"]].notna().sum()
print(stats[["N","mean","std","min","25%","50%","75%","max"]])
# 3b. Stratified Table 1 (treated vs. control + t-test / SMD)
from scipy import stats as sps
def table1(df, by, cols):
rows = []
for c in cols:
t, ctrl = df.loc[df[by]==1, c], df.loc[df[by]==0, c]
smd = (t.mean() - ctrl.mean()) / np.sqrt((t.var()+ctrl.var())/2)
p = sps.ttest_ind(t.dropna(), ctrl.dropna(), equal_var=False).pvalue
rows.append([c, t.mean(), t.std(), ctrl.mean(), ctrl.std(), smd, p])
return pd.DataFrame(rows, columns=["var","treat_mean","treat_sd",
"ctrl_mean","ctrl_sd","SMD","p"])
t1 = table1(df, by="training", cols=["log_wage","age","edu","tenure"])
print(t1.round(3))
# 3c. Correlation heatmap
corr = df[["log_wage","age","edu","tenure","training"]].corr()
sns.heatmap(corr, annot=True, cmap="RdBu_r", center=0, vmin=-1, vmax=1)
plt.tight_layout(); plt.savefig("fig_corr.pdf")
# 3d. Distribution plot — density by treatment status
fig, ax = plt.subplots(figsize=(6,4))
for g, sub in df.groupby("training"):
sub["log_wage"].plot.kde(ax=ax, label=f"training={g}")
ax.set_xlabel("log wage"); ax.legend(); plt.savefig("fig_dist.pdf")
# 3e. Time trends — treated vs. control (the DID motivation plot)
trend = df.groupby(["year","training"])["log_wage"].mean().unstack()
trend.plot(marker="o"); plt.ylabel("mean log wage"); plt.savefig("fig_trend.pdf")
# 3f. Panel balance diagnostic — how many units observed each year?
df.groupby("year")["worker_id"].nunique().plot(kind="bar")
plt.ylabel("# unique workers"); plt.savefig("fig_panel_balance.pdf")
步骤 3.5 — 识别图示（“识别，图示证据”一节）
AER 惯例：识别图应放在回归表之前。 在要求读者相信点估计之前，应先让他们看到平行趋势成立 / 第一阶段较强 / RD 清晰跳跃的图示证据。
3.5.1 事件研究图 + 数值前趋势检验（DID 识别）
处理前系数 ≈ 0（将 −1 期参考期规范化为零）是平行趋势的可视化证据。将图示与数值前趋势检验结合起来，免得审稿人只能凭目测判断。
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" not in 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-.")
3.5.2 第一阶段 F 统计量 + 散点图（IV 识别）
经验法则：OLS 风格推断中，第一阶段 F ≥ 10；AR 等效推断中，F ≥ 23（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 F
print(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 密度检验 + 标准 RD 图
RD 的标志性图形是 rdplot（CCT 风格的分箱散点图，并在两侧分别进行局部多项式拟合），并配合 McCrary 操纵检验。
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 匹配：love plot（匹配前后的标准化差异）
import causalml.matching as cm
from causalml.match import 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：合成控制轨迹 + 差距图
对于合成控制设计，标准的图 2 是处理组与合成组的时间序列，并标注处理时间。
# Pure Python: pysynth or sparseSC; for SDiD use rpy2 callout to R::synthdid
from 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)
识别专项检查（DID 的平行趋势、弱 IV 的 F 统计量、RD 的密度、匹配的共同支持）也会在步骤 5 的估计器中自动运行——此处不要重复数值检验，但务必生成图示：审稿人会先浏览图表。
第 4 步——诊断性统计检验
更深入的模式：references/04-statistical-tests.md——涵盖每种经典检验及其原假设／备择假设／判定规则（Shapiro–Wilk、Kolmogorov–Smirnov、Jarque–Bera、Breusch–Pagan、White、Goldfeld–Quandt、Durbin–Watson、Breusch–Godfrey、Ljung–Box、ADF、KPSS、Phillips–Perron、VIF、条件数、Hausman、Wu–Hausman、Sargan–Hansen）。
在将估计结果视为可靠之前，先进行诊断。 以下 5 类涵盖了 90% 的应用研究。
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}")
判定表（经典经验法则）：
| 检验 | 原假设 | 拒绝后的处理 |
|---|---|---|
| Jarque–Bera / Shapiro | 残差 ~ 正态分布 | N 较大时通常可忽略；若小样本推断很重要，则使用 Bootstrap 置信区间 |
| Breusch–Pagan / White | 误差同方差 | 使用 cov_type="HC3" 或聚类标准误 |
| Durbin–Watson / Breusch–Godfrey | 无自相关 | 使用 HAC（Newey–West）或按单位聚类 |
| VIF > 10 / cond# > 30 | — | 删除／合并共线回归变量 |
| ADF 拒绝 + KPSS 未拒绝 | 序列平稳 | 拟合水平值 |
| ADF 未拒绝 | 单位根 | 一阶差分或进行协整检验 |
第 5 步——基准实证建模（第 4 节：主要结果）
更深入的模式：references/05-modeling.md——涵盖每种经典估计量及其 API：OLS、WLS、GLS、logit/probit、面板 FE / RE / PO、聚类标准误、2SLS / LIML / GMM、DID（2×2、TWFE、事件研究、CS、SA、BJS、SDiD）、RD（锐断点、模糊断点、折点、多切点）、合成控制、PSM / IPW / EB、DML / 因果森林 / DR-Learner。
这是应用论文中最为密集的一节。现代 AER 第 4 节通常包含 2–3 张多元回归表和一张系数图：
- 表 2（主表）：逐步加入控制变量，4–6 列——见下文 模式 A
- 表 2-bis（设计方法对比）：OLS / IV / DID / DML 下的同一系数——模式 B
- 表 2-ter（多结果变量）：同一处理变量，多个结果变量并排展示——模式 C
- 图 3（系数图）：直观汇总不同设定下的 β̂ 和 95% 置信区间
估计量选择路径（请牢记——选错会悄无声息地产生荒谬结果）：
- 无 FE / 单一低维 FE →
smf.ols("y ~ X", df).fit(cov_type="cluster", cov_kwds={"groups": df["i"]})- 高维 FE →
pf.feols("y ~ X | fe1 + fe2", df, vcov={"CRV1":"i"})- 双向聚类 →
pf.feols(..., vcov={"CRV1":"firm_id+year"})- 2SLS / IV →
IV2SLS.from_formula("y ~ X + [D ~ Z]", df).fit(...)或pf.feols("y ~ X | D ~ Z", df)- DID / 事件研究 → SA 使用
pf.feols("y ~ sunab(G, t) | i + t", df)；CS 使用 R 调用did::att_gt
根据识别策略选择估计量（而非追逐“潮流”）：
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)
规范调用（详见 references/05-modeling.md）。下文的八种回归表模式 A–H 是 AER 回归表指南——pf.etable(*models, ...) 是常用工具，相当于 Stata 的 outreg2 / esttab 和 R 的 modelsummary。
5.A 模式 A——逐步加入控制变量（规范的表 2）
各列中的 β̂ 保持稳定 ⇒ 更不必担心估计结果是由可观测变量选择驱动的（Oster 2019 的选择稳定性逻辑；在第 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 惯例：展示所有控制变量（以及截距）。 不要传入
keep=或drop=，以便显示每个参数。仅当有意制作只展示核心系数的表格（交互项形式的异质性分析、IV 第一阶段三列表）时，才使用keep=["training"]；仅当希望出于论文排版考虑隐藏常数项时，才使用drop=["Intercept"]。
5.B 模式 B——设计方法对比（表 2-bis）
展示同一关注系数在多种识别策略下的估计结果。这是 AER 提升可信度的关键做法：不同设计各自采用不同识别假设，但证据相互印证。
from linearmodels.iv import IV2SLS
from econml.dml import LinearDML
from causalml.match import 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 模式 C——多结果变量表（相同的 X，多个 Y）
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 模式 D — 堆叠面板 A / 面板 B 表格
同一模型系列、两个时间跨度（短期 / 长期）或两个样本。使用两次 pf.etable 调用 + LaTeX 拼接，纵向堆叠。
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 模式 E — IV 报告三联表（第一阶段 / 约化式 / 2SLS）
教科书式的 AER IV 表格会将第一阶段、约化式和 2SLS 分列展示，以便读者验证 Wald 比率 = 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 三联表有意聚焦于核心结果：**仅展示 Z + 内生回归变量，以便读者一眼看出 Wald 比率。只有在审稿人要求提供完整系数列表时，才删除
keep=["Z"]。
5.F 模式 F — 通过 did::att_gt / synth / econml.DML 实现因果分析器主结果
对于 DID / SCM / 匹配主结果，现代 Python 估计器会返回自包含的估计值 + 自动安慰剂 / 前趋势 / 重叠诊断。通过精简适配器传入 pf.etable，或直接使用 Stargazer。
# DML with full diagnostics
from 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 模式 G — 子组 pf.etable（表 3，见步骤 7）
每个子组对应一列。详细代码见 §步骤 7 — 异质性。
5.H 模式 H — 稳健性汇总表（表 A1，见步骤 6）
将所有稳健性设定与基准设定并列。详细代码见 §步骤 6.f。
规范的估计器命令（底层原语）
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
步骤 6 — 稳健性检验组
更深入的模式：references/06-robustness.md — 替代设定（增加 / 删除控制变量）、子样本拆分、替代聚类层级、替代结果变量定义、安慰剂检验（虚假处理 / 虚假时点 / 置换）、设定曲线、Oster δ*、随机化推断、留一法、缩尾敏感性分析。
论文中的每个核心结果都需要有稳健性附录。规范的 6 项检验：
# 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))
步骤 7 — 进一步分析（机制 / 异质性 / 中介 / 调节）
更深入的模式：references/07-further-analysis.md — 子组 + Wald 交互检验、用于效应异质性的三重差分（DDD）、Baron–Kenny / Imai 中介分析、带敏感性分析的因果中介分析、调节中介、结果变量阶梯（短期 → 中间 → 最终）、通过因果森林估计 CATE。
# 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
步骤 8 — 发表用表格与图表
此步骤为必做项 — 每次分析运行都会生成 默认输出规范（位于本技能顶部）中定义的全部 5 张必需表格（T1–T5）和全部 4 张必需图表（F1–F4）。不要因为“回归已经跑完了”就跳过步骤 8。没有表格和图表的系数，并非应用经济学传达结果的方式。
更深入的模式：references/08-tables-plots.md — 使用 stargazer 和 pf.etable() 生成回归表格；带置信区间的系数图；事件研究图（带参考线的处理前 / 后系数）；分箱散点图；子组分析森林图；RD 图；导出为 LaTeX / Word / Excel。
# ============================================================
# 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.
交付物清单（在宣布运行完成前核验）：
[ ] tables/table1_balance.tex [ ] figures/fig1_trend.pdf
[ ] tables/table2_main.tex ★ [ ] figures/fig2_event_study.pdf
[ ] tables/table3_mechanism.tex [ ] figures/fig3_coefplot.pdf
[ ] tables/table4_heterogeneity.tex
[ ] tables/table5_robustness.tex [ ] figures/fig4_sensitivity.pdf
[ ] tables/tableA1_robustness.tex [ ] figures/fig5_spec_curve.pdf
[ ] artifacts/sample_construction.json (footnote 4)
[ ] artifacts/data_contract.json
[ ] artifacts/result.json (reproducibility stamp — see 8l)
8l. 可复现性标记
期刊复现办公室（或未来的合著者）复现核心数值所需的单一产物。保存 Python 版本、随机种子、数据集哈希值、基准系数 + 置信区间，以及协议 / 合约的引用路径：
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",
}
with open("artifacts/result.json", "w") as f:
json.dump(stamp, f, indent=2)
将 artifacts/result.json 与论文 PDF 一同提交。审稿人应能运行 python master.py，并逐位一致地复现此 JSON。
§A — 流行病学 / 公共健康模式
当用户的措辞表明属于模式 A（目标试验模拟 / IPTW / TMLE / MR / STROBE / 流行病学 / 公共健康 / RWE / 队列）时，8 个步骤仍然适用，但第 5 步会将 OLS 和固定效应方法替换为双重稳健 + 生存分析 + MR 三件套，交付内容则遵循 STROBE / TRIPOD-AI 规范。第 1–4 步（清理、构建、表 1、诊断）和第 8 步（导出表格/图表）与默认模式完全相同。
库配置（在默认配置之上安装）：
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 队列构建 + 目标试验方案
在接触数据之前编写方案。将其保存为 protocol.yml，并在论文中引用。
# 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 按暴露因素分组的表 1（与默认模式第 3 步相同）
使用第 3 步中相同的 tableone / gtsummary 风格调用，只需设置 groupby="A"（治疗指标）。未测量混杂的 E 值放在脚注中。
A.2 DAG + 倾向评分重叠（正值性检查）
# Estimate PS, plot overlap (positivity), love-plot SMDs
from 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 plot
import 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
A.3 IPTW + g 公式 + TMLE 双重稳健三件套（替换第 5 步）
流行病学中的“AER 表 2”：一张 3 列表格，每列对应 {IPTW-MSM、g 公式、TMLE} 之一，让读者确认双重稳健估计结果一致。
import zepid as ze
from zepid.causal.ipw import IPTW
from zepid.causal.gformula import TimeFixedGFormula
from zepid.causal.doublyrobust import TMLE
# IPTW (marginal structural model)
iptw = IPTW(df, treatment="A", outcome="event_5y")
iptw.treatment_model("age + edu + smoke + bmi + ldl + sbp", print_results=False)
iptw.outcome_model("A", print_results=False)
iptw.fit()
RD_iptw, CI_iptw = iptw.risk_difference, iptw.risk_difference_ci
# g-formula (parametric)
gf = TimeFixedGFormula(df, exposure="A", outcome="event_5y")
gf.outcome_model("A + age + edu + smoke + bmi + ldl + sbp")
gf.fit(treatment="all") ; r1 = gf.marginal_outcome
gf.fit(treatment="none") ; r0 = gf.marginal_outcome
RD_gf = r1 - r0
# TMLE (doubly robust)
tmle = TMLE(df, exposure="A", outcome="event_5y")
tmle.exposure_model("age + edu + smoke + bmi + ldl + sbp")
tmle.outcome_model("A + age + edu + smoke + bmi + ldl + sbp")
tmle.fit()
RD_tmle, CI_tmle = tmle.risk_difference, tmle.risk_difference_ci
# Stack the triplet into one paper table
import pandas as pd
tableA3 = pd.DataFrame({
"Estimator": ["IPTW-MSM", "g-formula", "TMLE"],
"RD": [RD_iptw, RD_gf, RD_tmle],
"95% CI": [CI_iptw, "—", CI_tmle],
})
tableA3.to_latex("tables/tableA3_dr_triplet.tex", index=False, float_format="%.3f")
A.4 生存结局 — KM / Cox / AFT / RMST
from lifelines import KaplanMeierFitter, CoxPHFitter, WeibullAFTFitter
from lifelines.utils import restricted_mean_survival_time
# KM by treatment
fig, ax = plt.subplots()
for a, sub in df.groupby("A"):
KaplanMeierFitter().fit(sub["time_at_risk"], sub["event_5y"], label=f"A={a}").plot_survival_function(ax=ax)
plt.savefig("figures/figA4_km.pdf")
# Cox HR (covariate-adjusted)
cox = CoxPHFitter().fit(df[["time_at_risk","event_5y","A","age","edu","smoke","bmi","ldl","sbp"]],
duration_col="time_at_risk", event_col="event_5y")
HR = cox.hazard_ratios_["A"]; HR_CI = cox.confidence_intervals_.loc["A"].values
# AFT (Weibull) for time-ratio interpretation
aft = WeibullAFTFitter().fit(df[["time_at_risk","event_5y","A","age","edu","smoke","bmi","ldl","sbp"]],
duration_col="time_at_risk", event_col="event_5y")
# RMST contrast at t=5y
rmst1 = restricted_mean_survival_time(df.query("A==1")["time_at_risk"], df.query("A==1")["event_5y"], t=5*365)
rmst0 = restricted_mean_survival_time(df.query("A==0")["time_at_risk"], df.query("A==0")["event_5y"], t=5*365)
A.5 孟德尔随机化（IVW / Egger / 加权中位数三件套）
通过 R 进行双样本 MR 最方便；在 Python 中，可使用 rpy2，或预先导出 beta / SE，再通过配套的 R 脚本调用 TwoSampleMR / MendelianRandomization。
# Pure-Python: pymr (or write the IVW/Egger formulas by hand — they're closed-form)
# Cross-language: rpy2 is the de-facto path for IVW/Egger/weighted median.
import rpy2.robjects as ro
ro.r('''
library(MendelianRandomization)
mr_input <- mr_input(bx=BX, bxse=BXSE, by=BY, byse=BYSE)
ivw <- mr_ivw(mr_input)
egger <- mr_egger(mr_input)
median <- mr_median(mr_input, weighting="weighted")
''')
# Stack the triplet (IVW, Egger, weighted-median) in one tableA5.
A.6 稳健性 — E 值 / 界限 / 主分层
# E-value from zepid (Linden & VanderWeele formula)
from zepid.sensitivity_analysis import e_value
ev = e_value(measure="RR", est=1.45, lcl=1.10, ucl=1.91) # → required strength of unmeasured confounding
print(ev)
A.7 STROBE / TRIPOD-AI 报告核对清单
保存为 replication/strobe_checklist.md，并在投稿前逐项勾选：
[ ] Eligibility criteria + dates (target-trial protocol)
[ ] Adjustment set with DAG justification (A.2)
[ ] Positivity / overlap diagnostic (A.2)
[ ] Doubly-robust triplet (IPTW + g-formula + TMLE) (A.3)
[ ] Risk difference + hazard ratio + RMST (A.3, A.4)
[ ] E-value for unmeasured confounding (A.6)
[ ] Loss-to-follow-up rate + censoring assumption (A.0)
[ ] Pre-registered protocol or analysis plan (A.0)
§B — 机器学习因果推断模式
当用户的措辞表明属于模式 B（DML / 元学习器 / 因果森林 / Dragonnet / BCF / CATE / 策略学习 / 共形因果 / 公平性 / 因果机器学习）时，流程保留默认模式的第 1–4 步和第 8 步，将第 5 步替换为机器学习估计器方法栈，并在第 7 步和第 8 步之间加入 CATE 分布 + 策略价值层。
库配置（在默认配置之上安装）：
pip install econml doubleml causalml dowhy # core estimators
pip install causal-learn cdt # causal discovery (PC / NOTEARS / GES)
pip install mapie # conformal prediction (incl. causal)
pip install fairlearn # fairness audit
pip install policytree-py # discrete policy learning (optional; also via econml.policy)
B.0 训练集/留出集拆分 + 干扰参数学习器方法栈
from sklearn.model_selection import train_test_split
from sklearn.ensemble import GradientBoostingRegressor, GradientBoostingClassifier
train, holdout = train_test_split(df, test_size=0.3, random_state=42, stratify=df["A"])
# Standard nuisance pair: outcome regression Q(X,A) and propensity score g(A|X)
Q_learner = GradientBoostingRegressor(n_estimators=300, max_depth=3)
g_learner = GradientBoostingClassifier(n_estimators=300, max_depth=3)
B.1 DAG / 估计目标声明（可选用 LLM 辅助）
# Causal discovery on observational data — start with PC algorithm
from causallearn.search.ConstraintBased.PC import pc
cg = pc(df[["A","Y","X1","X2","X3","X4"]].to_numpy())
cg.draw_pydot_graph(labels=["A","Y","X1","X2","X3","X4"])
# OR: NOTEARS for a fully differentiable DAG learner
# from cdt.causality.graph import NOTEARS
# notears = NOTEARS().predict(df[["A","Y","X1","X2","X3","X4"]])
# Always sanity-check the recovered DAG against domain knowledge before reading off the adjustment set.
B.2 估计器方法栈 — DML · 元学习器 · 因果森林 · 神经网络 · BCF（替换第 5 步）
机器学习因果推断中的“AER 表 2”：一张横向对比表，每列代表在同一组 (Y, A, X) 数据上使用的一类估计器——读者希望看到 DML、T-learner、X-learner、因果森林和 Dragonnet 对 ATE 得出一致（或不一致）的结果。
from econml.dml import LinearDML, NonParamDML
from econml.metalearners import SLearner, TLearner, XLearner, DomainAdaptationLearner
from econml.dr import DRLearner
from econml.grf import CausalForest
from causalml.inference.tf import DragonNet
# from bcf_py import BayesianCausalForest # if BCF needed
X_train, A_train, Y_train = train[["X1","X2","X3","X4"]], train["A"], train["Y"]
X_hold, A_hold, Y_hold = holdout[["X1","X2","X3","X4"]], holdout["A"], holdout["Y"]
# DML — semi-parametric, doubly robust to nuisance ML
dml = LinearDML(model_y=Q_learner, model_t=g_learner, discrete_treatment=True, cv=5).fit(Y_train, A_train, X=X_train, W=X_train)
ate_dml = dml.ate(X_hold); ci_dml = dml.ate_interval(X_hold, alpha=0.05)
# Meta-learners — drop the same Q/g into S/T/X
sl = SLearner(overall_model=Q_learner).fit(Y_train, A_train, X=X_train); ate_S = sl.ate(X_hold)
tl = TLearner(models=Q_learner).fit(Y_train, A_train, X=X_train); ate_T = tl.ate(X_hold)
xl = XLearner(models=Q_learner, propensity_model=g_learner).fit(Y_train, A_train, X=X_train); ate_X = xl.ate(X_hold)
# DR-Learner — orthogonal CATE
dr = DRLearner(model_propensity=g_learner, model_regression=Q_learner, model_final=GradientBoostingRegressor()).fit(Y_train, A_train, X=X_train)
ate_DR = dr.ate(X_hold)
# Causal forest — non-parametric CATE
cf = CausalForest(n_estimators=500).fit(X=X_train, T=A_train, y=Y_train); cate_cf = cf.predict(X_hold)
ate_cf = cate_cf.mean()
# Dragonnet — joint propensity + outcome neural net
dn = DragonNet(); dn.fit(X_train.values, A_train.values, Y_train.values); ate_dn = dn.predict(X_hold.values).mean()
# Stack the horse-race
import pandas as pd
tableB2 = pd.DataFrame({
"Estimator": ["DML (linear)", "S-learner", "T-learner", "X-learner", "DR-learner", "Causal Forest", "Dragonnet"],
"ATE": [ate_dml, ate_S, ate_T, ate_X, ate_DR, ate_cf, ate_dn],
})
tableB2.to_latex("tables/tableB2_ml_horserace.tex", index=False, float_format="%.4f")
B.3 CATE 分布 + 亚组 CATE 图（扩展第 7 步）
# CATE histogram and quantile plot — heterogeneity from causal forest
import numpy as np
cate = cf.predict(X_hold)
plt.figure(); plt.hist(cate, bins=30); plt.axvline(0, ls="--", color="k"); plt.xlabel("CATE")
plt.savefig("figures/figB3_cate_hist.pdf")
# CATE by quartile of a covariate
df_h = X_hold.assign(cate=cate, age_q=pd.qcut(X_hold["X1"], 4, labels=False))
df_h.groupby("age_q")["cate"].mean().plot(kind="bar"); plt.ylabel("Mean CATE")
plt.savefig("figures/figB3_cate_by_age_q.pdf")
B.4 策略学习 + 离策略评估
from econml.policy import DRPolicyTree # honest discrete policy tree
policy = DRPolicyTree(max_depth=3).fit(Y_train, A_train, X=X_train)
policy.plot() # prose-readable tree of "treat if X1<a and X2>b"
plt.savefig("figures/figB4_policy_tree.pdf")
# Off-policy evaluation (DR-style policy-value)
pred_policy = policy.predict(X_hold)
policy_value_DR = ((Y_hold * (pred_policy == A_hold).astype(int)).mean()
- (Y_hold * (pred_policy != A_hold).astype(int)).mean())
print(f"DR policy value (holdout): {policy_value_DR:.3f}")
B.5 不确定性（共形因果）+ 公平性 + 敏感性分析
from mapie.regression import MapieRegressor
# Conformal prediction interval around CATE — distribution-free coverage guarantee
mapie = MapieRegressor(estimator=GradientBoostingRegressor(), method="plus", cv=10).fit(X_train, cf.predict(X_train))
y_pred, y_pis = mapie.predict(X_hold, alpha=0.1) # 90% conformal PI
# Fairness audit — disparate-impact / equalized-odds for the learned policy
from fairlearn.metrics import MetricFrame, demographic_parity_difference, equalized_odds_difference
sens = X_hold["sensitive_attr"] # e.g. female
mf = MetricFrame(metrics={"acc": (lambda y,y_hat: (y == y_hat).mean())},
y_true=A_hold, y_pred=policy.predict(X_hold), sensitive_features=sens)
print("Demographic parity diff:", demographic_parity_difference(A_hold, policy.predict(X_hold), sensitive_features=sens))
B.6 机器学习因果推断专用报告核对清单
保存为 replication/ml_causal_checklist.md：
[ ] Nuisance learners listed (Q model, g model, hyperparameters, CV folds)
[ ] Cross-fitting / sample-splitting documented (DML K-fold)
[ ] Overlap / propensity diagnostics (B.0 + A.2-style overlap plot)
[ ] CATE summary (mean, SD, quartiles) + heterogeneity p-value
[ ] Policy value with confidence interval (B.4)
[ ] Conformal coverage rate on holdout (B.5)
[ ] Fairness gaps across sensitive attributes (B.5)
[ ] DAG / adjustment set + sensitivity to unmeasured confounding (E-value or Manski bounds)
库选择速查表
| 步骤 | 任务 | 首选库 | 备选方案 |
|---|---|---|---|
| 1 | 数据清洗 | pandas | 对于超过 10M 行的数据，使用 polars |
| 1 | 缺失数据可视化 | missingno | 手动绘制热力图 |
| 2 | 数据转换 | numpy、pandas、scipy.stats.mstats.winsorize | sklearn.preprocessing |
| 3 | 汇总统计量 + Table 1 | pandas、scipy.stats | tableone |
| 3 | 绘图 | matplotlib、seaborn | 交互式图表使用 plotly |
| 4 | 正态性 / 异方差 / 自相关 | statsmodels.stats.*、scipy.stats | pingouin |
| 4 | 平稳性 | statsmodels.tsa.stattools、arch | — |
| 5 | OLS / 面板固定效应 | pyfixest | statsmodels、linearmodels |
| 5 | IV | linearmodels | pyfixest（面板 IV） |
| 5 | DID（2×2、事件研究、TWFE） | pyfixest（+ 用于 CS/SA/BJS/SDiD 的 diff-diff） | 带交互项的 statsmodels |
| 5 | RD | rdrobust、rddensity | — |
| 5 | 合成控制 | pysynth、sdid | 手动使用 scipy.optimize |
| 5 | 匹配 / IPW / EB | causalml、econml、ebal | pymatch |
| 5 | DML / 因果森林 | econml | doubleml |
| 6 | 规格曲线 | 对公式进行自定义循环 | specurve |
| 7 | 中介分析 | pingouin.mediation_analysis | 手动使用 Baron–Kenny 方法 |
| 7 | CATE | econml 因果森林 | — |
| 8 | 回归表 | stargazer、pyfixest.etable | summary_col（statsmodels） |
| 8 | 系数图 / 事件研究 | pyfixest.iplot、matplotlib.errorbar | forestplot |
| 8 | 分箱散点图 | binsreg | 手动进行分位数分箱 |
常见错误（以及正确做法）
| 错误 | 正确做法 |
|---|---|
| 对面板数据进行 OLS，却不使用任何固定效应 | 使用 pyfixest.feols(... | unit + time, ...) |
| 对聚类数据 / 面板数据 / 聚类数较少的数据使用默认（iid）标准误 | 使用 vcov={"CRV1":"cluster_var"}；当聚类数 < 50 时使用野生自助法 |
| 在分期采用情形下使用 TWFE | 使用 CS / SA / BJS，避免负权重偏差 |
| 拟合模型时静默删除行 | 明确执行步骤 1 的清洗；打印数量 |
| 对缺失率较高的协变量一律用均值填补 | 考虑使用 MICE（statsmodels.imputation）或删除该变量 |
| 只报告系数，不报告置信区间 | 始终报告点估计 + 95% 置信区间 + N；绘制系数图 |
| 只用一种规格，不做稳健性检验 | 提供逐步规格（M1–M6）+ 替代标准误 |
| 只报告核心系数（没有 Table 2） | 务必提供包含多列的 M1→M6 主表——这是经济学论文的核心内容，而非摘要中的一句话 |
| 只有系数表，没有任何图 | 经济学研究结果至少需要 F1 趋势图 + F2 事件研究图 + F3 系数图 + F4 敏感性分析——参见默认输出规范 |
| RD 只使用一个带宽 | 展示带宽敏感性 + rdbwselect |
| IV 未报告第一阶段 F 值 | 始终打印 .first_stage（F > 10，理想情况下 > 104） |
| PSM 没有 SMD 平衡表 | 报告匹配前 / 后的 SMD，目标为 |SMD| < 0.1 |
| 未进行 Imai 敏感性分析就解读中介效应 | 对于因果中介分析，在严肃论文中切勿止步于 Baron–Kenny 方法 |
| “流程结束后在 Excel 中制作 Table 1” | 在进行任何回归之前，就在步骤 3 中制作 Table 1 |
| 通过选择规格进行 p 值操纵 | 公开完整的规格曲线，而非最喜欢的那一项 |
典型的 Agent 执行模式
# 1. Clean
df = load_and_clean(raw_path) # Step 1
# 2. Transform
df = construct_vars(df) # Step 2
# 3. Describe
table1 = build_table1(df, by="training") # Step 3
save_corr_heatmap(df, cols) # Step 3
# 4. Diagnose
run_diagnostics(df, y="log_wage", x=covariates) # Step 4
# 5. Model
results_main = fit_baseline(df) # Step 5
# 6. Robustify
results_robust = robustness_battery(df) # Step 6
# 7. Extend
results_het = heterogeneity_analysis(df) # Step 7
results_mech = mechanism_analysis(df) # Step 7
# 8. Export
export_latex_tables(results_main, results_robust, results_het) # Step 8
save_all_figures() # Step 8
每次运行都必须交付 Default Output Spec 中定义的经济学实证论文输出集：
- 5 张表 —
tables/table1_balance.{tex,docx}、tables/table2_main.{tex,docx}★、tables/table3_mechanism.tex、tables/table4_heterogeneity.tex、tables/table5_robustness.tex - 4 张图 —
figures/fig1_trend.{pdf,png}、figures/fig2_event_study.{pdf,png}、figures/fig3_coefplot.{pdf,png}、figures/fig4_sensitivity.{pdf,png} - 一份诊断日志，打印每项 Step 4 测试结果，并统计每条记录的排除情况
- 一份可复现的脚本 / 笔记本，可从原始数据重新生成整个输出集
如果任何交付物缺失或被跳过，必须打印原因（例如，“省略 F2 事件研究：设计为横截面”）。不得悄然省略。
回归表（pf.etable / Stargazer）指南（单页配方索引）
pf.etable(*models, ...) 和 Stargazer([...]) 是所有多回归表背后的两种基本工具。上述八种模式对应如下：
| 模式 | 列间变化的内容 | 步骤 |
|---|---|---|
| A. 逐步增加控制变量 | 协变量集 / 固定效应深度 | 5.A — 表 2 |
| B. 设计策略对比 | 识别策略（OLS / IV / DID / DML / PSM） | 5.B — 表 2-bis |
| C. 多结果变量 | 因变量 Y | 5.C — 表 2-ter |
| D. 堆叠面板 A / B | 时域 / 样本（面板行 × 规格列） | 5.D — 表 2-quater |
| E. IV 报告三联表 | 第一阶段 / 简化式 / 2SLS | 5.E — 表 2-quinto |
| F. 因果分析编排器 | 1 列，完整诊断（CausalForestDML / att_gt / synth） | 5.F |
| G. 子组表 | 子样本（总体 / 女性 / 男性 / Q1…Q4） | 7 — 表 3 |
| H. 稳健性总表 | 堆叠所有稳健性检验 | 6.g — 表 A1 |
AER 期刊风格的 pf.etable 默认设置：
pf.etable(models,
type="tex", file="tables/tableN.tex",
headers=["(1)", "(2)", ...], # column labels
digits=3, signif_code=[0.1, 0.05, 0.01],
notes=("Cluster-robust SE in parentheses. "
"* p<0.10, ** p<0.05, *** p<0.01."))
# Word version: same call with type="docx".
# For multi-source tables (mixing pyfixest + statsmodels + linearmodels), use Stargazer.
图表工厂（Python 中的 12 张标准 AER 图）
| # | 图表 | Python 命令 | 章节 |
|---|---|---|---|
| 1a | 原始趋势（DID 图 1） | df.groupby([time,treat])[y].mean().unstack().plot() | §1 |
| 1b | 处理实施热力图 | seaborn.heatmap(rollout_pivot) | §1 |
| 2a | 事件研究系数 | `pf.iplot(pf.feols("y ~ sunab(G,t) | i+t"))` |
| 2a' | Bacon 权重 | bacondecomp.bacon(...)（Python 移植版）或 rpy2 调用 | §3 |
| 2a'' | CS-DID 动态效应 | rpy2 → did::att_gt/ggdid | §3 |
| 2b | 第一阶段散点图 | binsreg(y=D, x=Z, w=X) | §3（Step 3.5.2） |
| 2c | RD 标准图 | rdrobust.rdplot(y, x, c=0) | §3（Step 3.5.3） |
| 2c' | McCrary 密度图 | rddensity.rdplotdensity(...) | §3 |
| 2d | 匹配平衡图 | causalml/MatchIt 的匹配前后标准化均值差（SMD） | §3（Step 3.5.4） |
| 2e | SCM 轨迹图 | pysynth · SyntheticControlMethods.Synth | §3（Step 3.5.5） |
| 3 | 主要规格的系数图 | pf.coefplot(models, coefs=["training"]) | §4 |
| 4a | 剂量—反应关系 | econml.dml.LinearDML(...).effect() + matplotlib | §5 |
| 4b | CATE 分布 | CausalForestDML(...).effect() + plt.hist | §5 |
| 5 | 规格曲线 | itertools.product + matplotlib.errorbar（见 6.h） | §7 |
| 6 | 敏感性分析仪表板 | HonestDiD + Oster δ + E-value 表 | §7（Step 6.i） |
| 7 | 最终主图 | 特定估计器（rdplot、iplot、Synth.plot） | §8 |
每张图都要同时导出为
.png ≥ 300 dpi（用于幻灯片 / 网页）。在master.py顶部调用一次plt.style.use("default")+seaborn.set_context("paper")，以保持样式一致。
方法目录
经典 OLS / 面板
import statsmodels.formula.api as smf
import pyfixest as pf
from linearmodels.panel import PanelOLS, RandomEffects, FirstDifferenceOLS
smf.ols ("y ~ X", df).fit(cov_type="cluster", cov_kwds={"groups": df["i"]}) # OLS
pf.feols ("y ~ X | fe1", df, vcov={"CRV1": "i"}) # OLS + 1 FE
pf.feols ("y ~ X | fe1 + fe2", df, vcov={"CRV1": "i"}) # HD FE workhorse
pf.feols ("y ~ X | fe1 + fe2", df, vcov={"CRV1": "fe1+fe2"}) # 2-way cluster
pf.fepois ("count ~ X | fe1 + fe2", df, vcov={"CRV1": "i"}) # Poisson + FE
pf.feglm ("y ~ X | fe1", df, family="logit", vcov={"CRV1": "i"}) # Logit + FE
PanelOLS.from_formula("y ~ X + EntityEffects + TimeEffects", df).fit(cov_type="clustered")
RandomEffects.from_formula("y ~ X", df).fit() # RE (Hausman)
双重差分
# 2×2 DID
pf.feols("y ~ i(treated, post, ref=0) | i + t", df, vcov={"CRV1":"i"})
# Staggered — modern stack (Sun-Abraham via pyfixest::sunab)
pf.feols("y ~ sunab(first_treat, year) | i + year", df, vcov={"CRV1":"i"})
# Callaway-Sant'Anna — R callout via rpy2
import rpy2.robjects as ro; from rpy2.robjects import pandas2ri; pandas2ri.activate()
ro.r('library(did); cs <- att_gt(yname="y", tname="t", idname="i", gname="G", data=df)')
# Borusyak-Jaravel-Spiess — also R callout (didimputation) or Python `did_imputation` if installed
# Synthetic DID — R callout to `synthdid`
工具变量 / 2SLS
from linearmodels.iv import IV2SLS, IVLIML, IVGMM
IV2SLS.from_formula("y ~ 1 + X + [D ~ Z]", df).fit(cov_type="clustered", clusters=df["firm_id"])
# Or via pyfixest:
pf.feols("y ~ X | D ~ Z", df, vcov={"CRV1":"firm_id"})
# Weak-IV diagnostics: ivreg2 has CD/KP, linearmodels reports first_stage().diagnostics
回归不连续
from rdrobust import rdrobust, rdplot, rdbwselect
from rddensity import rddensity
rdrobust(y=df["y"], x=df["x"], c=0, kernel="triangular", bwselect="mserd") # Sharp RD
rdrobust(y=df["y"], x=df["x"], c=0, fuzzy=df["D"]) # Fuzzy RD
rddensity(X=df["x"], c=0) # McCrary
rdplot (y=df["y"], x=df["x"], c=0)
匹配 / 重新加权
import causalml.match as cm
from causalml.match import NearestNeighborMatch
NearestNeighborMatch(replace=False, ratio=1).match(df, treatment_col="D", score_cols=["X1","X2"])
# Entropy balancing — econml or hand-rolled via scipy.optimize
合成控制
from SyntheticControlMethods import Synth
Synth(df, "y", "unit", "time", treatment_period=2015, treated_unit=1)
# SDiD via rpy2 callout to R::synthdid
机器学习因果推断（模式 B — 见 §B）
from econml.dml import LinearDML, CausalForestDML, NonParamDML
from econml.dr import DRLearner
from econml.metalearners import SLearner, TLearner, XLearner
from econml.iv.dml import OrthoIV
# Full pipeline in §B.
稳健性、敏感性与推断
# Two-way cluster: vcov={"CRV1":"firm_id+year"} in pyfixest
# Wild cluster bootstrap: rpy2 → R::fwildclusterboot, or pyfixest's wildboottest()
fit.wildboottest(param="training", B=9999, seed=42)
# Randomization inference: rpy2 → R::ri2 or hand-roll via permutation
# Romano-Wolf: rpy2 → R::wyoung; or hand-roll via stepdown
# Oster δ: see references/06-robustness.md §6.6
# HonestDiD: rpy2 → R::HonestDiD
# E-value: hand-rolled VanderWeele-Ding formula (see 6.i)
生存分析 / 流行病学（模式 A — 见 §A）
from lifelines import KaplanMeierFitter, CoxPHFitter, WeibullAFTFitter
KaplanMeierFitter().fit(df["time"], df["event"]).plot() # KM
CoxPHFitter().fit(df, duration_col="time", event_col="event") # Cox
WeibullAFTFitter().fit(df, duration_col="time", event_col="event") # AFT
# RMST: lifelines.utils.restricted_mean_survival_time
import zepid
zepid.causal.gformula.SurvivalGFormula(...) # g-formula
zepid.causal.ipw.IPTW(...) # IPTW
zepid.sensitivity_analyses.e_value(...) # E-value
何时转交给其他技能
- Agent 原生单次导入工作流（
import statspai as sp）→00-StatsPAI_skill。 - 因果推断：Mixtape 风格代码模板（Python/R/Stata 并列）→
10-Jill0099-causal-inference-mixtape。 - DSGE / HANK 数值宏观经济学 →
20-wenddymacro-python-econ-skill。 - 纯 pyfixest 参考（每个
feols关键字参数）→40-py-econometrics-pyfixest。 - 分析完成后的论文撰写 / LaTeX 起草 → 本仓库中面向写作的技能（
04-*-scientific-writer、08-*-web-latex等）。
本技能的职责止于步骤 8——润色后的表格和图表。论文起草不在职责范围内。
步骤 8.5 — 交接至论文手稿（退出契约）
本技能的职责仍止于步骤 8。不过，“止于此”不应意味着“走入死胡同”： 写作阶段需要一份机器可读的契约，而不是一堆零散文件。 在结束之前，请在表格和图表旁生成两项成果。
1. exhibits_index.md — 每个编号的展项占一行：
| 展项 | 文件 | 支持的论点 | 生成方式 |
|---|---|---|---|
| 表 2 | tables/table2_main.{tex,docx,xlsx} | 政策对产出对数的主要影响 | analysis.py §baseline |
| 图 2 | figures/fig2_event_study.{png,pdf} | 不存在差异化的前趋势 | analysis.py §event-study |
2. results_summary.json — 将每个核心估计值记录为扁平记录：
{"baseline": {"coef": 0.1234, "se": 0.0412, "p": 0.0031, "n": 12345, "spec": "TWFE + region-year FE"}}
为什么需要两项成果，又为什么采用这些格式：
- 手稿阶段通过包含指令引用展项，而不是重新输入数字——在 LaTeX 中使用
\input{tables/table2_main}，在 Markdown 中使用{{ include: table2_main }}。索引会告诉写作者哪个文件支持哪项论点。 - 论文中打印的每个数字都必须能追溯到结果文件，并且精度要与展示时一致。
results_summary.json能让这一点得到核验，而非仅凭信任。 如果某个系数只存在于渲染后的表格中，就无法验证引用它的正文是否准确。
然后进行交接：
- 完整论文 → 全程使用 Word
.docx/ LaTeX →skills/69-Paper-WorkFlow/。它负责协调第 5–9 阶段（起草 → 润色 → 去除 AI 生成痕迹 → 审阅 → 投稿），其scripts/assemble_manuscript_docx.py会根据这两项成果，组装出09_submission/main.docx——包含正文、上述每个展项、图表和参考文献列表的单个文件。在其第 0 阶段选择正文格式：交付物为 Word 时使用 Markdown（Markdown →.docx保真度高；LaTeX →.docx会有损），投稿渠道接受.tex时则使用 LaTeX。 - 仅起草正文 → 本仓库中的写作技能（
04-*-scientific-writer、35-*-academic-writing-skills）。 - 仅进行 LaTeX / Quarto 排版 →
08-ndpvt-web-latex-document-skill、60-regisely-superpapers。
评论
支持 Markdown