Full-empirical-analysis-skill について
応用経済学、疫学、因果機械学習を対象とする、Python向けの8段階の実証分析スキルです。ライブラリ呼び出しの可視性、診断、頑健性分析を重視し、論文向けの表や図を作成します。
主な機能
- 018段階の実証分析ワークフロー
- 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）を直接呼び出すため、エージェントも、エージェントのコードを読むユーザーも、処理の全体を把握でき、どのコンポーネントでも差し替えられます。
関連スキル：単一のインポートで使えるエージェント向け DSL（import statspai as sp）を希望する場合は、代わりに 00-StatsPAI_skill を案内してください。このスキルは正反対の考え方に基づいています。すべてを明示し、すべてを検査可能にし、診断はすべて手動で実行し、プロットはすべてユーザーが形作ります。
基本方針
- 従来型のスタックを使い、魔法は使いません。 エージェントはすべての行を読み、どのライブラリ、推定量、標準誤差の種類が使われているかを正確に把握できる必要があります。
- 推定だけでなく、パイプライン全体を扱います。 実際の論文に費やす時間の80％は、ステップ1〜4および6〜8です。このスキルでは、それらを後回しにせず、主要な作業として扱います。
- 豊富な出力を生成します。 各ステップで少なくとも1つの表または図を生成します。点推定値だけを単独で出すことはありません。
- 段階的に詳細を開示します。 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）、推定対象、識別デザイン、主なアウトカム／処置、パッケージスタック、必要な成果物を記録します。 - 変更範囲を限定：一度に変更する判断事項は1つにします（サンプルルール、推定量、クラスタリング、エクスポート形式、頑健性チェックなど）。検証に合格できる最小限のパッチを優先します。
- 選択用データ分割の規律：ユーザーが直面している問題、または要求した成果物を選択用データ分割として扱います。別のアウトカム、サンプル期間、推定量の種類、またはエクスポート先を少なくとも1つ、保留ゲート用に確保します。
- 保留ゲート：コードを実行する前にチェック項目を定義します。行数、キーの一意性、処置のサポート、欠損率のしきい値、想定される表／図のファイル、さらに変更によって破損してはならない、主な関心対象ではない頑健性チェックまたは仕様を含めます。
- 却下バッファー：候補仕様がゲートに失敗した場合は、その失敗、コードの差分、ゲートの出力を
analysis_log.mdに記録します。最後に合格した仕様に戻し、未検証の同じパターンを再試行しないでください。 - 低速／メタ更新：タスクの最後に、その経過から
accepted_patternsとrejected_patternsを書き出します。1回の合格だけを根拠に、標準プロジェクトテンプレートの適用範囲を広げないでください。 - 検証後にのみ昇格：単発の修正を再利用可能なプロジェクト用定型コードにするのは、現在のデータと、少なくとも1つの別のアウトカム／サンプル／仕様で検証に合格した後に限ります。
3 つの領域モード（デフォルト＝AER 経済学、代替モード＝疫学と ML 因果推論）
上記のデフォルトのプレイブックは、AER スタイルの応用計量経済学、つまり AEA の慣例に基づいています。推定式を明記し、識別仮定を示し、設計を競合させて比較し、あらゆる頑健性検証を行います。このスキルには、他の 2 つの主要な因果推論の伝統に対応する並行サブパイプラインも用意されています。それぞれ、同じ Steps 1–4（クリーニング／構築／記述統計／診断）と Step 8（表と図）を再利用し、ライブラリを切り替えるのは Step 5（推定量）と Step 6/7（頑健性／メカニズム）だけです。
| モード | 読者向けの慣例 | Step-5 推定量スタック | 報告スタック | 移動先 |
|---|---|---|---|---|
| デフォルト — 応用経済学（AER / QJE / AEJ） | 「式＋識別仮定＋設計の競合比較を示す。統制変数を明示し、標準誤差はクラスター化する」 | DID / IV / RD / SCM / マッチング / pyfixest.feols HDFE | AER の標準スタイルによる複数列の pf.etable / Stargazer ＋ 8 セクションの論文レイアウト | Steps 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 値の行 | §A. 疫学パイプライン |
| モード B — ML 因果推論 | 「DML / メタラーナー / 因果フォレスト / DR ラーナー、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 プロット ＋ポリシー価値表 ＋適合予測区間（mapie） | §B. ML 因果推論パイプライン |
デフォルト以外のモードを呼び出す方法（Claude / エージェントがユーザーの表現から判断します）：
| ユーザーの発言... | スキルが切り替えるモード |
|---|---|
| 「DID / IV / RD / イベントスタディを実行」「AER の表」「応用ミクロ」 | デフォルト（AER 経済学） — Steps 1 → 8 |
| 「ターゲットトライアルのエミュレーション」「g-formula」「IPTW」「TMLE」「メンデルランダム化」「STROBE / TRIPOD」「公共健康 / 流行病学」「疫学パイプライン」「RWE 研究」「コホート研究」「ケースコントロール」 | モード A（疫学） — §A |
| 「DML」「ダブル機械学習」「因果フォレスト」「メタラーナー」「CATE」「Dragonnet」「BCF」「ポリシー学習」「適合因果推論」「公平性監査」「ML 因果推論」「アップリフトモデリング」「因果機械学習」 | モード B（ML 因果推論） — §B |
| 「組み合わせ」（例：「DID を推定してから、異質性について ML で CATE を推定」） | デフォルト＋モード B を順に実行 — 各推定量は係数と標準誤差のペアを返すため、それらをすべて 1 つの pf.etable(...) に入れて競合比較の列を作成します |
3 つのモードでは、同じ Step 1–4 のクリーニング／Table 1／診断の枠組み、同じ Step 8 のエクスポートスタック、同じ DAG 優先の識別ロジックを共有しています。切り替えによって変わるのは、使用する Step 5 の推定量ファミリーだけであり、論文の構成は変わりません。記述統計／Table 1／バランスチェックだけが必要な場合、Step 3 の AER スタイルの tableone / gtsummary 呼び出しは、3 つのモードすべてで同じように機能します。
デフォルトの出力仕様 — 経済学の実証論文
このスキルは、デフォルトで応用経済学の論文形式を採用します。ユーザーが単一の点推定値を明示的に求めない限り、実行のたびに、以下の出版可能な完全な出力一式を生成します。これは Step 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／代替クラスター／代替サンプル／プラセボを1 つの表にまとめる | feols × 各種バリエーション → pf.etable（ステップ 6） | tables/table5_robustness.xlsx + .docx + .tex |
★ 表 2 は、すべての経済学論文の中心となる表です。 読者を生の相関（M1）から完全に特定された設計（M6：2 方向 FE + 交互作用を含む FE + クラスター頑健標準誤差）まで導く、複数列の回帰表です。これを単一列にまとめてはいけません。見出しとなる係数だけを報告してはいけません。この段階的な変化そのものが信頼性の論拠です。M1→M6 が単調で安定していれば、設計によって因果効果を特定できている可能性があります。FE を追加して結果が消失するなら、それ自体が結果です。
標準的な 6 列の順序：
- M1 生の二変量回帰（
y ~ treat）- M2 + 人口統計学的変数（
+ age + edu）- M3 + 産業セクターの統制変数（
+ tenure / firm_size / industry）- M4 + 単位 FE（
| worker_id）- M5 + 2 方向 FE（
| worker_id + year）- M6 + 交互作用を含む FE（
| worker_id + year + industry^year）、クラスター頑健標準誤差を使用
必須の図（常に作成）
| # | 図 | ソース／ライブラリ | 保存先 |
|---|---|---|---|
| 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、またはクラスター比較 | spec_curve のループまたは honest_did のプロット（ステップ 6） | 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） の 3 形式すべてで出力する必要があります
- PNG はスライド、Markdown 文書、メールなどに使用し、PDF は学術論文の組版に使用します
例外を認める場合
- 簡易的な推定を 1 回だけ行う場合 — 関連するセルだけを作成します。ただし、標準的な成果物は上記の一式であることを伝え、実行を提案してください。
- 設計上、図を作成できない場合（横断面データ → イベントスタディなし）— 理由を説明するメモを出力して省略してください。黙って省略してはいけません。
- 処置単位が N=1 の場合（合成コントロール） — F1／F2 を SCM の軌跡 + プラセボ分布に置き換えてください。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 論文の標準的なセクションに対応しています。各ステップは論文の 1 セクションに相当し、ディスク上に論文でそのまま使える成果物を出力します。
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 には 5 つの異なる推定量があり、異分散性には古典的な検定が 4 つあります）、SKILL.md では最初に使う方法を示し、残りについては
references/NN-<topic>.mdを参照するよう案内しています。ユーザーのケースが既定の方法に当てはまらない場合は、リファレンスファイルを読んでください。
セクション別の論文掲載用図表一覧（作成するもの）
現代のAER論文では、図が5～7点、主要な表が3～5点に加えて、付録の頑健性検証表が必要です。以下の各ステップでは、少なくとも1つの番号付き成果物がディスクに保存されます。デフォルトのファイル名は、.tex／.docx／.xlsxを並行して出力することを前提としています（エージェントは3種類すべてを作成してください。共著者は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：頑健性検証の総括（各検証を1列に配置） | 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 | 再現用バンドル：すべての表を1つの文書にまとめる | pf.etable([...] + extra=[...], type="docx") · pylatex／texdoc形式の複数パネル | replication/paper_tables.xlsx/.docx/.tex |
上記の各Python推定器（
pf.feols／IV2SLS／R呼び出し経由のatt_gt／CausalForestDML）は、pf.etable(...)／pf.coefplot(...)／Stargazer(...)にそのまま渡せる結果オブジェクトを返します。df.to_latex()からLaTeXを手作業で作成したり、python-docxでWordを直接生成したりしないでください。pf.etable／Stargazerは、ブックタブ形式の罫線、AERの有意性を示す星印、適切な標準誤差のラベルを自動で適用します。より詳しい出力方法（LaTeX／Word／Markdownの各形式、複数パネルの.docx、gtsummary形式の一連の処理）については、references/08-tables-plots.mdを参照してください。
エクスポートクックブック — LaTeX / Word / Excel を1ブロックで
重要なルール（必ず遵守）：すべての表は3形式すべてで同時にエクスポートする必要があります — Excel(.xlsx)、Word(.docx)、LaTeX(.tex)。すべての画像はPNG(≥300dpi)とPDFの2形式で同時に保存する必要があります。
スコープによって選ぶ3つのティア：
| ティア | 使用する場面 | API | 主要な kwargs |
|---|---|---|---|
| 1. 単一の複数列テーブル | 列を段階的に増やしながら 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 を1ファイルに） | 論文用テーブルブロック を生成する — メイン + 異質性 + 頑健性 + プラセボを単一ドキュメントとして | pf.etable(..., extralines=...) 呼び出しを繰り返してスタックするか、gtsummary 形式の連鎖テーブルを使用します。真の単一ファイル・マルチパネルにするには、.tex に書き込んでから連結します | 最初のパネル：書き込み、以降：追記、LaTeX の \section{} ヘッダーで囲む |
3. フルセッションバンドル（Stata の collect / R の gt に相当） | 要約統計量 + バランス + 複数の回帰テーブル + 見出し + 本文を 1つ のファイルに混在させた複製付録 | pylatex / python-docx / quarto でプログラム的に構成します — 一度レンダリングするだけです。または、テーブルを手早く連結するには statsmodels.iolib.summary2.summary_col を使用します。 | ジャーナル形式のテンプレート + セクションごとの見出し + 脚注マクロ |
ジャーナルスタイル — 適切な signif_code と SE ラベルを選ぶ。 AEA の慣例は [0.1, 0.05, 0.01] と SE ラベル「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 を参照してください。
Step −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)
Step 1 の 前に artifacts/pap.json をリポジトリにコミットします。AEA RCT Registry / OSF の事前登録ツールは、これを分析計画の展示として受理します。
Step 0 — サンプル構築ログ & 5項目チェックのデータ契約
AER の §1 Data セクションには3つの役割があります：(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 5項目チェックのデータ契約（go / no-go ゲート）
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 が Table 1 と Table 2 の間で 12,000 から 9,800 に減ったのか？」という査読者の質問に対する最も安価な保険です。
Step 1 — データクリーニング
より詳しいパターン：references/01-data-cleaning.md — 欠損値戦略（MCAR/MAR/MNAR）、外れ値検出（IQR / z-score / マハラノビス）、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 + 明示的な意思決定。推定器の内部で行を黙って落とすことは 決して しません — すべての行の除外は Step 1 で、出力されたカウントとともに行われます。
Step 2 — 変数の構築と変換
より詳しいパターン：references/02-data-transformation.md — log/ihs/Box–Cox、winsorizing 対 trimming、グループ内標準化、one-hot 対 ordinal 対 target エンコーディング、交互作用項、ラグ/リード演算子、一階差分、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"]
Step 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) |
| シャープRD | Y_i = α + β·1{X_i ≥ c} + f(X_i) + ε_i（局所多項式） | cにおけるE[Y(0)|X]の連続性、操作なし | 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)
Step 5 / Step 6を実行する前に、リポジトリ内で artifacts/strategy.md をコミットしてください。このファイルの git log が、そのまま分析計画になります。
Step 3 — 記述統計とTable 1
より詳細なパターンについては、references/03-descriptive-stats.mdを参照してください。SMDを含む層別Table 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")
Step 3.5 — 識別グラフィックス（セクション「識別、グラフィカルな証拠」）
AERの慣例では、識別図を回帰表より先に示します。 点推定値を信頼するよう求める前に、PTが成り立つ／第一段階が強い／RDで明確なジャンプがあることを示すグラフィカルな証拠を読者に見せてください。
3.5.1 イベントスタディ図＋数値による事前トレンド検定（DIDの識別）
介入前の係数がほぼ0（−1の参照期間を0に正規化）は、平行トレンドを示す視覚的な証拠です。査読者が目視で判断しなくて済むよう、図と数値による事前トレンド検定を組み合わせてください。
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 マッチング：ラブプロット（標準化差の実施前・実施後比較）
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：合成コントロールの推移＋ギャッププロット
合成コントロール設計における標準的なFigure 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のPT、弱いIVのF統計量、RDの密度、マッチングの共通サポート）は、Step-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 が大きい場合は通常無視してよい; 小標本での推論が重要な場合はブートストラップ信頼区間を使用 |
| Breusch–Pagan / White | 均一分散の誤差 | cov_type="HC3" またはクラスター標準誤差を使用 |
| Durbin–Watson / Breusch–Godfrey | 自己相関なし | HAC(Newey–West)または単位ごとのクラスタリングを使用 |
| VIF > 10 / 条件数 > 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(sharp、fuzzy、kink、マルチカットオフ)、Synthetic Control、PSM / IPW / EB、DML / causal forest / DR-Learner。
ここは応用論文の中で最も密度の高いセクションです。現代の AER の §4 は通常、2〜3 個の多重回帰テーブルと 1 個の係数プロット を含みます:
- テーブル 2(メイン): 段階的なコントロール変数、4〜6 列 — 以下の パターン A
- テーブル 2-bis(デザインの競争): OLS / IV / DID / DML の下での同じ係数 — パターン B
- テーブル 2-ter(複数アウトカム): 同じ処置、複数のアウトカムを並べて表示 — パターン C
- 図 3(coefplot): 各仕様にわたる β̂ と 95% CI の視覚的要約
推定量のルーティング(これは覚えてください — 間違えると気づかないうちに無意味な結果を生み出します):
- 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 にはdid::att_gtへの R 呼び出し
推定量は識別戦略によって選んでください(「流行っているから」ではなく):
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)。以下の 8 つの回帰テーブルパターン 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=も渡さないでください。keep=["training"]を使うのは、焦点となる係数のみのテーブルが意図的な場合(交互作用形式の異質性、IV の第一段階の三つ組)だけです;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 の表
同じモデルファミリーで、2つの期間（短期／長期）または2つのサンプルを扱います。2つの 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 報告用の3列構成（第一段階／誘導形／2SLS）
教科書的な AER の IV 表では、読者が Wald 比＝誘導形／第一段階を検証できるよう、第一段階、誘導形、2SLSを3列で示します。
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 の3列構成では、対象を意図的に絞ります。**読者が Wald 比を一目で確認できるよう、Z と内生変数のみを示します。査読者から全係数の一覧を求められた場合に限り、
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を参照）
サブグループごとに1列を設けます。詳しいコードは「§ステップ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 δ*、ランダム化推論、1つ抜き分析、ウィンザー化の感度分析を扱います。
論文で主要な結果として示すものにはすべて、ロバストネス検証の付録が必要です。標準となる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)
論文の PDF と一緒に artifacts/result.json をコミットしてください。査読者が python master.py を実行し、ビット単位で完全に一致するこの JSON を再現できる必要があります。
§A — 疫学／公衆衛生モード
ユーザーの表現からモード A（ターゲット試験のエミュレーション／IPTW／TMLE／MR／STROBE／流行病学／公衆衛生／RWE／コホート）が示される場合も、8 つのステップは引き続き適用されます。ただし、ステップ 5 では OLS と FE のスタックを、二重に頑健な手法＋生存時間解析＋MR の 3 点セットに置き換え、成果物は 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-formula＋TMLE の二重に頑健な 3 点セット（ステップ 5 の置き換え）
疫学における「AER 表 2」です。3 列の表を作成し、各列に {IPTW-MSM、g-formula、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／加重中央値の 3 点セット）
2 標本 MR は R を使うのが最も手軽です。Python からは rpy2 を使うか、ベータ値／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 — ML 因果推論モード
ユーザーの表現からモード B（DML／メタラーナー／因果フォレスト／Dragonnet／BCF／CATE／ポリシー学習／コンフォーマル因果推論／公平性／因果機械学習）が示される場合、パイプラインではデフォルトモードのステップ 1〜4 とステップ 8 を維持し、ステップ 5 を ML 推定器のスタックに置き換え、ステップ 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 の置き換え）
ML 因果推論における「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 ML 因果推論専用の報告チェックリスト
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 | 1,000万行を超える場合は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を必ず報告し、係数プロットを作成する |
| 頑健性検証なしで仕様が1つだけ | 段階的な仕様（M1～M6）と代替の標準誤差を提示する |
| 主要な係数だけを報告し、Table 2を示さない | 経済学論文の中心となるのは要旨の一文ではなく、複数列のM1→M6メインテーブルです。必ず提示してください |
| 係数表だけで図がない | 経済学の結果には、F1トレンド＋F2イベントスタディ＋F3係数プロット＋F4感度分析が最低限必要です。Default Output Specを参照してください |
| RDで単一のバンド幅だけを使う | バンド幅の感度分析とrdbwselectを示す |
| 第一段階のF統計量なしでIVを実施する | .first_stageを必ず出力する（F > 10、できれば > 104） |
| SMDバランス表なしでPSMを実施する | 事前／事後のSMDを報告し、|SMD| < 0.1を目標にする |
| Imaiの感度分析なしで媒介分析を解釈する | 因果媒介分析では、真剣な論文でBaron–Kenny法だけで終わらせてはいけません |
| 「パイプライン完了後にExcelでTable 1を作成する」 | 回帰の前に、ステップ3でTable 1を作成する |
| 仕様の選択によるpハッキング | お気に入りのセルではなく、仕様曲線の全体を開示する |
標準的なエージェント実行パターン
# 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）クックブック（1 ページのレシピ一覧）
pf.etable(*models, ...) と Stargazer([...]) は、複数回帰表すべての基盤となる 2 つのプリミティブです。上記の 8 つのパターンは、次のように対応します。
| パターン | 列ごとに異なる要素 | Step |
|---|---|---|
| A. 段階的な統制変数 | 共変量セット／固定効果の深さ | 5.A — Table 2 |
| B. 研究デザインの比較 | 識別戦略（OLS / IV / DID / DML / PSM） | 5.B — Table 2-bis |
| C. 複数のアウトカム | 従属変数 Y | 5.C — Table 2-ter |
| D. Panel A / B の積み重ね | 期間／サンプル（パネルの行 × 仕様の列） | 5.D — Table 2-quater |
| E. IV 報告の 3 点セット | 第一段階／縮約形／2SLS | 5.E — Table 2-quinto |
| F. 因果推論オーケストレーター | 1 列、完全な診断（CausalForestDML / att_gt / synth） | 5.F |
| G. サブグループ表 | サブサンプル（全体／女性／男性／Q1…Q4） | 7 — Table 3 |
| H. 頑健性チェックの一覧 | すべての頑健性チェックを積み重ねたもの | 6.g — Table 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 で作成する AER 標準図 12 点）
| # | 図 | Python コマンド | セクション |
|---|---|---|---|
| 1a | 生データのトレンド（DID Figure 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 値の表 | §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
他のスキルに引き継ぐタイミング
- エージェントネイティブな単一インポートのワークフロー（
import statspai as sp）→00-StatsPAI_skill。 - 因果推論：The 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 までです。ただし、「終了」が「行き止まり」を意味する必要はありません。 執筆段階では、ばらばらのファイルが入ったフォルダーではなく、機械可読な契約が必要です。 作業を終える前に、表と図の横に次の 2 つの成果物を出力してください。
1. exhibits_index.md — 番号付き成果物ごとに 1 行：
| 成果物 | ファイル | 裏付ける主張 | 作成元 |
|---|---|---|---|
| 表 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/。Stage 5～9（草稿 → 推敲 → AI 生成らしさの除去 → 査読 → 投稿）を統括します。また、scripts/assemble_manuscript_docx.pyはこの 2 つの成果物だけを元に、本文、上記のすべての成果物、図、参考文献一覧をまとめた09_submission/main.docxを作成します。Stage 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 対応