Full-empirical-analysis-skill-R について
応用経済学者などの実証研究者向けに、データ処理、因果分析、報告を8段階のワークフローにまとめたRスキルです。経済学、疫学、因果機械学習のモードを備え、再現可能なR・Quartoの実践と表・図の出力を定めています。
主な機能
- 018段階の実証分析ワークフロー
- 02経済学・疫学・因果機械学習のモード
- 03tidyverseとfixestを使うR手法
- 04再現可能なスクリプトとQuartoワークフロー
- 05論文用の表と図の出力
- 06分析変更の検証ゲート
活用シーン
- 01応用経済学論文の分析構築
- 02疫学または因果機械学習の実行
- 03再現可能なR・Quarto分析の準備
- 04実証分析の仕様変更の検証
インストール
コマンドをコピーしてターミナルで実行します。
$
npx skills add https:/インストール後、Claude Code で/full-empirical-analysis-skill-rを入力して読み込まれたか確認します。
同梱ファイル
SKILL.md と一緒にインストールされる参考資料やスクリプト。
README.mdreferences/03-descriptive-stats.md01-data-cleaning.md08-tables-plots.md06-robustness.md02-data-transformation.md04-statistical-tests.md05-modeling.md07-further-analysis.md詳細
完全な実証分析 — クラシックな R ワークフロー
このスキルは、応用経済学者がすべての実証論文で実行する標準的な8 ステップのパイプラインです。現代的な tidyverse と計量経済学の R エコシステムで記述されています。データには dplyr/tidyr/haven、パネル分析・IV・DID の主力パッケージには fixest、最新の DID には did/bacondecomp/HonestDiD、RD には rdrobust/rddensity、合成コントロールには Synth/gsynth/synthdid、マッチングには MatchIt/WeightIt/cobalt/ebal、ML による因果推論には grf/DoubleML、因果媒介分析には mediation、推定後分析には marginaleffects、出版用の表には modelsummary/kableExtra/gt、図には ggplot2/iplot/binsreg を使用します。
関連スキル：これは 00-StatsPAI_skill（Python DSL）、00.1-Full-empirical-analysis-skill（明示的な Python）、00.2-Full-empirical-analysis-skill_Stata（Stata .do）に対応する R 版です。4 つすべてが、それぞれのエコシステムで同じ 8 ステップを実装します。
基本方針
- Tidyverse + fixest による現代的な R の書き方。
lm(y ~ x + factor(unit) + factor(year))のような無理のある方法ではなく、feols(... | unit + year, cluster = ~unit)を使います。 - 再現可能なスクリプト / Quarto。 以下のすべての例は、コピーしてそのまま実行できます。パッケージのバージョン固定には
renv、説明文・コード・表・図をまとめるにはQuarto（.qmd）を使います。 - 8 ステップを重視。 R ユーザーは従来 Step 5 に過剰に注力しがちです。このスキルでは、Step 1〜4 と Step 6〜8 も中核として扱います。
- 豊富な出力。 各ステップで、少なくとも 1 つの表または図を生成します（tex/docx/png/pdf）。
- 段階的な情報開示。
SKILL.mdには各ステップの標準的な呼び出し方を記載し、references/にはバリエーション固有の詳細をまとめます。
SkillOpt 形式の実行ゲート
この詳細なプレイブックは、シードスキルとして使い、すべてを網羅的に適用するスクリプトとしては使わないでください。SkillOpt の原則に従い、R/Quarto に加える各変更を候補パッチとして扱います。再利用可能な定型コードにする前に、選択チェックを通過し、ホールドアウトチェックをクリアする必要があります。R スクリプトや Quarto ワークフローの作成・修正を始める前に、ユーザーの依頼をタスク固有の 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のいずれか）、推定対象、識別デザイン、主要なアウトカム／処置、R パッケージのファミリー、必要な成果物を記録します。 - 変更範囲の限定：一度に変更する判断事項は 1 つだけにします（サンプル規則、推定量、クラスタリング、エクスポート形式、またはロバストネスチェック）。検証に合格できる最小限のパッチを優先します。
- 選択用分割の規律：ユーザーが直面している問題、または依頼された成果物を選択用分割として扱います。別のアウトカム、サンプル期間、推定量ファミリー、またはエクスポート先を少なくとも 1 つ、ホールドアウトゲート用に確保します。
- ホールドアウトゲート：コードを実行する前にチェック項目を定義します。行数、
distinct()によるキーの一意性、処置のサポート、欠損値のしきい値、期待される表・図のファイル、および変更によって壊れてはならない主要対象外のロバストネスチェック／仕様を含めます。 - 却下バッファー：候補仕様がゲートに失敗した場合は、失敗内容、R/Quarto の差分、ゲートの出力を
analysis_log.mdに記録します。最後に合格した仕様に戻し、同じ未検証のパターンを再試行しないでください。 - スロー／メタ更新：タスクの最後に、作業の経緯に基づいて
accepted_patternsとrejected_patternsを書き留めます。1 回の実行で合格しただけでは、標準的なプロジェクトテンプレートの適用範囲を広げないでください。 - 検証後にのみ昇格：単発の修正を再利用可能なプロジェクトの定型コードにするのは、現在のデータと、少なくとも 1 つの別のアウトカム／サンプル／仕様で検証に合格した後に限ります。
3 つの分野モード（デフォルト＝AER 経済学、代替＝疫学と ML 因果推論）
上記のデフォルトのプレイブックは、AER 形式の応用計量経済学です。これは AEA の慣例に従い、推定式を明記し、識別仮定を示し、設計を複数手法で比較し、徹底的なロバスト性検証を行います。このスキルには、ほかの 2 つの主要な因果推論の伝統に対応する並行サブパイプラインも用意されています。どちらも同じステップ 1～4（クリーニング／構築／表 1／診断）とステップ 8（表／図）を再利用し、ステップ 5（推定量）とステップ 6／7 のパッケージのみを切り替えます。
| モード | 想定読者の慣例 | ステップ 5 の推定量スタック | 報告スタック | 参照先 |
|---|---|---|---|---|
| デフォルト — 応用経済学（AER / QJE / AEJ） | 「式と識別仮定を示し、設計を複数手法で比較する。統制変数を明示し、標準誤差はクラスター化する」 | DID / IV / RD / SCM / マッチング / fixest::feols HDFE | AER の書式に沿った複数列の modelsummary ＋ kableExtra / gt / flextable ＋ 8 セクションの論文構成 | ステップ 1 → 8（以下のプレイブック全体） |
| モード A — 疫学／公衆衛生 | 「STROBE / TRIPOD-AI、ターゲットトライアルプロトコル、二重ロバスト推定対象、絶対リスクと相対リスク、KM 生存分析」 | ターゲットトライアル・エミュレーション · IPTW（WeightIt / PSweight） · g-formula（gfoRmula） · TMLE（tmle / ltmle） · メンデルランダム化（MendelianRandomization / TwoSampleMR / MRPRESSO） · KM / Cox / AFT（survival / survminer / flexsurv） | 同じ modelsummary ＋ リスク差／ハザード比／E-value の行 | §A. 疫学パイプライン |
| モード B — ML 因果推論 | 「DML / メタラーナー / 因果フォレスト / DR ラーナー、CATE 分布、ポリシー値」 | DML（DoubleML） · S/T/X/R/DR-Learner（causalweight / grf） · GRF 因果フォレスト（grf::causal_forest） · BART/BCF（bartCause / bcf） · 行列補完（MCPanel） | modelsummary による ML 手法の比較 ＋ grf CATE プロット ＋ ポリシー値テーブル ＋ conformalInference PI | §B. ML 因果推論パイプライン |
デフォルト以外のモードの呼び出し方（Claude / エージェントはユーザーの表現からモードを判別します）：
| ユーザーの発言... | スキルが切り替えるモード |
|---|---|
| 「DID / IV / RD / イベントスタディを実行」「AER テーブル」「応用ミクロ」 | デフォルト（AER 経済学） — ステップ 1 → 8 |
| 「ターゲットトライアル・エミュレーション」「g-formula」「IPTW」「TMLE」「メンデルランダム化」「STROBE / TRIPOD」「公共衛生 / 疫学」「疫学パイプライン」「RWE 研究」「コホート研究」「症例対照研究」 | モード A（疫学） — §A |
| 「DML」「二重機械学習」「因果フォレスト」「メタラーナー」「CATE」「BCF」「policytree」「ポリシー学習」「コンフォーマル因果推論」「公平性監査」「ML 因果推論」「アップリフトモデリング」「因果機械学習」 | モード B（ML 因果推論） — §B |
| 「組み合わせ」（例：「DID を推定してから、異質性に対して ML CATE を適用」） | デフォルト ＋ モード B の順に実行 — すべての推定量が係数と標準誤差のペアを出力します。それらをすべて 1 つの modelsummary(...) に渡し、手法比較の列を作成します |
3 つのモードは、同じステップ 1～4 のクリーニング／表 1／診断の枠組み、同じステップ 8 のエクスポートスタック、同じ DAG 優先の識別ロジックを共有します。モードを切り替えると変わるのは、使用するステップ 5 の推定量ファミリーのみであり、論文全体の構成は変わりません。記述統計／表 1／バランス確認のみが必要な場合は、ステップ 3 の AER gtsummary::tbl_summary / modelsummary::datasummary_balance の呼び出しを 3 つすべてのモードで同じように使用できます。
デフォルト出力仕様 — 経済学実証論文
このスキルのデフォルトは、応用経済学論文の慣例です。ユーザーが単一の点推定値を明示的に求めない限り、実行のたびに以下の出版可能な完全な出力セットを生成します。これはステップ 8 の要件であり、必須であって、任意ではありません。
必須の表（常に作成）
| # | 表 | R のソース | 保存先 |
|---|---|---|---|
| T1 | 要約統計量とバランス（処置群と対照群、SMD／p 値を含む） | gtsummary::tbl_summary + add_p + add_difference（ステップ 3） | tables/table1_balance.xlsx + .docx + .tex |
| T2 ★ | 主な結果 — 複数列の回帰 M1→M6（段階的な統制変数 + FE） | fixest::feols × 6 仕様 → modelsummary（ステップ 5–6） | tables/table2_main.xlsx + .docx + .tex |
| T3 | メカニズム／アウトカムの階層 — 同じ処置、3 つ以上のアウトカムを横並びで表示 | ループで y ∈ {Y1, Y2, Y3, Y_main} に対して feols を実行 → modelsummary（ステップ 7） | tables/table3_mechanism.xlsx + .docx + .tex |
| T4 | 異質性 — サブグループ × 主な係数（性別、年齢、地域など） | サブグループの feols × linearHypothesis → modelsummary（ステップ 7） | tables/table4_heterogeneity.xlsx + .docx + .tex |
| T5 | 頑健性の検証 — 代替 SE／クラスター／サンプル／プラセボを1 つの表にまとめる | バリエーションごとに feols → modelsummary（ステップ 6） | tables/table5_robustness.xlsx + .docx + .tex |
★ 表 2 は、あらゆる経済学論文の中心となる表です。 生の相関（M1）から、完全に指定された設計（M6：2 方向 FE + 交互作用項を含む FE + クラスター頑健 SE）までを示す複数列の回帰表です。1 列にまとめないでください。 見出しとなる係数だけを報告しないでください。段階的な推移そのものが信頼性の根拠になります。M1→M6 が単調で安定していれば、その設計は識別に妥当性があると考えられます。FE を追加したときに崩れるなら、それ自体が結果です。
標準的な 6 列の順序：
- M1 生の二変量（
feols(y ~ treat, data)）- M2 + 人口統計学的変数（
+ age + edu）- M3 + 産業部門の統制変数（
+ tenure / firm_size）- M4 + 単位 FE（
| worker_id）- M5 + 2 方向 FE（
| worker_id + year）- M6 + 交互作用項を含む FE（
| worker_id + year + industry^year）、cluster = ~ worker_idを指定
必須の図（常に作成）
| # | 図 | R のソース | 保存先 |
|---|---|---|---|
| F1 | トレンド／動機づけ — 時間経過に伴う処置群と対照群の推移、政策開始時点の線を含む | dplyr のグループ平均 → ggplot + geom_line（ステップ 3） | figures/fig1_trend.png（300 dpi、PNG 形式で必ず出力）+ .pdf |
| F2 | 95% CI 付きイベントスタディ係数、基準期間は –1 | fixest::sunab()／did::ggdid／iplot（ステップ 5） | figures/fig2_event_study.png（300 dpi、PNG 形式で必ず出力）+ .pdf |
| F3 | 仕様 M1→M6 にわたる係数プロット | modelsummary::modelplot()（ステップ 8） | figures/fig3_coefplot.png（300 dpi、PNG 形式で必ず出力）+ .pdf |
| F4 | 頑健性／感度 — bacondecomp::bacon プロット、HonestDiD::createSensitivityPlot、または仕様曲線 | シナリオに応じて選択（ステップ 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 回だけ行う場合 — 該当するセルのみを作成します。ただし、標準の納品物は上記の一式であることを伝え、必要であれば一式を作成すると申し出てください。
- 設計上、図を作成できない場合（横断面データ → イベントスタディなし）— 理由を説明する
message()を出力して省略してください。何も説明せずに省略しないでください。 - 処置単位が N=1 の場合（
Synth／synthdid） — F1／F2 を SCM の軌跡 + プラセボ分布に置き換えてください。T1–T5 は引き続き作成します。
必須パッケージ
# Run once on a fresh R install:
install.packages(c(
# Data
"tidyverse", "haven", "readxl", "data.table", "janitor",
"naniar", "VIM", "mice", "validate",
# Description / tables
"gtsummary", "tableone", "modelsummary", "kableExtra", "gt",
"stargazer", "texreg", "flextable", "psych", "summarytools",
# Tests
"lmtest", "sandwich", "car", "tseries", "urca", "plm",
"clubSandwich", "fwildclusterboot",
# Modeling — workhorses
"fixest", # panel/IV/DID with HD FE — primary
"AER", # ivreg
"ivreg", # alternative IV
# Modern DID
"did", # Callaway–Sant'Anna
"didimputation", # Borusyak–Jaravel–Spiess
"fixest", # sunab() for Sun–Abraham
"synthdid", # Synthetic DID
"bacondecomp", "HonestDiD",
"DIDmultiplegtDYN", # de Chaisemartin–D'Haultfœuille
# RD
"rdrobust", "rddensity", "rdmulti",
# Synthetic control
"Synth", "gsynth", "tidysynth",
# Matching / weighting
"MatchIt", "WeightIt", "cobalt", "ebal",
# ML causal
"grf", "DoubleML",
# Mediation / SEM
"mediation", "lavaan",
# Robustness / inference
"robomit", # Oster delta
"ri2", "ritools", # randomization inference
"multcomp",
# Margins / post-estimation
"marginaleffects",
# Plotting
"ggplot2", "ggpubr", "cowplot", "patchwork",
"binsreg",
"ggdist", "ggrepel"
))
# fixest's iplot, esttex, etable are bundled.
8つのステップ — 標準パイプライン（AER論文の各セクションに対応）
┌──────────────────────────────────────────────────────────────────────┐
│ Step −1 Pre-Analysis Plan (PAP) pwr / WebPower / DeclareDesign │
│ Step 0 Sample log + data contract sample_log/stopifnot/jsonlite │
│ Step 1 Data import & cleaning read_csv/read_dta/janitor/naniar/mice│
│ Step 2 Variable construction mutate/across/winsorize/lag/group_by │
│ Step 2.5 Empirical strategy equation × ID assumption + pre-reg │
│ Step 3 Descriptive statistics gtsummary/datasummary_balance/cor_pmat│
│ Step 3.5 Identification graphics iplot/binsreg/rdplot/cobalt/Synth │
│ Step 4 Diagnostic tests shapiro/bptest/dwtest/vif/adf/kpss │
│ Step 5 Baseline modeling feols/ivreg/att_gt/synthdid/MatchIt │
│ Step 6 Robustness battery bacondecomp/HonestDiD/fwildclusterboot│
│ Step 7 Further analysis marginaleffects/mediation/grf │
│ Step 8 Tables & figures modelsummary/iplot/ggplot2/cowplot │
└──────────────────────────────────────────────────────────────────────┘
8つのステップは、応用経済学分野のAER／QJE／AEJ論文の標準的なセクションに対応しています。各ステップは論文の1つのセクションにあたり、論文にそのまま使える成果物をディスク上に出力します。
Paper section Step R moves
─────────────────────────── ───── ────────────────────────────────────────────────
Pre-Analysis Plan −1 pwr / WebPower / DeclareDesign + freeze pap.json
§1. Data 0 sample_log + 5-check stopifnot → JSON via jsonlite
§1. Data 1 haven::read_dta · janitor::clean_names · naniar/mice
§1. Data 2 mutate/across/Winsorize/lag/lead/diff · CPI deflate
§1.1 Descriptives (Table 1) 3 gtsummary::tbl_summary · datasummary_balance
§2. Empirical Strategy 2.5 write equation + ID assumption → strategy.md
§3. Identification graphics 3.5 fixest::iplot · binsreg · rdplot · cobalt::love.plot · Synth
§3.5 Diagnostics 4 bptest · dwtest · car::vif · urca::ur.df · phtest
§4. Main Results (Table 2) 5 fixest::feols progressive (m1...m6) · modelsummary
§5. Heterogeneity (Table 3) 7 feols(... + i(.):X) · marginaleffects::avg_slopes
§6. Mechanisms / Channels 7 mediation::mediate · lavaan · outcome ladder
§7. Robustness gauntlet 6 bacondecomp · HonestDiD · robomit · fwildclusterboot · ri2
§8. Replication package 8 modelsummary("...tex") · gt → docx · result.json
以下は、各ステップでの標準的な呼び出しです。すべての例で、労働経済学のパネルデータを使った同じ一連のシナリオを扱います。training（処置）がlog_wage（アウトカム）に与える影響を、共変量age、edu、tenure、パネルキーworker_id／firm_id／yearを使って分析します。変数名とパラメーター値は例示用です。
ステップに多くのバリエーションがある場合（段階的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 の両方の形式で保存する必要があります。
| § | 成果物 | R のプリミティブ | ファイル名 |
|---|---|---|---|
| §1 | 図1：未加工のトレンド／処置の導入 | df %>% group_by(year, treat) %>% summarise(mean(y)) %>% ggplot() | figures/fig1_trend.png(300dpi)+.pdf |
| §1 | 表1：要約統計量（全体／処置群／対照群＋Δ＋SMD） | gtsummary::tbl_summary · modelsummary::datasummary_balance | tables/table1_balance.xlsx/.docx/.tex |
| §3 | 図2：識別の図（イベントスタディ／第一段階／McCrary／RD 散布図／SCM 軌跡） | fixest::iplot(es) · binsreg · rdrobust::rdplot · rddensity · Synth::path.plot | figures/fig2_event_study.png(300dpi)+.pdf |
| §4 | 表2：主な結果 — 段階的な統制変数 M1→M6 | modelsummary(list("(1)"=m1,...,"(6)"=m6)) · fixest::etable | tables/table2_main.xlsx/.docx/.tex |
| §4 | 表2-bis：設計の比較（OLS／IV／DID／DML） | modelsummary(list("OLS"=ols, "2SLS"=iv, "CS-DID"=cs, "DML"=dml)) | tables/table2b_designs.xlsx/.docx/.tex |
| §4 | 図3：仕様別の係数プロット | modelplot(list(m1,...,m6), coef_map="training") | figures/fig3_coefplot.png(300dpi)+.pdf |
| §5 | 表3：サブグループ別の異質性 | modelsummary(g_full, g_male, g_fem, g_q1, ..., g_q4) | tables/table3_heterogeneity.xlsx/.docx/.tex |
| §5 | 図4：用量反応／CATE | marginaleffects::plot_predictions · grf::plot.causal_forest | figures/fig4_cate.png(300dpi)+.pdf |
| §6 | 表4：メカニズム／アウトカムの段階的分析 | アウトカムを順に feols で推定 → modelsummary | tables/table4_mechanism.xlsx/.docx/.tex |
| §7 | 表A1：頑健性検証の一覧（検証ごとに1列） | modelsummary(list(base, no99, balpan, dropearly, wfe, cl2way, logy, ihsy, psm, ebal)) | tables/tableA1_robustness.xlsx/.docx/.tex |
| §7 | 図5：仕様曲線 | specr::specr() + plot_specs（または purrr::pmap で独自に作成） | figures/fig5_spec_curve.png(300dpi)+.pdf |
| §7 | 図6：感度分析（HonestDiD／Oster／E-value） | HonestDiD::createSensitivityPlot · robomit::o_test · EValue | figures/fig6_sensitivity.png(300dpi)+.pdf |
| §8 | 再現用バンドル：すべての表を1つの文書にまとめる | modelsummary(..., output="docx") · gt::gtsave() · Quarto / Rmd | replication/paper_tables.xlsx/.docx/.tex |
上記の各 R 推定器（
fixest::feols／AER::ivreg／did::att_gt／grf::causal_forest／synthdid_estimate）は、modelsummary(...)／modelplot(...)／etable(...)にそのまま渡せる結果オブジェクトを返します。kable()から LaTeX を手作業で作成したり、flextableを直接使って Word に出力したりしないでください。modelsummary、etable、gtsummaryは、表の罫線、AER 形式の有意性を示す星印、適切な標準誤差ラベルを自動的に適用します。より詳しい出力方法については、references/08-tables-plots.mdを参照してください。
エクスポートの手引き — LaTeX／Word／Excelを1つのブロックで
重要なルール（必ず守ってください）：各表は、Excel（.xlsx）、Word（.docx）、LaTeX（.tex）の3形式すべてでエクスポートする必要があります。各画像は、PNG（300dpi以上）とPDFの両形式で保存する必要があります。
Rは、3つの言語の中で最も優れた出版用テーブルのエコシステムを備えています。対象範囲に応じて、次の3つのレベルから選びます。
| レベル | 使用する場面 | API | 主な引数 |
|---|---|---|---|
| 1. 単一の複数列テーブル | 段階的な列を持つTable 2／Table 3／Table A1を1つエクスポートする場合 | `modelsummary(list("(1)"=m1,...,"(N)"=mN), output="tables/tab.tex", stars=c(""=.1,""=.05,""=.01), gof_omit="BIC | AIC |
| 2. 複数パネル形式の論文（Table 2 + 3 + A1 + A2を1ファイルにまとめる） | メイン分析＋異質性分析＋頑健性分析＋プラセボ分析を1つの文書にまとめた論文の表ブロックを作成する場合 | modelsummaryをgt::gt_group()と組み合わせ、セクション見出し付きの1つの文書を作成するか、Quartoの.qmdで文章の間に複数のmodelsummary呼び出しをレンダリングする | gt_group(modelsummary(...), modelsummary(...)) · quarto render paper.qmd |
3. セッション全体のバンドル（Stataのcollect／PythonのStargazer + pylatexに相当） | 要約統計＋バランス検定＋複数の回帰表＋見出し＋文章を1つのファイルにまとめる再現用付録 | 現代のRネイティブな選択肢はQuartoです。master.qmdに文章と、modelsummary／gtsummary／ggplot2の出力を生成するチャンクを交互に配置します。quarto renderを1回実行すると、.pdf／.docx／.htmlが生成されます | YAMLフロントマターでformat: [pdf, docx, html]を設定し、3形式を出力します |
ジャーナル向けのスタイル — starsとSEラベルを適切に選びます。 AEAの慣例では、c("*"=.1, "**"=.05, "***"=.01)と、notes = "Cluster-robust standard errors in parentheses..."を使用します。master.Rの冒頭でラッパーを一度定義します。
# Top of master.R — journal house-style wrapper
# 输出三格式：.xlsx（编辑）、.docx（Word）、.tex（LaTeX）
aer_table <- function(models, output, headers = NULL, coef_map = NULL) {
base <- tools::file_path_sans_ext(output)
for (ext in c(".xlsx", ".docx", ".tex")) {
output_file <- paste0(base, ext)
fmt <- if (ext == ".xlsx") "html" else if (ext == ".docx") "docx" else "latex"
modelsummary(
models,
output = output_file,
stars = c("*" = 0.1, "**" = 0.05, "***" = 0.01),
gof_omit = "BIC|AIC|F|Log|Adj",
coef_map = coef_map,
notes = paste("Cluster-robust standard errors in parentheses.",
"* p<0.10, ** p<0.05, *** p<0.01."),
output_format = fmt
)
}
}
複数パネルの.docx／.xlsxとQuartoの手引き（論文の表を1ファイルにまとめる方法）については、references/08-tables-plots.mdを参照してください。
Step −1 — 事前分析計画（データを見る前；AEA RCT Registry形式）
データに触れる前に、（a）対象母集団、（b）研究デザイン、（c）予定サンプルサイズとα=0.05、β=0.20のもとでの**最小検出可能効果（MDE）**を記録します。結果をpap.jsonとして保存すれば、データを見る前に必要な検出力を確保したデザインだったことを査読者が確認できます。
library(pwr) # classical power calculations
library(WebPower) # cluster RCT, longitudinal, mixed designs
library(jsonlite)
# Two-sample MDE for a continuous outcome (Cohen's d framing)
pwr.t.test(d = 0.20, power = 0.80, sig.level = 0.05,
type = "two.sample", alternative = "two.sided")
# → required n per arm
# Solve for MDE given fixed n
pwr.t.test(n = 2000, power = 0.80, sig.level = 0.05,
type = "two.sample")$d
# → minimum detectable Cohen's d
# Cluster-randomized RCT — design effect
# Solve via WebPower::wp.crt2arm(...) for clusters / per-cluster size / power triangle
WebPower::wp.crt2arm(f = 0.20, J = NULL, n = 50, icc = 0.05, power = 0.80,
alpha = 0.05, alternative = "two.sided")
# → required clusters per arm
# DID power (Frison-Pocock / Bloom 1995): use WebPower::wp.kanova() or simulate
# RD power: simulate via DeclareDesign — see references/05-modeling.md §5.5
# Persist the protocol — referee will ask whether design was powered ex ante
pap <- list(
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>"
)
write_json(pap, "artifacts/pap.json", pretty = TRUE, auto_unbox = TRUE)
DAGを考慮した、より詳細な検出力分析（DAGを記述し、推定対象を宣言し、デザインをシミュレーションする）には、**DeclareDesign**を使用します。これは、EGAPの事前分析ワークフローに相当するRネイティブなツールです。
Step 1の前に、リポジトリでartifacts/pap.jsonをコミットします。AEA RCT Registry／OSFの事前登録ツールでは、これを分析計画の添付資料として使用できます。
Step 0 — サンプル構築ログと5項目のデータ契約
AERの§1「Data」セクションには、3つの役割があります。（a）データソースを説明すること、（b）すべてのサンプル制限（「脚注4」のサンプルログ）を記録すること、（c）パネル構造を確定することです。
0.1 サンプル構築ログ（脚注4）
library(tidyverse); library(jsonlite)
sample_log <- tibble::tibble(step = character(), n = integer())
df_raw <- read_dta("raw/panel.dta") %>% janitor::clean_names()
sample_log <- sample_log %>% add_row(step = "0. raw", n = nrow(df_raw))
df1 <- df_raw %>% drop_na(wage)
sample_log <- sample_log %>% add_row(step = "1. drop missing wage", n = nrow(df1))
df2 <- df1 %>% filter(between(age, 18, 65))
sample_log <- sample_log %>% add_row(step = "2. drop age outside 18-65", n = nrow(df2))
df3 <- df2 %>% filter(industry %in% c("manuf", "construction", "transport"))
sample_log <- sample_log %>% add_row(step = "3. keep target industries", n = nrow(df3))
df <- df3
print(sample_log)
write_json(sample_log, "artifacts/sample_construction.json", pretty = TRUE)
出力されたtibbleをそのまま論文の脚注4に貼り付けます。
0.2 5項目のデータ契約（実行可否の判定ゲート）
library(validate); library(assertr)
data_contract <- function(df, y, treatment, id = NULL, time = NULL, covariates = c()) {
keys <- c(y, treatment, id, time, covariates)
contract <- list(
n_obs = nrow(df), # 1. shape
dtypes = sapply(df[keys], function(x) class(x)[1]), # 2. dtypes
n_missing = sapply(df[keys], function(x) sum(is.na(x))), # 3. missingness
n_dupes_on_keys = if (!is.null(id) && !is.null(time))
sum(duplicated(df[, c(id, time)])) else 0, # 4. duplicates
panel_balanced = NULL,
cohort_sizes = NULL
)
if (!is.null(id) && !is.null(time)) {
bal <- df %>% count(.data[[id]])
contract$panel_balanced <- all(bal$n == max(bal$n)) # 5. balance
contract$n_dropped_by_balance <- sum(bal$n != max(bal$n))
if ("first_treat" %in% names(df)) {
contract$cohort_sizes <- df %>% distinct(.data[[id]], .keep_all = TRUE) %>%
count(first_treat) %>% deframe()
}
}
contract$y_range <- range(df[[y]], na.rm = TRUE)
contract$treatment_share <- mean(df[[treatment]], na.rm = TRUE)
# MCAR sniff test (Rubin) — if missing(y) is associated with covariates,
# listwise deletion biases the estimate. Use mice / IPW instead.
miss_y <- is.na(df[[y]])
contract$mcar_hint <- "likely MCAR (listwise OK)"
if (any(miss_y) && any(!miss_y)) {
for (cov in covariates) {
if (is.numeric(df[[cov]])) {
p <- t.test(df[[cov]][miss_y], df[[cov]][!miss_y])$p.value
if (p < 0.05) {
contract$mcar_hint <- sprintf("NOT MCAR (y-miss differs on %s, p=%.3f) → use mice / IPW",
cov, p)
break
}
}
}
}
contract
}
contract <- data_contract(df, y = "wage", treatment = "training",
id = "worker_id", time = "year",
covariates = c("age", "edu", "tenure"))
stopifnot(contract$n_dupes_on_keys == 0)
stopifnot(all(contract$n_missing == 0))
write_json(contract, "artifacts/data_contract.json",
pretty = TRUE, auto_unbox = TRUE)
stopifnotが1つでも実行された場合は、処理を止めて、まずdplyrで修正します。Rの推定器は後続処理でNA行を黙って除外します。このデータ契約は、「なぜTable 1からTable 2の間でNが12,000から9,800に減ったのか」という査読者からの質問に備える、最も低コストの対策です。
Step 1 — データのインポートとクリーニング
より詳しいパターンについては、references/01-data-cleaning.mdを参照してください。各種形式（haven／readxl／data.table::fread／arrow::read_parquet／DBI）、janitor::clean_names、欠損状況の可視化に使うnaniar、miceによるMCAR／MAR／MNARのトリアージ、validate／assertrによる検証、パネル構造の確認を扱っています。
library(tidyverse)
library(haven) # .dta / .sav / .sas7bdat
library(janitor) # clean_names()
library(naniar) # missing-data viz
library(skimr) # one-line dataset summary
# 1a. Load + first look
df <- read_dta("raw/panel.dta") %>%
clean_names() # standardize to snake_case
skim(df) # rich one-line-per-var summary
naniar::miss_var_summary(df)
naniar::vis_miss(df) # missingness heatmap
# 1b. Dtypes
df <- df %>%
mutate(
year = as.integer(year),
wage = as.numeric(wage),
gender = as.factor(gender),
date = as.Date(date)
)
# 1c. Missing values — decide PER VARIABLE
key_vars <- c("wage", "training", "worker_id", "year")
df <- df %>%
drop_na(all_of(key_vars))
cat("After dropping NA on keys:", nrow(df), "rows\n")
df <- df %>%
mutate(
tenure_missing = is.na(tenure),
tenure = if_else(is.na(tenure), median(tenure, na.rm = TRUE), tenure),
union = fct_explicit_na(as.factor(union), na_level = "unknown")
)
# 1d. Outliers — flag, don't drop yet
df <- df %>%
mutate(wage_z = scale(wage)[,1],
outlier_z4 = abs(wage_z) > 4)
cat("|z|>4 on wage:", sum(df$outlier_z4, na.rm = TRUE), "\n")
# 1e. Deduplicate panel key
stopifnot(nrow(df %>% distinct(worker_id, year)) == nrow(df))
# 1f. Merge with assertion
firm_chars <- read_dta("raw/firm_chars.dta")
n_before <- nrow(df)
df <- df %>%
left_join(firm_chars, by = "firm_id", relationship = "many-to-one")
stopifnot(nrow(df) == n_before) # no row inflation
# 1g. Panel structure
df %>% count(year) # per-year
df %>% count(worker_id) %>% summary() # per-unit
重要な原則：dplyrと明示的なstopifnot()アサーションを使用します。後続処理で行を黙って除外しないでください。
Step 2 — 変数の構築と変換
より詳しいパターンについては、references/02-data-transformation.md を参照してください。MASS::boxcox を使った対数変換／IHS変換／Box–Cox変換、dplyr を使ったグループごとのウィンザー化、scale() と bestNormalize、因子の処理、dplyr::lag を使ったラグ／リード、パネルデータのタイミングを扱います。
library(DescTools) # Winsorize()
df <- df %>%
mutate(
# 2a. Log / IHS
log_wage = log(pmax(wage, 1)),
ihs_assets = asinh(assets),
# 2b. Winsorize 1/99
wage_w1 = DescTools::Winsorize(wage, probs = c(0.01, 0.99), na.rm = TRUE),
# 2c. Standardize
age_std = as.numeric(scale(age)),
# 2d. Polynomial / interaction (or use formula syntax in fixest)
age_sq = age^2,
trt_x_edu = training * edu
) %>%
# 2e. Within-group winsorize
group_by(industry, year) %>%
mutate(wage_w1_iy = DescTools::Winsorize(wage, probs = c(0.01, 0.99),
na.rm = TRUE)) %>%
ungroup() %>%
# 2f. Panel operators (always arrange first to make lag deterministic)
arrange(worker_id, year) %>%
group_by(worker_id) %>%
mutate(
log_wage_l1 = lag(log_wage, 1),
log_wage_f1 = lead(log_wage, 1),
d_log_wage = log_wage - lag(log_wage, 1),
wage_mean_i = mean(log_wage, na.rm = TRUE),
log_wage_dm = log_wage - wage_mean_i
) %>%
ungroup() %>%
# 2g. Staggered-DID timing
group_by(worker_id) %>%
mutate(first_treat = ifelse(any(training == 1),
min(year[training == 1]), NA_real_)) %>%
ungroup() %>%
mutate(rel_time = year - first_treat,
never_treated = is.na(first_treat))
# 2h. CPI deflation
cpi <- read_csv("raw/cpi.csv")
df <- df %>%
left_join(cpi, by = "year") %>%
mutate(cpi_base = cpi[year == 2010][1],
wage_real = wage * cpi_base / cpi,
log_wage_real = log(pmax(wage_real, 1)))
Step 2.5 — 実証戦略（式と識別仮定を記述する）
ここがAER論文の核心です。コードを書く前に、式を明示的に記述し、識別仮定を述べてください。識別に関する曖昧な説明は、査読者が実証研究論文を却下する最も一般的な理由です。戦略は strategy.md として保存し、日付入りのバージョン管理された成果物にしてください。係数を見た後に事後的に正当化するために書くものではありません。
式 × 識別仮定 × R推定量（判断表）
| デザイン | 推定式 | 識別仮定 | R推定量 |
|---|---|---|---|
| 2×2 DID | Y_it = α_i + λ_t + β·D_it + X'γ + ε_it | X を条件とした平行トレンド | `feols(y ~ i(treated, post, ref=0) |
| イベントスタディ（CS / SA） | Y_it = α_i + λ_t + Σ_{e≠-1} β_e · 1{t-G_i = e} + ε_it | 予期効果なし＋グループ・時点別PT | `feols(y ~ sunab(G, t) |
| 2SLS | Y_i = α + β·D_i + X'γ + ε_i; D_i = π·Z_i + X'δ + u_i | 除外制約＋関連性＋単調性 | `feols(y ~ X |
| シャープRD | Y_i = α + β·1{X_i ≥ c} + f(X_i) + ε_i （局所多項式） | c における E[Y(0)|X] の連続性、操作なし | rdrobust::rdrobust(y, x, c=0) （＋ rddensity） |
| SCM | Ŷ_1t(0) = Σ_j ŵ_j Y_jt, τ_t = Y_1t − Ŷ_1t(0)（t≥T_0） | 介入前期間の適合＋補間の妥当性 | Synth::synth · gsynth::gsynth · synthdid::synthdid_estimate · tidysynth |
| 観測可能な共変量に基づく選択（マッチング/IPW/DML） | Y_i = m(X_i) + β·D_i + ε_i（Robinsonの部分化） | 非交絡性＋オーバーラップ | MatchIt::matchit ＋ lm · WeightIt · DoubleML::DoubleMLPLR · grf::causal_forest |
デザイン選択ガイド（ユーザーが判断に迷っている場合）
┌─ running var + cutoff ───────────────── RDD (rdrobust)
│
├─ exogenous instrument Z ─────────────── IV/2SLS (feols / AER::ivreg)
data + question ─┤
├─ pre/post × treat/control ─┬ 2 periods ── 2×2 DID (feols + i())
│ └ staggered ── CS / SA / BJS (att_gt / sunab / did_imputation)
│
├─ 1 treated unit + donor pool + long pre ── SCM (Synth / gsynth / synthdid)
│
├─ high-dim X, selection-on-observables ── ML causal (DoubleML / grf — see §B)
│
└─ none of the above ──────────────────── matching + sensitivity (MatchIt + EValue)
事前登録用 strategy.md テンプレート
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 `feols(y ~ sunab(G, t) | i + t, data)`
- Borusyak-Jaravel-Spiess via `didimputation::did_imputation`
- Synthetic DID via `synthdid::synthdid_estimate`
"
writeLines(strategy, "artifacts/strategy.md")
Step 5／Step 6 を実行する前に、リポジトリに artifacts/strategy.md をコミットしてください。このファイルの git log が分析計画そのものです。
Step 3 — 記述統計とTable 1
より詳しいパターンについては、references/03-descriptive-stats.md を参照してください。gtsummary::tbl_summary（現代のTable 1の標準）、SMD付きの modelsummary::datasummary_balance、tableone::CreateTableOne、corrplot ／ psych::corr.test による有意性付き相関行列、ggplot2 による分布プロットを扱います。
library(gtsummary)
library(modelsummary)
# 3a. Full-sample summary — one line, publication ready
df %>%
select(log_wage, age, edu, tenure, training) %>%
datasummary_skim()
# Or
df %>%
select(log_wage, age, edu, tenure, training) %>%
tbl_summary(
type = list(all_continuous() ~ "continuous2"),
statistic = all_continuous() ~ c("{N_nonmiss}", "{mean} ({sd})",
"{min} – {median} – {max}")
) %>%
bold_labels() %>%
as_kable_extra() %>%
kableExtra::save_kable("tables/table1_full.tex")
# 3b. Stratified Table 1 (treated vs control, with SMDs + p-values)
df %>%
select(log_wage, age, edu, tenure, female, training) %>%
tbl_summary(by = training, missing = "ifany") %>%
add_p() %>%
add_difference() %>%
add_n() %>%
modify_header(label = "**Variable**") %>%
bold_labels() %>%
as_gt() %>%
gt::gtsave("tables/table1_balance.html")
# Or via modelsummary (writes LaTeX/Word/HTML)
datasummary_balance(~ training,
data = df %>% select(training, age, edu, tenure, female),
output = "tables/table1_balance.tex")
# 3c. Correlation matrix with stars
library(corrplot); library(psych)
corr_obj <- corr.test(df %>% select(log_wage, age, edu, tenure, training),
method = "pearson")
corrplot(corr_obj$r, method = "color", type = "upper",
p.mat = corr_obj$p, sig.level = 0.05, insig = "blank",
addCoef.col = "black", number.cex = 0.7,
tl.col = "black", tl.srt = 45,
col = colorRampPalette(c("#B2182B","white","#2166AC"))(200))
# 3d. Distribution plots
library(ggplot2)
p1 <- ggplot(df, aes(log_wage, fill = factor(training))) +
geom_density(alpha = 0.5) +
scale_fill_manual(values = c("0" = "darkred", "1" = "navy"),
labels = c("Control", "Treated"), name = "") +
labs(x = "Log wage", y = "Density",
title = "Log-wage density by treatment") +
theme_classic()
p2 <- ggplot(df, aes(sample = log_wage)) +
stat_qq() + stat_qq_line() +
labs(title = "Normal Q-Q") + theme_classic()
cowplot::plot_grid(p1, p2, labels = "auto") %>%
ggsave("figures/distributions.pdf", plot = ., width = 10, height = 4)
# 3e. Time-trend (DID motivation)
df %>%
group_by(year, training) %>%
summarise(mean_log_wage = mean(log_wage, na.rm = TRUE), .groups = "drop") %>%
ggplot(aes(year, mean_log_wage, color = factor(training))) +
geom_line(linewidth = 1) + geom_point(size = 2) +
geom_vline(xintercept = policy_year, linetype = "dashed") +
scale_color_manual(values = c("0" = "darkred", "1" = "navy"),
labels = c("Control","Treated"), name = "") +
labs(x = "Year", y = "Mean log wage") + theme_classic()
ggsave("figures/trend_did.pdf", width = 7, height = 4)
Step 3.5 — 識別のグラフィック（「識別、グラフィカルな証拠」セクション）
AERの慣例では、識別を示す図を回帰表より先に提示します。 点推定値を信頼するよう求める前に、PT が成立している／第一段階が強い／RDで明確なジャンプがあることを示すグラフィカルな証拠を読者に見せてください。
3.5.1 イベントスタディ図＋数値による事前トレンド検定（DIDの識別）
介入前の係数がおよそ0（−1の参照期間は0に正規化）であることが、平行トレンドを示す視覚的な証拠になります。査読者が目視で判断しなくて済むよう、図と数値による事前トレンド検定を組み合わせてください。
library(fixest); library(ggplot2)
# (a) Sun-Abraham via fixest::sunab — the modern primary for staggered DID
es <- feols(log_wage ~ sunab(first_treat, year) | worker_id + year,
data = df, cluster = ~ worker_id)
# (b) Coefficient figure
iplot(es,
xlab = "Years relative to treatment",
ylab = "Coefficient (ATT, 95% CI)",
main = "Figure 2a. Event-study coefficients (95% CI; ref. e = -1)")
ggsave("figures/fig2a_event_study.pdf", width = 7, height = 4)
ggsave("figures/fig2a_event_study.png", width = 7, height = 4, dpi = 300)
# (c) Numerical pre-trends Wald test (joint zero on the leads)
pre_idx <- grep("year::-", names(coef(es)))[!grepl("ref", names(coef(es)))]
W <- wald(es, names(coef(es))[pre_idx])
cat(sprintf("Pre-trends Wald χ² = %.2f, p = %.3f\n", W$stat, W$p))
# (d) Bacon decomposition (Goodman-Bacon 2021) — TWFE diagnostic
library(bacondecomp)
bd <- bacon(log_wage ~ training, data = df,
id_var = "worker_id", time_var = "year")
ggplot(bd, aes(weight, estimate, color = type, shape = type)) +
geom_point(size = 2) +
labs(title = "Figure 2a-bis. Goodman-Bacon decomposition",
x = "Weight", y = "Estimate")
ggsave("figures/fig2a_bacon.pdf", width = 7, height = 4)
# (e) Callaway-Sant'Anna dynamic ATT (when att_gt is the main estimator)
library(did)
cs <- att_gt(yname = "log_wage", tname = "year", idname = "worker_id",
gname = "first_treat", data = df,
control_group = "nevertreated", est_method = "dr",
clustervars = "firm_id")
ggdid(aggte(cs, type = "dynamic")) +
labs(title = "Figure 2a-ter. Dynamic ATT (Callaway-Sant'Anna)")
ggsave("figures/fig2a_csdid.pdf", width = 7, height = 4)
3.5.2 第一段階のF統計量＋散布図（IVの識別）
経験則として、OLS型の推論では第一段階の F ≥ 10、AR同等の推論では F ≥ 23 です（Stock–Yogo／Lee 2022）。fixest::feols はF値を自動で報告します。AER::ivreg では summary(..., diagnostics = TRUE) が必要です。
iv <- feols(log_wage ~ age + edu | training ~ Z1 + Z2,
data = df, cluster = ~ firm_id)
summary(iv, stage = 1)
fitstat(iv, ~ ivf + ivwald + sargan + cd) # CD / KP / Sargan / first-stage F
# Binscatter for the first-stage scatter (residualized on age + edu)
library(binsreg)
binsreg(y = df$training, x = df$Z1, w = df[, c("age","edu")],
nbins = 20, polyreg = 2, ci = c(3, 3))
ggsave("figures/fig2b_first_stage.pdf", width = 7, height = 4)
3.5.3 RD：McCrary密度検定＋標準的なRDプロット
代表的なRD図は rdplot（各側の局所多項式フィットを伴うCCT形式のビン散布図）です。McCrary操作検定と組み合わせてください。
library(rdrobust); library(rddensity)
# (a) Canonical RD plot — binned means + local poly on each side
rdplot(y = df$outcome, x = df$running_var, c = 0,
p = 4, kernel = "triangular", binselect = "esmv",
title = "Figure 2c. RD plot")
ggsave("figures/fig2c_rdplot.pdf", width = 7, height = 4)
# (b) McCrary density (Cattaneo-Jansson-Ma 2018)
rdd <- rddensity(X = df$running_var, c = 0)
print(summary(rdd))
rdplotdensity(rdd, X = df$running_var,
title = "Figure 2c-bis. McCrary density (manipulation test)")
ggsave("figures/fig2c_mccrary.pdf", width = 7, height = 4)
3.5.4 マッチング：ラブプロット（マッチング前後の標準化差）
library(MatchIt); library(cobalt)
m.out <- matchit(training ~ age + edu + tenure + firm_size,
data = df, method = "nearest", ratio = 1)
love.plot(m.out, threshold = 0.10,
var.order = "unadjusted", abs = TRUE,
title = "Figure 2d. Love plot — |SMD| pre vs post matching")
ggsave("figures/fig2d_loveplot.pdf", width = 7, height = 4)
3.5.5 SCM: 合成対照の軌跡 + ギャッププロット
合成対照デザインでは、標準的な Figure 2 は処置時点を注記した処置群対合成群の時系列です。
library(tidysynth)
sc <- df %>%
synthetic_control(outcome = log_wage, unit = unit_id, time = year,
i_unit = "treated_unit_name", i_time = 2015) %>%
generate_predictor(time_window = 2010:2014,
mean_age = mean(age, na.rm = TRUE),
mean_edu = mean(edu, na.rm = TRUE)) %>%
generate_weights() %>% generate_control()
plot_trends(sc); ggsave("figures/fig2e_synth_trajectory.pdf", width = 7, height = 4)
plot_differences(sc); ggsave("figures/fig2e_synth_gap.pdf", width = 7, height = 4)
# Synthetic DID
library(synthdid)
sdid_setup <- panel.matrices(df, unit = "worker_id", time = "year",
outcome = "log_wage", treatment = "training")
sdid_fit <- synthdid_estimate(sdid_setup$Y, sdid_setup$N0, sdid_setup$T0)
plot(sdid_fit, control.name = "Synthetic DiD")
ggsave("figures/fig2e_sdid.pdf", width = 7, height = 4)
識別戦略固有のチェック（DID の PT、弱操作変数の F、RD の密度、マッチングの共通サポート）は Step-5 の推定量の内部でも自動実行されます — ここで数値を重複して計算しないでください。ただし、図は必ず作成してください。査読者はまず図に目を通します。
Step 4 — 診断的統計検定
より詳しいパターン: references/04-statistical-tests.md — あらゆる古典的検定を網羅しています。lmtest/sandwich/car/tseries/urca/plm。
library(lmtest)
library(sandwich)
library(car)
library(tseries)
library(urca)
# Fit baseline OLS for diagnostics
ols <- lm(log_wage ~ training + age + edu + tenure, data = df)
# 4a. Normality of residuals
shapiro.test(sample(residuals(ols), min(5000, length(residuals(ols)))))
tseries::jarque.bera.test(residuals(ols))
# 4b. Heteroskedasticity
bptest(ols) # Breusch-Pagan
bptest(ols, ~ I(fitted(ols)^2) + ., data = df) # White-style
# 4c. Autocorrelation (time series / panel)
dwtest(ols) # Durbin-Watson
bgtest(ols, order = 4) # Breusch-Godfrey
Box.test(residuals(ols), lag = 8, type = "Ljung-Box")
# Panel-specific
library(plm)
pdata <- pdata.frame(df, index = c("worker_id", "year"))
plm_fe <- plm(log_wage ~ training + age + edu, data = pdata, model = "within")
pbgtest(plm_fe) # Wooldridge serial correlation
pcdtest(plm_fe, test = "cd") # Pesaran cross-sectional dependence
# 4d. Multicollinearity
vif(ols) # VIFs
kappa(model.matrix(ols), exact = TRUE) # condition number
# 4e. Stationarity (time series — assumes a single y over time)
adf.test(df$log_wage, k = 4) # ADF
kpss.test(df$log_wage, null = "Level") # KPSS
# 4f. Hausman (FE vs RE)
plm_re <- plm(log_wage ~ training + age + edu, data = pdata, model = "random")
phtest(plm_fe, plm_re)
# 4g. Specification — RESET
resettest(ols, power = 2:3, type = "fitted")
判断表:
| 検定 | 帰無仮説 | 棄却された場合の対応 |
|---|---|---|
shapiro.test / jarque.bera.test | 残差が正規分布 | N が小さい場合はブートストラップ信頼区間 |
bptest | 等分散 | coeftest(ols, vcov = vcovHC(ols, "HC3")) による HC3、またはクラスター |
dwtest / bgtest | 自己相関なし | HAC 標準誤差（vcovHAC）、または単位によるクラスター |
pbgtest（パネル） | パネル自己相関なし | エンティティによるクラスター |
pcdtest | 横断面依存なし | Driscoll–Kraay（vcovDC） |
vif > 10 | — | 除外 / 統合 |
| ADF が棄却 + KPSS が棄却しない | 定常 | 水準 |
| ADF が棄却しない | 単位根 | 一階差分 |
phtest | RE が一致 | FE を使用 |
Step 5 — ベースラインの実証モデリング（セクション 4: 主結果）
より詳しいパターン: references/05-modeling.md — あらゆる推定量を網羅しています。fixest が主力です。
ここは実証論文で最も密度の高いセクションです。現代の AER §4 には通常、2〜3 個の重回帰表と 1 個の係数プロットが含まれます:
- Table 2（メイン）: 段階的なコントロール変数、4〜6 列 — 以下の Pattern A
- Table 2-bis（デザインの競合比較）: OLS / IV / DID / DML の下での同じ係数 — Pattern B
- Table 2-ter（複数アウトカム）: 同じ処置、複数のアウトカムを並べて表示 — Pattern C
- Figure 3（coefplot）: 各仕様にわたる β̂ と 95% 信頼区間の視覚的要約
推定量のルーティング（これを暗記してください — 間違えると気づかないうちにナンセンスな結果が出力されます）:
- FE なし / 単一の低基数 FE →
feols(y ~ X, data, cluster = ~i)- 高次元 FE →
feols(y ~ X | fe1 + fe2, data, cluster = ~i)- 二方向クラスター →
feols(..., cluster = ~ firm_id + year)- 2SLS / IV →
feols(y ~ X | D ~ Z, data, cluster = ~ firm_id)（診断にはAER::ivregも可）- DID / イベントスタディ →
feols(y ~ sunab(G, t) | i + t, data)（SA）·did::att_gt（CS）·didimputation::did_imputation（BJS）
識別戦略による選択:
Cross-section, selection on observables → feols | MatchIt + lm | WeightIt
Panel + policy shock + parallel trends → feols / did::att_gt / sunab / didimputation / synthdid
Exogenous instrument → feols(... | endog ~ z) | AER::ivreg
Discontinuity → rdrobust + rddensity + rdmc
N=1 treated, long panel → Synth / gsynth / synthdid
Selection on observables + heterogeneity → WeightIt + cobalt; grf::causal_forest
Binary outcome → feglm or glm(family=binomial)
Count outcome → fepois
標準的な呼び出し（以下の 8 つのパターン A〜H は AER の表のクックブックです — modelsummary(...) と fixest::etable(...) が 2 つの主力であり、Stata の outreg2/esttab や Python の pf.etable/Stargazer に相当します）。
5.A Pattern A — 段階的コントロール変数（標準的な Table 2）
列をまたいで β̂ が安定している ⇒ 観測可能変数による選択が推定値を駆動しているという懸念が少なくなります（Oster 2019 の選択安定性の論理。Step 6 で定量化）。
library(fixest); library(modelsummary)
m1 <- feols(log_wage ~ training, data = df, cluster = ~ firm_id)
m2 <- feols(log_wage ~ training + age + edu, data = df, cluster = ~ firm_id)
m3 <- feols(log_wage ~ training + age + edu + tenure + firm_size, data = df, cluster = ~ firm_id)
m4 <- feols(log_wage ~ training + age + edu + tenure + firm_size | industry + year,
data = df, cluster = ~ firm_id)
m5 <- feols(log_wage ~ training + age + edu + tenure + firm_size | worker_id + year,
data = df, cluster = ~ firm_id)
m6 <- feols(log_wage ~ training + age + edu + tenure + firm_size | worker_id + year + industry^year,
data = df, cluster = ~ firm_id)
modelsummary(
list("(1) Baseline" = m1,
"(2) +Demog" = m2,
"(3) +Labor-mkt" = m3,
"(4) Ind×Yr FE" = m4,
"(5) Worker FE" = m5,
"(6) Worker FE+Ind×Yr" = m6),
output = "tables/table2_main.tex",
stars = c("*" = 0.1, "**" = 0.05, "***" = 0.01),
gof_omit = "BIC|AIC|F|Log|Adj",
coef_map = c("training" = "Job training",
"age" = "Age", "edu" = "Education",
"tenure" = "Tenure", "firm_size" = "Firm size"),
notes = c("Cluster-robust SE in parentheses, clustered at firm_id.",
"* p<0.10, ** p<0.05, *** p<0.01.")
)
modelsummary(list("(1)"=m1,"(2)"=m2,"(3)"=m3,"(4)"=m4,"(5)"=m5,"(6)"=m6),
output = "tables/table2_main.docx")
AER の慣習: すべてのコントロール変数（および切片）を表示する。
keep =もcoef_omit =も渡さないでください。そうすればすべてのパラメータが可視化されます。coef_map = c("training" = "Training")（単一マッピング）を使うのは、特定の関心のある係数のみの表が意図されている場合（交互作用形式の異質性、IV 第1段階の三点セット）だけにしてください。coef_omit = "Intercept"を使うのは、論文の見た目のために定数を抑制したい場合だけにしてください。
5.B Pattern B — デザインの競合比較（Table 2-bis）
複数の識別戦略の下で同じ関心のある係数を示します。これは まさに AER の信頼性を示す手法です。それぞれ異なる識別仮定を置く各デザインにわたる収束的な証拠です。
library(fixest); library(AER); library(did); library(MatchIt); library(WeightIt)
ols <- feols(log_wage ~ training + age + edu + tenure | industry + year,
data = df, cluster = ~ firm_id)
iv <- feols(log_wage ~ age + edu + tenure | training ~ Z1 + Z2,
data = df, cluster = ~ firm_id)
cs <- att_gt(yname = "log_wage", tname = "year", idname = "worker_id",
gname = "first_treat", data = df,
control_group = "nevertreated", est_method = "dr",
clustervars = "firm_id")
psm <- matchit(training ~ age + edu + tenure, data = df,
method = "nearest", ratio = 1)
psm_lm <- lm(log_wage ~ training + age + edu + tenure,
data = match.data(psm), weights = weights)
ebal <- weightit(training ~ age + edu + tenure, data = df, method = "ebal")
ebal_lm <- lm(log_wage ~ training + age + edu + tenure,
data = df, weights = ebal$weights)
modelsummary(
list("(1) OLS+FE" = ols,
"(2) 2SLS" = iv,
"(3) CS-DID" = aggte(cs, type = "simple"),
"(4) PSM" = psm_lm,
"(5) Entropy bal." = ebal_lm),
output = "tables/table2b_designs.tex",
stars = c("*" = 0.1, "**" = 0.05, "***" = 0.01),
coef_map = c("training" = "Job training (β̂)"),
gof_omit = "BIC|AIC|F|Log|Adj",
notes = "Convergent evidence: same β̂ under five identification strategies."
)
5.C Pattern C — 複数アウトカム表（同じ X、複数の Y）
ys <- c("log_wage", "weeks_employed", "left_firm", "promoted")
multi_y <- lapply(ys, function(y)
feols(as.formula(paste(y, "~ training + age + edu + tenure | industry + year")),
data = df, cluster = ~ firm_id))
names(multi_y) <- ys
modelsummary(multi_y,
output = "tables/table2c_multi_outcome.tex",
stars = c("*" = 0.1, "**" = 0.05, "***" = 0.01),
coef_map = c("training" = "Training"),
notes = "Each column is a separate regression on the labelled outcome.")
5.D Pattern D — Panel A / Panel B を積み重ねた表
同じモデルファミリー、2 つの期間（短期 / 長期）または 2 つのサンプル。gt::gt_group() を使って 2 つの modelsummary ブロックをパネルヘッダー付きで積み重ねます。
library(gt)
panelA <- list(
"(1) Industry FE" = feols(wage_t1 ~ training + X | industry + year, data = df, cluster = ~ firm_id),
"(2) Worker FE" = feols(wage_t1 ~ training + X | worker_id + year, data = df, cluster = ~ firm_id))
panelB <- list(
"(1) Industry FE" = feols(wage_t5 ~ training + X | industry + year, data = df, cluster = ~ firm_id),
"(2) Worker FE" = feols(wage_t5 ~ training + X | worker_id + year, data = df, cluster = ~ firm_id))
ms_A <- modelsummary(panelA, output = "gt") %>%
tab_header(title = "Panel A. Short-run (1 year)")
ms_B <- modelsummary(panelB, output = "gt") %>%
tab_header(title = "Panel B. Long-run (5 years)")
gt_group(ms_A, ms_B) %>%
gtsave("tables/table2d_horizons.tex")
gt_group(ms_A, ms_B) %>%
gtsave("tables/table2d_horizons.docx")
5.E パターン E — IV 報告トリプレット（第一段階／縮約形／2SLS）
教科書的な AER の IV 表では、第一段階、縮約形、2SLSを3列で提示し、読者が Wald 比 = 縮約形 / 第一段階であることを確認できるようにします。
fs <- feols(training ~ Z + age + edu | industry + year, data = df, cluster = ~ firm_id)
rf <- feols(log_wage ~ Z + age + edu | industry + year, data = df, cluster = ~ firm_id)
iv2 <- feols(log_wage ~ age + edu | training ~ Z, data = df, cluster = ~ firm_id)
modelsummary(
list("(1) First stage" = fs,
"(2) Reduced form" = rf,
"(3) 2SLS" = iv2),
output = "tables/table2e_iv_triplet.tex",
stars = c("*" = 0.1, "**" = 0.05, "***" = 0.01),
coef_map = c("Z" = "Instrument Z", "training" = "Training (endog.)"),
gof_map = list(list(raw = "ivf", clean = "First-stage F", fmt = 2)),
notes = "Wald ratio: $\\hat\\beta_{2SLS} = \\hat\\beta_{RF} / \\hat\\pi_{FS}$."
)
IV トリプレットは意図的に焦点を絞っています： 読者が Wald 比を一目で確認できるよう、Z + 内生変数のみを表示します。査読者から係数一覧をすべて表示するよう求められた場合に限り、
coef_map=を削除します。
5.F パターン F — did::att_gt / synthdid / grf::causal_forest を使った因果推論オーケストレーターのメイン分析
DID / SCM / マッチング / フォレストのメイン分析では、最新の R 推定器が、推定値と自動生成されるプラセボ検定／事前トレンド／オーバーラップ診断をまとめて返します。自動 tidy 関数を介して modelsummary にパイプします。
# CS-DID with full diagnostics
cs <- att_gt(yname = "log_wage", tname = "year", idname = "worker_id",
gname = "first_treat", data = df,
control_group = "nevertreated", est_method = "dr",
clustervars = "firm_id")
print(aggte(cs, type = "group")) # ATT(g) summary
print(aggte(cs, type = "dynamic", min_e = -4, max_e = 4)) # event-study aggregation
# Synthetic DID
library(synthdid)
sdid_setup <- panel.matrices(df, unit="worker_id", time="year",
outcome="log_wage", treatment="training")
sdid_fit <- synthdid_estimate(sdid_setup$Y, sdid_setup$N0, sdid_setup$T0)
print(summary(sdid_fit))
# Causal forest with overlap + variable importance
library(grf)
cf <- causal_forest(X = as.matrix(df[, c("age","edu","tenure","firm_size")]),
Y = df$log_wage, W = df$training, num.trees = 4000)
average_treatment_effect(cf, target.sample = "treated")
test_calibration(cf)
variable_importance(cf)
5.G パターン G — サブグループ別 modelsummary（表3、ステップ7を参照）
サブグループごとに1列を設けます。詳細なコードは §ステップ7 — 異質性を参照してください。
5.H パターン H — 頑健性検証マスター（表A1、ステップ6を参照）
すべての頑健性仕様をベースラインの隣に並べます。詳細なコードは §ステップ6を参照してください。
標準的な推定器コマンド（基盤となるプリミティブ）
library(fixest)
# 5a. OLS with cluster-robust SEs — feols is the modern primary
ols <- feols(log_wage ~ training + age + edu + tenure,
data = df, cluster = ~ firm_id)
summary(ols)
# 5b. Two-way FE — single line
fe <- feols(log_wage ~ training + age + edu + tenure | worker_id + year,
data = df, cluster = ~ worker_id)
# Multi-way clustering
fe_mw <- feols(log_wage ~ training | worker_id + year,
data = df, cluster = ~ worker_id + firm_id)
# High-dim interaction FE
fe_hd <- feols(log_wage ~ training | worker_id + industry^year,
data = df, cluster = ~ firm_id)
# 5c. 2×2 DID
did22 <- feols(log_wage ~ i(treated, post, ref = 0) + age + edu,
data = df, cluster = ~ worker_id)
# Or with absorbed FE:
did22 <- feols(log_wage ~ i(treated, post, ref = 0) | worker_id + year,
data = df, cluster = ~ worker_id)
# 5d. Event study — base period at -1
es <- feols(log_wage ~ i(rel_time, ref = -1) | worker_id + year,
data = df %>% filter(!is.na(first_treat)),
cluster = ~ worker_id)
iplot(es,
xlab = "Years relative to treatment",
main = "Event study")
# 5e. Staggered DID — modern estimators (see references/05-modeling.md §5.4)
library(did)
cs <- att_gt(yname = "log_wage", tname = "year", idname = "worker_id",
gname = "first_treat", data = df,
control_group = "nevertreated",
est_method = "dr",
clustervars = "firm_id")
ggdid(cs) # event-study plot
# Sun & Abraham via fixest::sunab
sa <- feols(log_wage ~ sunab(first_treat, year) | worker_id + year,
data = df, cluster = ~ worker_id)
iplot(sa, sub.title = "Sun-Abraham (2021)")
# Borusyak–Jaravel–Spiess (didimputation)
library(didimputation)
bjs <- did_imputation(data = df, yname = "log_wage", gname = "first_treat",
tname = "year", idname = "worker_id",
horizon = 0:5, pretrends = -5:-1,
cluster_var = "worker_id")
# Synthetic DID
library(synthdid)
sdid_setup <- synthdid::panel.matrices(df, unit = "worker_id", time = "year",
outcome = "log_wage", treatment = "training")
sdid_fit <- synthdid_estimate(sdid_setup$Y, sdid_setup$N0, sdid_setup$T0)
# 5f. IV / 2SLS
iv <- feols(log_wage ~ age + edu | training ~ draft_lottery + z2,
data = df, cluster = ~ firm_id)
summary(iv, stage = 1)
fitstat(iv, ~ ivf + ivwald + sargan) # first-stage F + Wald + overid
# Or via AER:
library(AER)
iv_aer <- ivreg(log_wage ~ training + age + edu |
draft_lottery + z2 + age + edu, data = df)
summary(iv_aer, vcov. = sandwich, diagnostics = TRUE)
# 5g. Sharp RD
library(rdrobust); library(rddensity)
rd <- rdrobust(y = df$outcome, x = df$running_var, c = 0,
kernel = "triangular", bwselect = "mserd")
summary(rd)
rdplot(y = df$outcome, x = df$running_var, c = 0)
rddensity(X = df$running_var, c = 0) # manipulation test
# 5h. Binary outcome
logit <- feglm(employed ~ training + age + edu | firm_id + year,
data = df, family = binomial(link = "logit"),
cluster = ~ firm_id)
library(marginaleffects)
avg_slopes(logit, variables = "training") # AME
# 5i. Count w/ HD FE
pois <- fepois(citations ~ training + age | firm_id + year,
data = df, cluster = ~ firm_id)
ステップ6 — 頑健性検証一式
より詳細なパターン：references/06-robustness.md — M1～M6 には modelsummary、clubSandwich/fwildclusterboot、bacondecomp/HonestDiD/robomit、ランダム化推論には ri2 を使用します。
library(modelsummary)
library(fixest)
# 6a. Progressive specs (M1 → M6)
m1 <- feols(log_wage ~ training, data = df, cluster = ~ firm_id)
m2 <- feols(log_wage ~ training + age + edu, data = df, cluster = ~ firm_id)
m3 <- feols(log_wage ~ training + age + edu + tenure | worker_id,
data = df, cluster = ~ worker_id)
m4 <- feols(log_wage ~ training + age + edu + tenure | worker_id + year,
data = df, cluster = ~ worker_id)
m5 <- feols(log_wage ~ training + age + edu + tenure | worker_id + year + region,
data = df, cluster = ~ worker_id)
m6 <- feols(log_wage ~ training + age + edu + tenure | worker_id + year + industry^year,
data = df, cluster = ~ worker_id)
modelsummary(list("(1)" = m1, "(2)" = m2, "(3)" = m3,
"(4)" = m4, "(5)" = m5, "(6)" = m6),
stars = c('*' = .1, '**' = .05, '***' = .01),
gof_omit = "BIC|AIC|F|Log",
coef_map = c("training" = "Training",
"age" = "Age", "edu" = "Education", "tenure" = "Tenure"),
output = "tables/table_main.tex")
# 6b. Alternative cluster levels
for (cl in c("worker_id", "firm_id", "industry", "state")) {
fit <- feols(log_wage ~ training | worker_id + year, data = df,
cluster = as.formula(paste0("~", cl)))
cat(cl, ": b=", coef(fit)["training"], " se=", se(fit)["training"], "\n")
}
# 6c. Wild cluster bootstrap (when few clusters)
library(fwildclusterboot)
boot <- boottest(m4, param = "training", clustid = "state",
B = 9999, seed = 42)
summary(boot)
# 6d. Subsample splits
splits <- list(
"Female=0" = df %>% filter(female == 0),
"Female=1" = df %>% filter(female == 1),
"Young (<40)" = df %>% filter(age < 40),
"Old (>=40)" = df %>% filter(age >= 40)
)
sub_fits <- imap(splits, ~ feols(log_wage ~ training | worker_id + year,
data = .x, cluster = ~ worker_id))
modelsummary(sub_fits, stars = TRUE)
# 6e. Placebo — fake timing
df_placebo <- df %>%
mutate(fake_first = first_treat - 3,
fake_post = year >= fake_first) %>%
filter(year < first_treat)
feols(log_wage ~ fake_post | worker_id + year,
data = df_placebo, cluster = ~ worker_id)
# 6f. Randomization inference
library(ri2)
ri_out <- conduct_ri(formula = log_wage ~ training + age + edu,
declaration = randomizr::declare_ra(N = nrow(df),
prob = mean(df$training)),
assignment = "training",
sharp_hypothesis = 0,
data = df,
sims = 1000)
summary(ri_out); plot(ri_out)
# 6g. TWFE bias diagnosis
library(bacondecomp)
bacon_out <- bacon(log_wage ~ training,
data = df, id_var = "worker_id", time_var = "year")
ggplot(bacon_out, aes(weight, estimate, color = type)) + geom_point()
ggsave("figures/bacon.pdf")
# 6h. Parallel-trends sensitivity
library(HonestDiD)
honest_out <- createSensitivityResults(betahat = es$coefficients,
sigma = vcov(es),
numPrePeriods = 5, numPostPeriods = 5,
Mbarvec = seq(0, 0.5, by = 0.05))
createSensitivityPlot(honest_out, originalResults = honest_out$mainResult)
ggsave("figures/honestdid.pdf")
# 6i. Oster (2019) δ*
library(robomit)
o_test(y = "log_wage", x = "training",
con = "age + edu + tenure | worker_id + year",
id = "worker_id", time = "year",
data = df, R2max = 1.3 * fitstat(m6, "r2"), beta = 0)
# ============================================================
# 6j. Pattern H — Robustness master table (Table A1, one column per check)
# ============================================================
library(modelsummary); library(MatchIt); library(WeightIt)
base <- feols(log_wage ~ training + age + edu + tenure | industry + year,
data = df, cluster = ~ firm_id)
no99 <- feols(log_wage ~ training + age + edu + tenure | industry + year,
data = df %>% filter(wage < quantile(wage, 0.99, na.rm = TRUE)),
cluster = ~ firm_id)
balpan <- feols(log_wage ~ training + age + edu + tenure | industry + year,
data = df %>% group_by(worker_id) %>%
filter(n_distinct(year) == max(n_distinct(year))) %>% ungroup(),
cluster = ~ firm_id)
dropearly <- feols(log_wage ~ training + age + edu + tenure | industry + year,
data = df %>% filter(first_treat > 2008), cluster = ~ firm_id)
wfe <- feols(log_wage ~ training + age + edu + tenure | worker_id + year,
data = df, cluster = ~ firm_id)
cl2way <- feols(log_wage ~ training + age + edu + tenure | industry + year,
data = df, cluster = ~ firm_id + year)
logy <- feols(log(wage + 1) ~ training + age + edu + tenure | industry + year,
data = df, cluster = ~ firm_id)
ihsy <- feols(asinh(wage) ~ training + age + edu + tenure | industry + year,
data = df, cluster = ~ firm_id)
m_psm <- matchit(training ~ age + edu + tenure + firm_size, data = df, method = "nearest")
psm_lm <- lm(log_wage ~ training + age + edu + tenure, data = match.data(m_psm), weights = weights)
ebal_w <- weightit(training ~ age + edu + tenure + firm_size, data = df, method = "ebal")
ebal_lm <- lm(log_wage ~ training + age + edu + tenure, data = df, weights = ebal_w$weights)
modelsummary(
list("(1) Baseline" = base,
"(2) Drop top 1%" = no99,
"(3) Balanced" = balpan,
"(4) Drop early" = dropearly,
"(5) Worker FE" = wfe,
"(6) 2-way cluster" = cl2way,
"(7) log Y" = logy,
"(8) IHS Y" = ihsy,
"(9) PSM" = psm_lm,
"(10) Entropy bal." = ebal_lm),
output = "tables/tableA1_robustness.tex",
stars = c("*" = 0.1, "**" = 0.05, "***" = 0.01),
coef_map = c("training" = "Training (β̂)"),
gof_omit = "BIC|AIC|F|Log|Adj",
notes = "Each column is one robustness check. β̂ on training is the focal coefficient."
)
# ============================================================
# 6k. Specification curve (Simonsohn-Simmons-Nelson 2020) via `specr`
# ============================================================
library(specr); library(ggplot2)
specs <- setup(data = df,
y = c("log_wage", "ihs_wage"),
x = "training",
model = c("feols"),
controls = c("age", "edu", "tenure", "firm_size"),
subsets = list(industry = c("manuf", "construction", "transport")))
results <- specr(specs)
plot(results, choices = c("x", "y", "controls", "subsets"))
ggsave("figures/fig5_spec_curve.pdf", width = 10, height = 6)
ggsave("figures/fig5_spec_curve.png", width = 10, height = 6, dpi = 300)
# Hand-rolled alternative when `specr` doesn't fit (with custom FE / SE):
# spec_grid <- expand.grid(controls = list(c("age"), c("age","edu"), c("age","edu","tenure")),
# ytrans = c("log_wage", "ihs_wage"),
# sample = c("all", "manuf", "no99"),
# cluster = c("firm_id", "firm_id+year"))
# Loop, run feols, collect b/se, ggplot::geom_pointrange.
# ============================================================
# 6l. Sensitivity dashboard — HonestDiD + Oster + E-value
# ============================================================
# (a) HonestDiD — Rambachan-Roth (2023): bound on β̂ under bounded PT violation
library(HonestDiD)
es_pre <- coef(es)[grep("year::-", names(coef(es)))]
es_post <- coef(es)[grep("year::[0-9]", names(coef(es)))]
honest_out <- createSensitivityResults(betahat = c(es_pre, es_post),
sigma = vcov(es)[c(names(es_pre), names(es_post)),
c(names(es_pre), names(es_post))],
numPrePeriods = length(es_pre),
numPostPeriods = length(es_post),
Mbarvec = seq(0, 0.5, by = 0.05))
createSensitivityPlot(honest_out, originalResults = honest_out$mainResult)
ggsave("figures/fig6_honestdid.pdf", width = 7, height = 4)
# (b) Oster δ — `robomit::o_test` (already shown in 6i)
# (c) E-value (VanderWeele-Ding 2017) — for risk-ratio outcomes
library(EValue)
evalue(RR(1.45), lo = 1.10, hi = 1.91)
# → reports the minimum strength of unmeasured confounding to nullify the result.
ステップ7 — 追加分析
より詳細なパターン：references/07-further-analysis.md — 事後推定には marginaleffects が主力です。Imai の媒介分析には mediation::mediate、SEM には lavaan、CATE には grf::causal_forest を使用します。
library(marginaleffects)
library(fixest)
# 7a. Heterogeneity via interaction
het <- feols(log_wage ~ i(female, training, ref = 0) + age + edu | worker_id + year,
data = df, cluster = ~ worker_id)
summary(het)
iplot(het) # visualize interaction
# Continuous moderator + marginsplot
het_c <- feols(log_wage ~ training * tenure + age + edu | worker_id + year,
data = df, cluster = ~ worker_id)
plot_slopes(het_c, variables = "training",
condition = list(tenure = seq(0, 20, by = 2))) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(x = "Tenure", y = "Marginal effect of training")
ggsave("figures/het_tenure.pdf", width = 6, height = 4)
# 7b. Triple difference
ddd <- feols(log_wage ~ treated * post * high_exposure | worker_id + year,
data = df, cluster = ~ firm_id)
# 7c. Outcome ladder
out_ladder <- list()
for (y in c("hours_worked", "productivity", "log_wage")) {
out_ladder[[y]] <- feols(as.formula(paste(y, "~ training | worker_id + year")),
data = df, cluster = ~ worker_id)
}
modelsummary(out_ladder, stars = TRUE,
coef_map = c("training" = "Training"),
output = "tables/outcome_ladder.tex")
# 7d. Mediation — Imai et al. (2010)
library(mediation)
med_M <- lm(hours_worked ~ training + age + edu, data = df)
med_Y <- lm(log_wage ~ training + hours_worked + age + edu, data = df)
med <- mediate(med_M, med_Y, treat = "training", mediator = "hours_worked",
boot = TRUE, sims = 1000)
summary(med); plot(med)
# Sensitivity to unobserved M-Y confounding
medsens <- medsens(med, rho.by = 0.05, effect.type = "indirect")
plot(medsens)
# 7e. CATE via causal forest
library(grf)
cf <- causal_forest(X = as.matrix(df %>% select(age, edu, tenure, firm_size)),
Y = df$log_wage, W = df$training,
num.trees = 2000, min.node.size = 5)
df$tau_hat <- predict(cf)$predictions
variable_importance(cf)
average_treatment_effect(cf, target.sample = "all")
# Plot CATE by a moderator
ggplot(df, aes(tenure, tau_hat)) +
geom_smooth(method = "loess", se = TRUE) +
labs(x = "Tenure", y = "Estimated CATE")
ggsave("figures/cate_tenure.pdf")
# 7f. Dose-response — splines
library(splines)
dr <- feols(log_wage ~ ns(training_hours, df = 4) + age + edu | worker_id + year,
data = df, cluster = ~ worker_id)
plot_predictions(dr, condition = "training_hours")
ステップ8 — 論文掲載用の表と図
このステップは必須です — すべての分析実行で、このスキルの冒頭にある デフォルト出力仕様 で定義された必須表5つ（T1～T5）と必須図4つ（F1～F4）をすべて作成します。「回帰分析はすでに実行した」という理由でステップ8を省略しないでください。表や図のない係数は、応用経済学における結果の伝え方ではありません。
より詳細なパターン：references/08-tables-plots.md — modelsummary は最新の標準ツールです（1回の呼び出しで LaTeX／Word／HTML／Excel を出力）。LaTeX の詳細なスタイル設定には kableExtra、HTML／Word には gt、図には ggplot2 + iplot + ggpubr + cowplot + binsreg を使用します。
library(modelsummary)
library(kableExtra)
library(gt)
library(fixest)
library(ggplot2)
# ============================================================
# 8a. ★ TABLE 2 — Main results, multi-column regression M1→M6
# (the centerpiece of every economics paper)
# ============================================================
modelsummary(
list("(1) Raw" = m1,
"(2) +Demog" = m2,
"(3) +Tenure" = m3,
"(4) +Unit FE" = m4,
"(5) +2-way FE" = m5,
"(6) +Ind×Yr FE" = m6),
stars = c('*' = .1, '**' = .05, '***' = .01),
coef_map = c("training" = "Training",
"age" = "Age", "edu" = "Education", "tenure" = "Tenure"),
gof_map = list(
list("raw" = "nobs", "clean" = "N", "fmt" = 0),
list("raw" = "r.squared", "clean" = "R²", "fmt" = 3),
list("raw" = "adj.r.squared","clean" = "Adj. R²", "fmt" = 3)
),
notes = "Cluster-robust SE at worker_id in parentheses. * p<0.10, ** p<0.05, *** p<0.01.",
output = "tables/table2_main.tex"
)
modelsummary(list("(1)"=m1, "(2)"=m2, "(3)"=m3, "(4)"=m4, "(5)"=m5, "(6)"=m6),
stars = TRUE, output = "tables/table2_main.docx")
# ============================================================
# 8b. TABLE 1 — Summary statistics & balance
# ============================================================
library(gtsummary)
tbl1 <- df %>%
select(log_wage, age, edu, tenure, female, training) %>%
tbl_summary(by = training, missing = "ifany",
statistic = all_continuous() ~ "{mean} ({sd})") %>%
add_p() %>% add_difference() %>% add_n() %>% bold_labels()
tbl1 %>% as_kable_extra(format = "latex", booktabs = TRUE) %>%
kableExtra::save_kable("tables/table1_balance.tex")
tbl1 %>% as_flex_table() %>%
flextable::save_as_docx(path = "tables/table1_balance.docx")
# ============================================================
# 8c. TABLE 3 — Mechanism / outcome ladder (3+ outcomes)
# ============================================================
ladder <- list()
for (y in c("hours_worked", "productivity", "log_wage")) {
ladder[[y]] <- feols(as.formula(paste(y, "~ training + age + edu + tenure | worker_id + year")),
data = df, cluster = ~ worker_id)
}
modelsummary(ladder,
stars = c('*' = .1, '**' = .05, '***' = .01),
coef_map = c("training" = "Training"),
notes = "Each column is a separate regression on the labelled outcome. Cluster-robust SE at worker_id.",
output = "tables/table3_mechanism.tex")
# ============================================================
# 8d. TABLE 4 — Heterogeneity (subgroup × main coef)
# ============================================================
het_specs <- list(
"All" = df,
"Female=0" = df %>% filter(female == 0),
"Female=1" = df %>% filter(female == 1),
"Age<40" = df %>% filter(age < 40),
"Age≥40" = df %>% filter(age >= 40),
"Manuf." = df %>% filter(industry == "manufacturing")
)
het_models <- imap(het_specs,
~ feols(log_wage ~ training + age + edu + tenure | worker_id + year,
data = .x, cluster = ~ worker_id))
modelsummary(het_models,
stars = c('*' = .1, '**' = .05, '***' = .01),
coef_map = c("training" = "Training"),
notes = "Cluster-robust SE at worker_id. Wald p-values for cross-subgroup equality should accompany this table — see references/07.",
output = "tables/table4_heterogeneity.tex")
# ============================================================
# 8e. TABLE 5 — Robustness battery (alt SE / cluster / sample / placebo)
# ============================================================
rob <- list(
"Baseline" = feols(log_wage ~ training | worker_id + year, data = df,
cluster = ~ worker_id),
"Cluster=Firm" = feols(log_wage ~ training | worker_id + year, data = df,
cluster = ~ firm_id),
"2-way Cluster" = feols(log_wage ~ training | worker_id + year, data = df,
cluster = ~ worker_id + firm_id),
"Winsor 1/99" = feols(log_wage ~ training | worker_id + year,
data = df %>% mutate(log_wage = DescTools::Winsorize(log_wage,
probs = c(.01,.99),
na.rm = TRUE)),
cluster = ~ worker_id),
"Drop Manuf." = feols(log_wage ~ training | worker_id + year,
data = df %>% filter(industry != "manufacturing"),
cluster = ~ worker_id),
"Placebo (-3)" = feols(log_wage ~ fake_post | worker_id + year,
data = df %>% filter(year < first_treat),
cluster = ~ worker_id)
)
modelsummary(rob,
stars = c('*' = .1, '**' = .05, '***' = .01),
output = "tables/table5_robustness.tex")
# ============================================================
# 8f. ★ FIGURE 3 — Coefficient plot across M1→M6
# ============================================================
modelplot(list("(1)"=m1, "(2)"=m2, "(3)"=m3, "(4)"=m4, "(5)"=m5, "(6)"=m6),
coef_map = c("training" = "Training"),
conf_level = 0.95) +
geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) +
labs(x = "Coefficient on training (95% CI)", y = "Specification",
title = "Effect of training across specifications") +
theme_classic(base_size = 11)
ggsave("figures/fig3_coefplot.pdf", width = 6, height = 4)
ggsave("figures/fig3_coefplot.png", width = 6, height = 4, dpi = 300)
# ============================================================
# 8g. FIGURE 2 — Event-study plot (dynamic DID, base period = -1)
# ============================================================
pdf("figures/fig2_event_study.pdf", width = 7, height = 4)
iplot(es,
xlab = "Years relative to treatment",
ylab = "Coefficient (ATT, 95% CI)",
main = "Event study: dynamic effect of training",
ref.line = -0.5)
dev.off()
png("figures/fig2_event_study.png", width = 2100, height = 1200, res = 300)
iplot(es,
xlab = "Years relative to treatment",
ylab = "Coefficient (ATT, 95% CI)",
main = "Event study: dynamic effect of training",
ref.line = -0.5)
dev.off()
# ============================================================
# 8h. FIGURE 4 — Sensitivity / robustness curve
# (HonestDiD / spec curve / forest of robustness battery)
# ============================================================
# HonestDiD example (after the event study with stored b/V):
library(HonestDiD)
honest_out <- createSensitivityResults(betahat = es$coefficients,
sigma = vcov(es),
numPrePeriods = 5, numPostPeriods = 5,
Mbarvec = seq(0, 0.5, by = 0.05))
sens_plot <- createSensitivityPlot(honest_out, originalResults = honest_out$mainResult)
ggsave("figures/fig4_sensitivity.pdf", plot = sens_plot, width = 7, height = 4)
ggsave("figures/fig4_sensitivity.png", plot = sens_plot, width = 7, height = 4, dpi = 300)
# Alternative — robustness forest plot:
# rob_summary <- imap_dfr(rob, ~ tibble(
# group = .y,
# est = coef(.x)[1],
# se = se(.x)[1]
# ))
# ggplot(rob_summary, aes(est, fct_rev(factor(group)))) +
# geom_point(size = 3, color = "navy") +
# geom_errorbarh(aes(xmin = est - 1.96*se, xmax = est + 1.96*se),
# height = 0.2, color = "navy") +
# geom_vline(xintercept = 0, linetype = "dashed") +
# labs(x = "Coefficient on training (95% CI)", y = NULL,
# title = "Robustness forest plot")
# ggsave("figures/fig4_sensitivity.pdf", width = 7, height = 4)
# ============================================================
# 8i. FIGURE 1 — Trend / motivation (treated vs control over time)
# ============================================================
df %>%
group_by(year, training) %>%
summarise(mean_log_wage = mean(log_wage, na.rm = TRUE), .groups = "drop") %>%
ggplot(aes(year, mean_log_wage, color = factor(training))) +
geom_line(linewidth = 1) + geom_point(size = 2) +
geom_vline(xintercept = policy_year, linetype = "dashed", color = "gray40") +
scale_color_manual(values = c("0" = "darkred", "1" = "navy"),
labels = c("Control", "Treated"), name = "") +
labs(x = "Year", y = "Mean log wage",
title = "Treated vs control trend") +
theme_classic(base_size = 11) +
theme(legend.position = "bottom")
ggsave("figures/fig1_trend.pdf", width = 7, height = 4)
ggsave("figures/fig1_trend.png", width = 7, height = 4, dpi = 300)
# ============================================================
# 8j. Auxiliary plots (optional — produce when relevant)
# ============================================================
library(binsreg)
binsreg(y = df$log_wage, x = df$tenure, w = df %>% select(age, edu, female))
ggsave("figures/figA_binscatter.pdf", width = 6, height = 4)
# RD plot (only when running_var exists)
# rdplot(y = df$outcome, x = df$running_var, c = 0,
# title = "RD plot", x.label = "Running variable", y.label = "Outcome")
# ============================================================
# 8k. Multi-panel combined (optional, for slides / appendix)
# ============================================================
library(cowplot)
# plot_grid(p_trend, p_event, p_coef, p_sens, ncol = 2, labels = "AUTO") %>%
# ggsave("figures/combined.pdf", plot = ., width = 10, height = 8)
# ============================================================
# 8l. Theme — set once at top of script for consistency
# ============================================================
theme_set(theme_classic(base_size = 11) +
theme(legend.position = "bottom",
plot.title = element_text(face = "bold")))
納品物チェックリスト（実行完了を宣言する前に確認してください）：
[ ] 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 8m)
8m. 再現性スタンプ
ジャーナルの再現性検証担当者（または将来の共著者）が、主要な数値を再現するために必要な唯一の成果物です。R のバージョン、シード、データセットのハッシュ、ベースラインの係数 + CI、プロトコル／契約への参照先を保存します：
library(jsonlite); library(digest)
# Get baseline result (assumes `base` is the headline feols object)
b_hat <- coef(base)["training"]
se_b <- se(base)["training"]
ci <- c(b_hat - 1.96 * se_b, b_hat + 1.96 * se_b)
stamp <- list(
R_version = R.version.string,
fixest_version = as.character(packageVersion("fixest")),
modelsummary_version = as.character(packageVersion("modelsummary")),
seed = 42,
dataset_sha256 = substr(digest::digest(df, algo = "sha256"), 1, 16),
n_obs = base$nobs,
estimand = "ATT",
estimator = "fixest::feols",
estimate = unname(b_hat),
se_cluster = unname(se_b),
ci95 = unname(ci),
pre_registration = "artifacts/strategy.md",
data_contract = "artifacts/data_contract.json",
sample_log = "artifacts/sample_construction.json",
paper_bundle = "tables/table2_main.tex"
)
write_json(stamp, "artifacts/result.json", pretty = TRUE, auto_unbox = TRUE)
論文の PDF とともに artifacts/result.json をコミットします。査読者が Rscript master.R を実行し、この JSON をビット単位で完全に再現できる必要があります。
§A — 疫学／公衆衛生モード
ユーザーの表現からモード A（ターゲット・トライアル・エミュレーション／IPTW／TMLE／MR／STROBE／疫学／公衆衛生／RWE／コホート）が示される場合も、8つのステップは引き続き適用されます。ただしステップ5では、OLS と固定効果の構成を、二重ロバスト法 + 生存時間分析 + MR のトリプレットに置き換え、納品物は STROBE／TRIPOD-AI の規約に従います。ステップ1～4（クリーニング、変数作成、表1、診断）とステップ8（表／図のエクスポート）は、デフォルトモードと同じです。
パッケージ構成（デフォルトの構成に追加してインストールします）：
install.packages(c(
"WeightIt", "PSweight", "cobalt", # IPTW / propensity weighting + balance
"gfoRmula", # parametric g-formula (time-varying)
"tmle", "ltmle", # TMLE / longitudinal TMLE
"survival", "survminer", "flexsurv", # KM / Cox / AFT / RMST
"MendelianRandomization", "TwoSampleMR", # IVW, Egger, weighted-median MR
"MRPRESSO", # outlier-robust MR
"EValue" # E-value sensitivity (VanderWeele)
))
A.0 コホートの構築 + ターゲット・トライアル・プロトコル
データに触れる前にプロトコルを作成します。protocol.yml として保存し、論文内で引用します。
# protocol.yml — target-trial emulation skeleton
# eligibility: age 40-75, no_prior_event, ascertained_at t0
# treatment: A=1 statin initiation; A=0 no initiation
# assignment: emulated random at t0 via IPTW on baseline covariates
# outcome: incident MI within 5 years
# estimand: ITT ATE on risk difference + hazard ratio
library(dplyr)
cohort <- df |>
filter(age >= 40, age <= 75, prior_MI == 0) |>
mutate(
t0 = coalesce(statin_initiation_date, enrollment_date),
event_5y = as.integer((MI_date - t0) <= 365 * 5 & !is.na(MI_date)),
time_at_risk = pmin(as.numeric(censor_date - t0), 365 * 5)
)
A.1 曝露別の表1（デフォルトのステップ3と同じ）
ステップ3と同じ gtsummary::tbl_summary を使用し、by = A とします。未測定交絡に関する E 値はフッターに記載します。
library(gtsummary)
cohort |>
select(A, age, edu, smoke, bmi, ldl, sbp) |>
tbl_summary(by = A, missing = "ifany") |>
add_difference() |>
add_p() |>
bold_labels()
A.2 DAG + 傾向スコアのオーバーラップ（ポジティビティの確認）
library(WeightIt); library(cobalt)
# Estimate PS + IPTW weights
w_out <- weightit(A ~ age + edu + smoke + bmi + ldl + sbp,
data = cohort, method = "glm", estimand = "ATE")
# Overlap density (positivity)
bal.plot(w_out, var.name = "prop.score", which = "both")
ggsave("figures/figA2_ps_overlap.pdf")
# Love plot (SMDs before vs after IPTW)
love.plot(w_out, threshold = 0.1, abs = TRUE)
ggsave("figures/figA2_love.pdf")
A.3 IPTW + g-formula + TMLE の二重ロバストな三つの推定法（ステップ 5 の差し替え）
疫学の「AER Table 2」です。3 列の表で、各列に {IPTW-MSM, g-formula, TMLE} のいずれかを配置し、二重ロバストな推定結果が一致することを確認できます。
# IPTW marginal structural model
library(survey)
des <- svydesign(ids = ~1, data = cohort, weights = w_out$weights)
msm <- svyglm(event_5y ~ A, design = des, family = quasibinomial())
RD_iptw <- coef(msm)["A"]; CI_iptw <- confint(msm)["A", ]
# g-formula (parametric, time-fixed)
library(gfoRmula)
gf <- gformula_binary_eof(
obs_data = cohort,
id = "subject_id", time_name = "t", outcome_name = "event_5y",
covnames = c("age","edu","smoke","bmi","ldl","sbp"),
intvars = list("A"), interventions = list(list(c(static, 1)), list(c(static, 0))),
ref_int = 1, time_points = 1, basecovs = c("age","edu","smoke","bmi","ldl","sbp")
)
# TMLE (doubly robust)
library(tmle)
fit_tmle <- tmle(
Y = cohort$event_5y, A = cohort$A,
W = cohort[, c("age","edu","smoke","bmi","ldl","sbp")],
family = "binomial",
Q.SL.library = c("SL.glm","SL.glmnet","SL.ranger"),
g.SL.library = c("SL.glm","SL.glmnet","SL.ranger")
)
RD_tmle <- fit_tmle$estimates$ATE$psi
CI_tmle <- fit_tmle$estimates$ATE$CI
# Stack the triplet into one paper table
library(modelsummary)
tableA3 <- tibble::tribble(
~Estimator, ~RD, ~`95% CI`,
"IPTW-MSM", RD_iptw, sprintf("[%.3f, %.3f]", CI_iptw[1], CI_iptw[2]),
"g-formula", gf$result[2,"mean"] - gf$result[1,"mean"], "—",
"TMLE", RD_tmle, sprintf("[%.3f, %.3f]", CI_tmle[1], CI_tmle[2])
)
modelsummary::datasummary_df(tableA3, output = "tables/tableA3_dr_triplet.tex")
A.4 生存アウトカム — KM / Cox / AFT / RMST
library(survival); library(survminer); library(flexsurv)
# KM by treatment
fit_km <- survfit(Surv(time_at_risk, event_5y) ~ A, data = cohort)
ggsurvplot(fit_km, conf.int = TRUE, pval = TRUE, risk.table = TRUE)
ggsave("figures/figA4_km.pdf")
# Cox HR (covariate-adjusted)
fit_cox <- coxph(Surv(time_at_risk, event_5y) ~ A + age + edu + smoke + bmi + ldl + sbp,
data = cohort, weights = w_out$weights)
HR <- exp(coef(fit_cox)["A"]); HR_CI <- exp(confint(fit_cox)["A", ])
# AFT (Weibull) for time-ratio interpretation
fit_aft <- flexsurvreg(Surv(time_at_risk, event_5y) ~ A + age + edu + smoke + bmi + ldl + sbp,
data = cohort, dist = "weibull")
# RMST contrast at t = 5 years
library(survRM2)
rmst <- rmst2(cohort$time_at_risk, cohort$event_5y, cohort$A, tau = 365 * 5)
A.5 メンデルランダム化（IVW / Egger / 重み付き中央値の三つの推定法）
library(MendelianRandomization)
mri <- mr_input(bx = BX, bxse = BXSE, by = BY, byse = BYSE,
exposure = "Statin use", outcome = "MI")
ivw <- mr_ivw(mri)
egger <- mr_egger(mri) # pleiotropy intercept test
wmedian<- mr_median(mri, weighting = "weighted")
# Or harmonized two-sample workflow
# library(TwoSampleMR); harmonised <- harmonise_data(exposure_dat, outcome_dat)
# res <- mr(harmonised, method_list = c("mr_ivw", "mr_egger_regression", "mr_weighted_median"))
# Sensitivity to outliers
library(MRPRESSO)
mr_presso(BetaOutcome = "by", BetaExposure = "bx", SdOutcome = "byse", SdExposure = "bxse",
OUTLIERtest = TRUE, DISTORTIONtest = TRUE, data = data.frame(bx, by, bxse, byse), NbDistribution = 1000)
A.6 頑健性 — E-value / 境界 / 主層別化
library(EValue)
ev <- evalue(RR(1.45), lo = 1.10, hi = 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 / メタラーナー / 因果フォレスト / BCF / CATE / ポリシー学習 / 共形因果推論 / 公平性 / 因果機械学習）が示される場合、パイプラインではデフォルトモードのステップ 1〜4 とステップ 8 を維持し、ステップ 5 を ML 推定器スタックに差し替え、ステップ 7 とステップ 8 の間に CATE 分布とポリシー価値のレイヤーを追加します。
パッケージ構成（デフォルトスタックに追加してインストール）：
install.packages(c(
"DoubleML", "mlr3", "mlr3learners", # DML + ML nuisance learners
"grf", # causal forest, GRF, instrumental forest
"causalweight", # IPW / DR / sensitivity for CATE
"bartCause", "bcf", # BART / Bayesian causal forest
"policytree", # honest policy trees
"conformalInference", # conformal prediction (general)
# cfcausal — install via devtools::install_github("lihualei71/cfcausal")
"fairmodels", # fairness audit
"pcalg", "bnlearn" # causal discovery (PC / GES / Bayesian net)
))
B.0 訓練／ホールドアウト分割 + ナイサンス学習器スタック
library(mlr3); library(mlr3learners); library(DoubleML)
set.seed(42)
idx <- sample(seq_len(nrow(df)), size = 0.7 * nrow(df))
train <- df[idx, ]; holdout <- df[-idx, ]
# Standard nuisance pair: outcome regression Q(X,A) and propensity g(A|X)
ml_g <- lrn("regr.ranger", num.trees = 500, mtry = 5) # outcome
ml_m <- lrn("classif.ranger", num.trees = 500, mtry = 5) # propensity
B.1 DAG / 推定対象量の宣言（必要に応じて LLM を活用）
library(pcalg)
# PC algorithm — constraint-based DAG discovery
suffStat <- list(C = cor(df[, c("A","Y","X1","X2","X3","X4")]), n = nrow(df))
pc.fit <- pc(suffStat, indepTest = gaussCItest,
alpha = 0.01, labels = c("A","Y","X1","X2","X3","X4"))
plot(pc.fit, main = "PC-recovered DAG")
# OR: bnlearn for hill-climbing GES
# library(bnlearn); hc.fit <- hc(df[, c("A","Y","X1","X2","X3","X4")]); plot(hc.fit)
B.2 推定器スタック — DML · メタラーナー · 因果フォレスト · BCF（ステップ 5 の差し替え）
ML 因果推論の「AER Table 2」です。同じ (Y, A, X) データに対し、各列に異なる推定器ファミリーを配置した比較表です。読者は、DML、T-learner、因果フォレスト、BCF のすべてで ATE の推定結果が一致するか（または一致しないか）を確認したいと考えています。
# DML — partially linear or interactive regression model
dml_data <- DoubleMLData$new(train, y_col = "Y", d_cols = "A",
x_cols = c("X1","X2","X3","X4"))
dml_plr <- DoubleMLPLR$new(dml_data, ml_g = ml_g, ml_m = ml_m, n_folds = 5)
dml_plr$fit()
ate_dml <- dml_plr$coef; ci_dml <- dml_plr$confint()
# Causal forest (GRF) — non-parametric CATE
library(grf)
cf <- causal_forest(X = as.matrix(train[, c("X1","X2","X3","X4")]),
Y = train$Y, W = train$A, num.trees = 2000)
ate_cf <- average_treatment_effect(cf, target.sample = "all")
cate_cf <- predict(cf, newdata = as.matrix(holdout[, c("X1","X2","X3","X4")]))$predictions
# T-learner / DR-learner (use causalweight or hand-rolled with grf::*)
library(causalweight)
dr <- treatDML(y = train$Y, d = train$A, x = as.matrix(train[, c("X1","X2","X3","X4")]),
MLmethod = "lasso")$effect
ate_DR <- mean(dr)
# Bayesian Causal Forest — separate prognostic + treatment functions
library(bcf)
bcf_fit <- bcf(y = train$Y, z = train$A,
x_control = as.matrix(train[, c("X1","X2","X3","X4")]),
x_moderate = as.matrix(train[, c("X1","X2","X3","X4")]),
pihat = predict(glm(A ~ ., data = train[, c("A","X1","X2","X3","X4")], family = binomial), type = "response"),
nburn = 1000, nsim = 1000)
ate_bcf <- mean(bcf_fit$tau)
# Stack the horse-race
library(modelsummary)
tableB2 <- tibble::tribble(
~Estimator, ~ATE,
"DML (PLR)", ate_dml[1],
"Causal Forest", ate_cf[1],
"DR-learner", ate_DR,
"Bayesian Causal Forest", ate_bcf
)
modelsummary::datasummary_df(tableB2, fmt = 4, output = "tables/tableB2_ml_horserace.tex")
B.3 CATE 分布 + サブグループ CATE プロット（ステップ 7 の拡張）
library(ggplot2)
# CATE histogram
data.frame(cate = cate_cf) |>
ggplot(aes(x = cate)) +
geom_histogram(bins = 30, fill = "grey70", colour = "black") +
geom_vline(xintercept = 0, lty = 2) +
labs(x = "CATE", y = "Count")
ggsave("figures/figB3_cate_hist.pdf")
# CATE by quartile of a covariate
holdout |>
mutate(cate = cate_cf, age_q = ntile(X1, 4)) |>
group_by(age_q) |>
summarise(mean_cate = mean(cate)) |>
ggplot(aes(age_q, mean_cate)) + geom_col() + labs(y = "Mean CATE")
ggsave("figures/figB3_cate_by_age_q.pdf")
B.4 ポリシー学習 + オフポリシー評価
library(policytree)
# Honest discrete policy tree on doubly-robust scores from causal forest
dr_scores <- double_robust_scores(cf)
ptree <- policy_tree(X = as.matrix(train[, c("X1","X2","X3","X4")]),
Gamma = dr_scores, depth = 3)
print(ptree) # human-readable tree of "treat if X1<a and X2>b"
plot(ptree)
ggsave("figures/figB4_policy_tree.pdf")
# Off-policy evaluation — DR policy value on holdout
holdout_X <- as.matrix(holdout[, c("X1","X2","X3","X4")])
pred_pol <- predict(ptree, holdout_X)
DR_holdout <- double_robust_scores(cf, newdata = holdout_X)
policy_value_DR <- mean(DR_holdout[cbind(seq_len(nrow(DR_holdout)), pred_pol)])
cat(sprintf("DR policy value (holdout): %.3f\n", policy_value_DR))
B.5 不確実性（共形因果推論）+ 公平性 + 感度分析
# Conformal prediction interval around CATE (split conformal via cfcausal)
# devtools::install_github("lihualei71/cfcausal")
library(cfcausal)
ci90 <- conformalIte(X = as.matrix(train[, c("X1","X2","X3","X4")]),
Y = train$Y, T = train$A,
alpha = 0.1,
algo = "nest",
type = "CQR",
X.test = as.matrix(holdout[, c("X1","X2","X3","X4")]))
# Fairness audit — disparate impact / equalised odds
library(fairmodels)
fobject <- fairness_check(model_treated = predict(ptree, holdout_X),
data = holdout, protected = holdout$sensitive_attr,
privileged = "majority")
plot(fobject)
ggsave("figures/figB5_fairness.pdf")
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 (grf::test_calibration)
[ ] 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 | データの読み込み | haven / readr / readxl / data.table::fread | Parquet には arrow |
| 1 | 列名の整理 | janitor::clean_names | 手作業 |
| 1 | 欠損値 | naniar / mice | Hmisc |
| 2 | ウィンザー化 | DescTools::Winsorize | 手作業の pmin/pmax |
| 2 | パネルデータのラグ | dplyr::lag（arrange+group_byと併用） | data.table::shift |
| 3 | Table 1 | gtsummary / modelsummary::datasummary_balance | tableone |
| 3 | 相関 | psych::corr.test + corrplot | Hmisc::rcorr |
| 4 | 不均一分散／自己相関 | lmtest::bptest / dwtest / bgtest | car |
| 4 | パネル検定 | plm::pbgtest / pcdtest / phtest | — |
| 4 | 定常性 | tseries::adf.test / tseries::kpss.test | urca |
| 5 | OLS／パネル固定効果 | fixest::feols | lfe::felm（旧式） |
| 5 | IV | fixest::feols(| ~ ) | AER::ivreg / ivreg::ivreg |
| 5 | DID — 2×2 | i(treated, post)を使ったfeols | — |
| 5 | DID — CS | did::att_gt | — |
| 5 | DID — SA | fixest::sunab | — |
| 5 | DID — BJS | didimputation::did_imputation | — |
| 5 | DID — SDID | synthdid | — |
| 5 | RD | rdrobust / rddensity / rdmulti | — |
| 5 | SC | Synth / gsynth / tidysynth | — |
| 5 | PSM | MatchIt::matchit | — |
| 5 | IPW | WeightIt::weightit | — |
| 5 | エントロピーバランシング | ebal | — |
| 5 | DML | DoubleML | — |
| 5 | CATE（因果フォレスト） | grf::causal_forest | — |
| 5 | 媒介分析 | mediation::mediate | lavaan |
| 6 | クラスターワイルドブートストラップ | fwildclusterboot::boottest | clubSandwich |
| 6 | ランダム化推論 | ri2::conduct_ri | 手作業の boot |
| 6 | 多重検定 | multcomp / Romano-Wolf法を手作業で実装 | — |
| 6 | TWFEの診断 | bacondecomp::bacon | — |
| 6 | 平行トレンド感度分析 | HonestDiD | — |
| 6 | Oster δ* | robomit::o_test / o_beta | — |
| 7 | 限界効果／傾き | marginaleffects::avg_slopes / plot_slopes | — |
| 7 | 感度分析付き媒介分析 | mediation::mediate + medsens | — |
| 7 | SEM | lavaan::sem | — |
| 8 | 回帰表（任意の形式） | modelsummary | texreg / stargazer |
| 8 | Word表 | flextable / gt::gtsave | officer |
| 8 | LaTeX表のスタイル設定 | kableExtra | — |
| 8 | 係数プロット／イベントスタディ | modelplot / fixest::iplot | ggplot2で手作業 |
| 8 | ビン散布図 | binsreg | — |
| 8 | マルチパネル | cowplot::plot_grid / patchwork | gridExtra |
よくある間違い（および代わりにすべきこと）
| 間違い | 正しい方法 |
|---|---|
大規模なパネルデータにlm(y ~ x + factor(unit) + factor(year))を使う | `feols(y ~ x |
| クラスターデータにデフォルトのiid標準誤差を使う | feols(..., cluster = ~ id)を使います。クラスター数が50未満ならboottestを使います |
| 段階的導入にTWFEを使う | did::att_gt / fixest::sunab / didimputation::did_imputation |
arrange() + group_by()なしでlag(x)を使う | 常にarrange(id, time) %>% group_by(id) %>% mutate(x_l1 = lag(x))を使います |
| 行数を確認せずに結合する | dplyr::*_joinのrelationship引数を使い、その後stopifnot(nrow(df) == n_before)を実行します |
| ロジット係数をそのまま解釈する | AMEにはmarginaleffects::avg_slopes(model)を使います |
| 点推定値だけを報告する | 常にプロットします — modelplot、iplot、plot_slopes |
| 回帰表を手作業で書式設定する | modelsummaryなら1回の呼び出しでLaTeX／Word／HTMLを書き出せます |
| 主要な係数だけを報告する（Table 2なし） | 必ずM1→M6の多列のメイン表を掲載します。これが経済学論文の中心であり、要旨の一文ではありません |
| 図が一つもない係数表 | 経済学の結果には最低限、F1のトレンド + F2のイベントスタディ + F3の係数プロット + F4の感度分析が必要です — Default Output Specを参照してください |
プロットを.pngだけで保存する | LaTeX投稿用に.pdfも保存します |
| スクリプト内でデータセットのパスをハードコードする | here::here()とrenv::init()を使います |
| テストを毎回手動で実行する | targets::tar_make()またはQuartoに組み込みます |
一般的なプロジェクトの骨組み
project/
├── R/
│ ├── 01_clean.R # produces data/analysis.rds
│ ├── 02_transform.R
│ ├── 03_describe.R
│ ├── 04_diagnose.R
│ ├── 05_model.R # saves models to estimates/
│ ├── 06_robust.R
│ ├── 07_further.R
│ └── 08_tables_figures.R
├── data/
│ ├── raw/
│ └── analysis.rds
├── tables/
├── figures/
├── estimates/ # saved fixest objects via saveRDS
├── logs/
├── renv.lock # package versions locked
├── _targets.R # or main.qmd / main.R
└── README.md
最上位に _targets.R（targets パッケージを使用）または main.qmd（Quarto）を置くことで、パイプライン全体を再現可能にできます。
# main.R — minimal driver
source("R/01_clean.R")
source("R/02_transform.R")
source("R/03_describe.R")
source("R/04_diagnose.R")
source("R/05_model.R")
source("R/06_robust.R")
source("R/07_further.R")
source("R/08_tables_figures.R")
Quarto での執筆（説明文・コード・表／図を組み合わせ、PDF／HTML／Word にレンダリング）については、references/08-tables-plots.md §12 を参照してください。
回帰表（modelsummary / etable）クックブック（一ページのレシピ一覧）
modelsummary(...) と fixest::etable(...) は、複数回帰表の基盤となる 2 つの基本要素です。上記の 8 つのパターンは、次のように対応します。
| パターン | 列間で異なるもの | 手順 |
|---|---|---|
| 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 報告の 3 点セット | 第 1 段階／縮約形／2SLS | 5.E — 表 2-quinto |
| F. 因果分析オーケストレーター | 1 列、完全な診断（att_gt / synthdid / causal_forest） | 5.F |
| G. サブグループ表 | サブサンプル（全体／女性／男性／Q1…Q4） | 7 — 表 3 |
| H. 頑健性の総括 | すべての頑健性チェックを積み上げ | 6.j — 表 A1 |
AER の標準スタイルに合わせた modelsummary のデフォルト設定：
modelsummary(
list("(1)" = m1, ..., "(N)" = mN),
output = "tables/tableN.tex", # or .docx / .html
stars = c("*" = 0.1, "**" = 0.05, "***" = 0.01), # AER stars
gof_omit = "BIC|AIC|F|Log|Adj",
coef_map = c("training" = "Training"), # pretty names
notes = c("Cluster-robust SE in parentheses.",
"* p<0.10, ** p<0.05, *** p<0.01.")
)
# For multi-panel paper bundles, use gt::gt_group(modelsummary(...), modelsummary(...))
# or render via Quarto for a single .pdf / .docx / .html target.
図の作成（R による AER 標準図 12 種）
| # | 図 | R コマンド | セクション |
|---|---|---|---|
| 1a | 生のトレンド（DID 図 1） | df %>% group_by(year, treat) %>% summarise(mean(y)) %>% ggplot() | §1 |
| 1b | 処置導入のヒートマップ | panelView::panelview(...) · ggplot + geom_tile | §1 |
| 2a | イベントスタディ係数 | `fixest::iplot(feols(y ~ sunab(G, t) | i + t))` |
| 2a' | Bacon の重み | bacondecomp::bacon + ggplot | §3 |
| 2a'' | CS-DID の動学効果 | did::ggdid(aggte(cs, type="dynamic")) | §3 |
| 2b | 第 1 段階の散布図 | binsreg::binsreg(y=D, x=Z, w=X) | §3（手順 3.5.2） |
| 2c | RD の標準プロット | rdrobust::rdplot(y, x, c=0) | §3（手順 3.5.3） |
| 2c' | McCrary 密度 | rddensity::rdplotdensity(rdd, X) | §3 |
| 2d | マッチングのラブプロット | cobalt::love.plot(MatchIt::matchit(...)) | §3（手順 3.5.4） |
| 2e | SCM の軌跡 | tidysynth::plot_trends · synthdid::plot · Synth::path.plot | §3（手順 3.5.5） |
| 3 | 主な仕様の係数プロット | modelsummary::modelplot(list(m1,...,m6), coefs="training") | §4 |
| 4a | 用量反応 | marginaleffects::plot_predictions(model, condition="dose") | §5 |
| 4b | CATE の分布 | grf::causal_forest(...) + ggplot::geom_histogram(predict(cf)$predictions) | §5 |
| 5 | 仕様曲線 | specr::plot(specr(...))（6.k を参照） | §7 |
| 6 | 感度分析ダッシュボード | HonestDiD::createSensitivityPlot · EValue::evalue | §7（手順 6.l） |
| 7 | 最終的な主図 | 推定量に応じたもの（rdplot、iplot、Synth::path.plot） | §8 |
すべての図は
ggsave()で.png ≥ 300 dpi（スライド／ウェブ用）の両方にエクスポートします。一貫したスタイルにするため、master.Rの先頭でtheme_set(theme_classic(base_size = 11))を一度設定してください。
手法カタログ
古典的 OLS / パネル
library(fixest); library(plm); library(sandwich); library(lmtest)
feols(y ~ X, data = df, cluster = ~ i) # OLS (modern primary)
feols(y ~ X | fe1, data = df, cluster = ~ i) # OLS + 1 FE
feols(y ~ X | fe1 + fe2, data = df, cluster = ~ i) # HD FE workhorse
feols(y ~ X | fe1 + fe2, data = df, cluster = ~ fe1 + fe2) # 2-way cluster
fepois(count ~ X | fe1 + fe2, data = df, cluster = ~ i) # Poisson + FE
feglm (y ~ X | fe1, data = df, family = binomial(link = "logit"),
cluster = ~ i) # Logit + FE
plm (y ~ X, data = df, model = "within", index = c("i","t")) # panel FE
plm (y ~ X, data = df, model = "random", index = c("i","t")) # RE (Hausman: phtest)
差分の差分法
library(fixest); library(did); library(didimputation); library(synthdid); library(bacondecomp); library(HonestDiD); library(DIDmultiplegtDYN)
feols(y ~ i(treated, post, ref = 0) | i + t, df, cluster = ~ i) # 2×2
feols(y ~ sunab(first_treat, year) | i + year, df, cluster = ~ i) # SA event study
att_gt(yname="y", tname="t", idname="i", gname="G", data=df,
control_group="nevertreated", est_method="dr", clustervars="i") # CS-DID
did_imputation(data=df, yname="y", gname="G", tname="t", idname="i",
horizon=0:5, pretrends=-5:-1, cluster_var="i") # BJS imputation
DIDmultiplegtDYN(df, "y", "i", "t", "training", effects=5, placebo=3) # de Chaisemartin
synthdid_estimate(panel.matrices(df,"i","t","y","training"), ...) # synthetic DID
bacon(y ~ training, data=df, id_var="i", time_var="t") # TWFE diagnostic
HonestDiD::createSensitivityResults(...) # PT sensitivity
操作変数法 / 2SLS
library(fixest); library(AER); library(ivreg)
feols(y ~ X | D ~ Z, df, cluster = ~ firm_id) # workhorse w/ HD FE
fitstat(iv, ~ ivf + ivwald + sargan + cd) # CD/KP/Sargan/F
AER::ivreg(y ~ D + X | Z + X, data = df) # classic API
summary(iv, vcov. = sandwich, diagnostics = TRUE) # with diagnostics
回帰不連続デザイン
library(rdrobust); library(rddensity); library(rdmulti)
rdrobust(y, x, c = 0, kernel = "triangular", bwselect = "mserd") # Sharp RD
rdrobust(y, x, c = 0, fuzzy = D) # Fuzzy RD
rddensity(X = x, c = 0) # McCrary density
rdplot(y, x, c = 0)
rdmc(y, x, cutoffs = c(0, 5, 10)) # multi-cutoff
マッチング / 再重み付け
library(MatchIt); library(WeightIt); library(cobalt)
matchit (D ~ X1 + X2, data = df, method = "nearest", ratio = 1) # PSM
matchit (D ~ X1 + X2, data = df, method = "cem") # Coarsened EM
weightit(D ~ X1 + X2, data = df, method = "ebal") # entropy balancing
weightit(D ~ X1 + X2, data = df, method = "ps", estimand = "ATE") # IPW
love.plot(matchit_obj, threshold = 0.10) # SMD diagnostic
合成コントロール法
library(Synth); library(gsynth); library(tidysynth); library(synthdid)
Synth::synth(...) # ADH SCM
gsynth(y ~ training, data = df, index = c("i","t"), force = "two-way") # generalized SC
synthdid_estimate(panel.matrices(...)) # synthetic DID
tidysynth::synthetic_control(df, ...) %>% generate_predictor(...) %>%
generate_weights() %>% generate_control()
機械学習による因果推論（モード B — §B を参照）
library(grf); library(DoubleML); library(mlr3); library(causalDML)
causal_forest(X, Y, W, num.trees = 4000, honesty = TRUE) # GRF causal forest
DoubleML::DoubleMLPLR$new(data, ml_l = lrn("regr.ranger"),
ml_m = lrn("regr.ranger")) # DML PLR
DoubleML::DoubleMLIRM$new(data, ...) # DML interactive
predict(cf)$predictions # CATE per row
average_treatment_effect(cf, target.sample = "treated")
test_calibration(cf); variable_importance(cf)
policytree::policy_tree(X, gamma, depth = 3) # policy tree
頑健性、感度分析、推論
library(fwildclusterboot); library(ri2); library(multcomp); library(robomit); library(EValue)
boottest(model, param = "training", clustid = "state", B = 9999) # wild cluster bootstrap
ri2::conduct_ri(...) # randomization inference
robomit::o_test(...) # Oster δ
EValue::evalue(RR(1.45), lo = 1.10, hi = 1.91) # E-value
fwildclusterboot::boottest(..., type = "rademacher") # alt bootstrap dist
生存時間分析 / 疫学（モード A — §A を参照）
library(survival); library(survminer); library(survRM2); library(ipw); library(tmle); library(zelig)
survfit(Surv(time, event) ~ A, data = df) # KM
coxph (Surv(time, event) ~ A + X, data = df) # Cox
survreg(Surv(time, event) ~ A + X, data = df, dist = "weibull") # AFT
rmst2 (time, status, arm, tau = 1825) # RMST contrast
ipw::ipwpoint(...) # IPTW
tmle (Y, A, W = X, ...) # TMLE
gfoRmula::gformula_survival(...) # parametric g-formula
TwoSampleMR::mr(...) # Mendelian randomization
他のスキルへの引き継ぎが必要な場合
- エージェントネイティブな単一インポートの Python ワークフロー（
import statspai as sp）→00-StatsPAI_skill。 - Python の従来型スタックを明示的に使う場合 →
00.1-Full-empirical-analysis-skill。 - Stata の
.doパイプライン →00.2-Full-empirical-analysis-skill_Stata。 - 複数言語対応の Mixtape テンプレート（Python/R/Stata を並べて表示）→
10-Jill0099-causal-inference-mixtape。 - ベイズ統計の R ワークフロー（
brms/rstan/cmdstanr）→23-Learning-Bayesian-Statistics-baygent-skills。 - 分析後の論文執筆 → このリポジトリの執筆スキル。
このスキルはステップ 8 で終了します（.tex / .docx の表と .pdf の図まで）。論文の執筆は対象外です。
ステップ 8.5 — 原稿への引き継ぎ（終了時の取り決め）
このスキルの範囲は、引き続きステップ 8 までです。ただし、「終了」は「行き止まり」を意味するべきではありません。 執筆段階では、ばらばらのファイルが入ったフォルダーではなく、機械可読な取り決めが必要です。 終了する前に、表と図の隣に 2 つの成果物を出力してください。
1. exhibits_index.md — 番号付きの各図表につき 1 行：
| 図表 | ファイル | 裏付ける主張 | 作成元 |
|---|---|---|---|
| 表 2 | tables/table2_main.{tex,docx,xlsx} | 政策が対数アウトプットに及ぼす主効果 | analysis.R §baseline |
| 図 2 | figures/fig2_event_study.{png,pdf} | 差のある事前トレンドがないこと | analysis.R §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により、本文、上記のすべての図表、図、参考文献一覧を 1 つのファイルにまとめた09_submission/main.docxを、まさにこの 2 つの成果物から組み立てます。ステージ 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 対応