Full Empirical Analysis — Classical R Workflow

This skill is the canonical 8-step pipeline an applied economist runs on every empirical paper, written in the modern tidyverse + econometrics R ecosystem — dplyr/tidyr/haven for data, fixest as the panel/IV/DID workhorse, did/bacondecomp/HonestDiD for modern DID, rdrobust/rddensity for RD, Synth/gsynth/synthdid for synthetic control, MatchIt/WeightIt/cobalt/ebal for matching, grf/DoubleML for ML causal, mediation for causal mediation, marginaleffects for post-estimation, modelsummary/kableExtra/gt for publication tables, ggplot2/iplot/binsreg for figures.

Companion skills: this is the R sibling of 00-StatsPAI_skill (Python DSL), 00.1-Full-empirical-analysis-skill (explicit Python), and 00.2-Full-empirical-analysis-skill_Stata (Stata .do). All four implement the same 8 steps, in their respective ecosystems.

Philosophy

1. Tidyverse + fixest, the modern R idioms. feols(... | unit + year, cluster = ~unit), not Frankenstein-y lm(y ~ x + factor(unit) + factor(year)).
2. Reproducible scripts / Quarto. Every example below is paste-runnable. renv for package locking; Quarto (.qmd) for combined narrative + code + tables/figures.
3. 8 steps, first-class. R users historically over-invest in Step 5; this skill treats Steps 1–4 and 6–8 as core.
4. Rich outputs. Every step yields at least one table or figure — tex/docx/png/pdf.
5. Progressive disclosure. SKILL.md gives the canonical call per step; references/ holds variant-specific depth.

SkillOpt-style execution gate

Use this long playbook as a seed skill, not as a script to exhaustively apply. SkillOpt discipline: treat each local R/Quarto change as a candidate patch that must beat a selection check and survive a held-out check before it becomes reusable boilerplate. Before writing or revising an R script/Quarto workflow, compress the user's request into a task-local best_skill card:

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>

1. Route card: record the mode (econ, epi, or ml-causal), estimand, identification design, focal outcome/treatment, R package family, and required artifacts.
2. Bounded edit: change one decision at a time (sample rule, estimator, clustering, export format, or robustness check). Prefer the smallest patch that can pass validation.
3. Selection split discipline: treat the user's immediate failure or requested artifact as the selection split. Reserve at least one alternate outcome, sample window, estimator family, or export target as the held-out gate.
4. Held-out gate: define checks before running code: row counts, distinct() key uniqueness, treatment support, missingness thresholds, expected table/figure files, and one non-focal robustness/specification that the change must not break.
5. Reject buffer: if a candidate spec fails the gate, log the failure, R/Quarto diff, and gate output in analysis_log.md; revert to the last passing spec and do not retry the same unchecked pattern.
6. Slow/meta update: at the end of the task, write down accepted_patterns and rejected_patterns from the trajectory. Do not widen the canonical project template from a single passing run.
7. Promote only after validation: only turn a one-off fix into reusable project boilerplate after it passes the current data and at least one alternate outcome/sample/specification.

---

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

The default playbook above is AER-style applied econometrics — the AEA convention: written-out estimating equation, identifying assumption, design horse-race, full robustness gauntlet. The skill also ships two parallel sub-pipelines for the other two big causal-inference traditions, each reusing the same Steps 1–4 (cleaning / construction / Table 1 / diagnostics) and Step 8 (tables/figures) — only Step 5 (estimator) and Step 6/7 swap packages:

Mode	Reader convention	Step-5 estimator stack	Reporting stack	Jump to
Default — Applied Econ (AER / QJE / AEJ)	"Show the equation + identifying assumption + design horse-race; controls visible; clustered SE"	DID / IV / RD / SCM / matching / fixest::feols HDFE	AER house-style multi-column modelsummary + kableExtra / gt / flextable + 8-section paper layout	Steps 1 → 8 (entire playbook below)
Mode A — Epidemiology / Public Health	"STROBE / TRIPOD-AI; target trial protocol; doubly-robust estimand; absolute & relative risk; KM survival"	Target-trial emulation · IPTW (WeightIt / PSweight) · g-formula (gfoRmula) · TMLE (tmle / ltmle) · Mendelian randomization (MendelianRandomization / TwoSampleMR / MRPRESSO) · KM / Cox / AFT (survival / survminer / flexsurv)	Same modelsummary + risk-difference / hazard-ratio / E-value rows	§A. Epidemiology pipeline
Mode B — ML Causal Inference	"DML / meta-learners / causal forest / DR-learner; CATE distribution; policy value"	DML (DoubleML) · S/T/X/R/DR-Learner (causalweight / grf) · GRF causal forest (grf::causal_forest) · BART/BCF (bartCause / bcf) · matrix completion (MCPanel)	modelsummary ML horse-race + grf CATE plot + policy-value table + conformalInference PI	§B. ML causal pipeline

How to invoke a non-default mode (Claude / agent picks this up from the user's wording):

User says...	Mode the skill switches to
"Run a DID / IV / RD / event study", "AER table", "applied micro"	Default (AER econ) — Steps 1 → 8
"Target trial emulation", "g-formula", "IPTW", "TMLE", "Mendelian randomization", "STROBE / TRIPOD", "公共健康 / 流行病学", "epi pipeline", "RWE study", "cohort study", "case-control"	Mode A (Epi) — §A
"DML", "double machine learning", "causal forest", "meta-learner", "CATE", "BCF", "policytree", "policy learning", "conformal causal", "fairness audit", "ML causal", "uplift modeling", "因果机器学习"	Mode B (ML causal) — §B
"Mix" (e.g. "estimate DID + then ML CATE on the heterogeneity")	Default + Mode B in sequence — every estimator yields a coefficient + SE pair, drop them all into one modelsummary(...) for the horse-race column

The three modes share the same Step 1–4 cleaning / Table 1 / diagnostics scaffolding, the same Step 8 export stack, and the same DAG-first identification logic — switching modes only changes which Step-5 estimator family you reach for, not the surrounding paper structure. If you only want descriptive stats / Table 1 / a balance check, the AER gtsummary::tbl_summary / modelsummary::datasummary_balance calls in Step 3 work identically across all three modes.

---

Default Output Spec — Economics Empirical Paper

This skill defaults to the applied-economics paper convention. Unless the user explicitly asks for a single point estimate, every run produces the full publication-ready output set below. Treat it as the contract of Step 8 — mandatory, not opt-in.

Required tables (always produced)

#	Table	R source	Saves to
T1	Summary statistics & balance (treated vs control, with SMD / p-values)	gtsummary::tbl_summary + add_p + add_difference (Step 3)	tables/table1_balance.xlsx + .docx + .tex
T2 ★	Main results — multi-column regression M1→M6 (progressive controls + FE)	fixest::feols × 6 specs → modelsummary (Step 5–6)	tables/table2_main.xlsx + .docx + .tex
T3	Mechanism / outcome ladder — same treatment, 3+ outcomes side-by-side	loop feols over y ∈ {Y1, Y2, Y3, Y_main} → modelsummary (Step 7)	tables/table3_mechanism.xlsx + .docx + .tex
T4	Heterogeneity — subgroup × main coef (gender, age, region, …)	subgroup feols × linearHypothesis → modelsummary (Step 7)	tables/table4_heterogeneity.xlsx + .docx + .tex
T5	Robustness battery — alt SE / cluster / sample / placebo, in one table	feols × variants → modelsummary (Step 6)	tables/table5_robustness.xlsx + .docx + .tex

★ Table 2 is the centerpiece of every economics paper. It is the multi-column regression table that walks the reader from raw correlation (M1) to the fully-specified design (M6: 2-way FE + interacted FE + cluster-robust SE). Do not collapse it into a single column. Do not report only the headline coefficient. The progression is the credibility argument: if M1→M6 is monotone and stable, the design is plausibly identifying; if it collapses on adding FE, that is the result.

Canonical 6 columns, in order:

1. M1 raw bivariate (feols(y ~ treat, data))
2. M2 + demographics (+ age + edu)
3. M3 + sector controls (+ tenure / firm_size)
4. M4 + unit FE (| worker_id)
5. M5 + 2-way FE (| worker_id + year)
6. M6 + interacted FE (| worker_id + year + industry^year) with cluster = ~ worker_id

Required figures (always produced)

#	Figure	R source	Saves to
F1	Trend / motivation — treated vs control over time, with policy line	dplyr group means → ggplot + geom_line (Step 3)	figures/fig1_trend.png (300 dpi, 必须导出 PNG) + .pdf
F2	Event-study coefficients with 95% CI, base period at –1	fixest::sunab() / did::ggdid / iplot (Step 5)	figures/fig2_event_study.png (300 dpi, 必须导出 PNG) + .pdf
F3	Coefficient plot across specs M1→M6	modelsummary::modelplot() (Step 8)	figures/fig3_coefplot.png (300 dpi, 必须导出 PNG) + .pdf
F4	Robustness / sensitivity — bacondecomp::bacon plot, HonestDiD::createSensitivityPlot, or spec curve	scenario-specific (Step 6)	figures/fig4_sensitivity.png (300 dpi, 必须导出 PNG) + .pdf

Output file layout (default)

project/
├── tables/    table1_balance.xlsx/.docx/.tex  table2_main.xlsx/.docx/.tex
│              table3_mechanism.xlsx/.docx/.tex table4_heterogeneity.xlsx/.docx/.tex
│              table5_robustness.xlsx/.docx/.tex
└── figures/   fig1_trend.png(300dpi)+.pdf      fig2_event_study.png(300dpi)+.pdf
               fig3_coefplot.png(300dpi)+.pdf   fig4_sensitivity.png(300dpi)+.pdf

关键输出规则（必须遵守）：

• 图片格式：所有图片必须同时导出 PNG 格式（≥300 dpi） 和 PDF 格式（用于 LaTeX 排版）
• 表格格式：所有回归表格必须同时导出 Excel（.xlsx）、Word（.docx） 和 LaTeX（.tex） 三种格式
• PNG 用于幻灯片、Markdown 文档、邮件等场景；PDF 用于学术论文排版

When to deviate

• Single quick estimate — produce only the relevant cell, but warn that the standard deliverable is the full set above and offer to run it.
• Design does not support a figure (cross-section → no event study) — skip with a printed message() explaining why; do not silently drop.
• N=1 treated unit (Synth / synthdid) — replace F1/F2 with the SCM trajectory + placebo distribution; T1–T5 still apply.

---

Required packages

# 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.

---

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

┌──────────────────────────────────────────────────────────────────────┐
│ 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   │
└──────────────────────────────────────────────────────────────────────┘

The 8 steps mirror the canonical sections of an applied AER / QJE / AEJ paper. Each step is one paper section and emits a paper-ready artifact on disk:

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

Below is the canonical call at each step. All examples share one running narrative — labor-econ panel where training (treatment) affects log_wage (outcome), with covariates age, edu, tenure, panel keys worker_id/firm_id/year. Variable names and parameter values are illustrative.

When a step has many variants (5 staggered-DID estimators; 4 hetero tests), SKILL.md shows the one you reach for first; deeper variants live in references/NN-<topic>.md.

---

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

A modern AER paper has 5–7 figures and 3–5 main tables + an appendix robustness table. Every step below leaves at least one numbered artifact on disk. Default file names assume parallel .tex / .docx / .xlsx exports (the agent should produce all three so co-authors can edit in Word, the build system can use LaTeX, and editors can edit raw numbers in Excel). 所有图片必须同时保存 PNG（≥300 dpi）和 PDF 两种格式。

§	Artifact	R primitive	Filenames
§1	Figure 1: raw trends / treatment rollout	df %>% group_by(year, treat) %>% summarise(mean(y)) %>% ggplot()	figures/fig1_trend.png(300dpi)+.pdf
§1	Table 1: summary stats (full / treated / control + Δ + SMD)	gtsummary::tbl_summary · modelsummary::datasummary_balance	tables/table1_balance.xlsx/.docx/.tex
§3	Figure 2: identification graphic (event-study / first-stage / McCrary / RD scatter / SCM trajectory)	fixest::iplot(es) · binsreg · rdrobust::rdplot · rddensity · Synth::path.plot	figures/fig2_event_study.png(300dpi)+.pdf
§4	Table 2: main results — progressive controls M1→M6	modelsummary(list("(1)"=m1,...,"(6)"=m6)) · fixest::etable	tables/table2_main.xlsx/.docx/.tex
§4	Table 2-bis: design horse-race (OLS / IV / DID / DML)	modelsummary(list("OLS"=ols, "2SLS"=iv, "CS-DID"=cs, "DML"=dml))	tables/table2b_designs.xlsx/.docx/.tex
§4	Figure 3: coefficient plot across specs	modelplot(list(m1,...,m6), coef_map="training")	figures/fig3_coefplot.png(300dpi)+.pdf
§5	Table 3: heterogeneity by subgroup	modelsummary(g_full, g_male, g_fem, g_q1, ..., g_q4)	tables/table3_heterogeneity.xlsx/.docx/.tex
§5	Figure 4: dose-response / CATE	marginaleffects::plot_predictions · grf::plot.causal_forest	figures/fig4_cate.png(300dpi)+.pdf
§6	Table 4: mechanism / outcome ladder	loop feols over outcomes → modelsummary	tables/table4_mechanism.xlsx/.docx/.tex
§7	Table A1: robustness master (one column per check)	modelsummary(list(base, no99, balpan, dropearly, wfe, cl2way, logy, ihsy, psm, ebal))	tables/tableA1_robustness.xlsx/.docx/.tex
§7	Figure 5: spec curve	specr::specr() + plot_specs (or hand-rolled purrr::pmap)	figures/fig5_spec_curve.png(300dpi)+.pdf
§7	Figure 6: sensitivity (HonestDiD / Oster / E-value)	HonestDiD::createSensitivityPlot · robomit::o_test · EValue	figures/fig6_sensitivity.png(300dpi)+.pdf
§8	Replication bundle: all tables in one document	modelsummary(..., output="docx") · gt::gtsave() · Quarto / Rmd	replication/paper_tables.xlsx/.docx/.tex

Every R estimator above (fixest::feols / AER::ivreg / did::att_gt / grf::causal_forest / synthdid_estimate) returns a result object that can be passed straight into modelsummary(...) / modelplot(...) / etable(...). Don't hand-roll LaTeX from kable(), and don't render Word via flextable directly — modelsummary, etable, and gtsummary apply book-tab borders, AER stars, and the right SE label automatically. For deeper export recipes, see references/08-tables-plots.md.

---

Export cookbook — LaTeX / Word / Excel in one block

关键规则（必须遵守）：每个表格必须同时导出三种格式——Excel(.xlsx)、Word(.docx)、LaTeX(.tex)。每个图片必须同时保存PNG(≥300dpi)和PDF两种格式。

R has the best publication-table ecosystem of the three languages. Three tiers, picked by scope:

Tier	Use when	API	Hot args
1. Single multi-column table	Exporting one Table 2 / Table 3 / Table A1 with progressive columns	`modelsummary(list("(1)"=m1,...,"(N)"=mN), output="tables/tab.tex", stars=c(""=.1,""=.05,""=.01), gof_omit="BIC	AIC
2. Multi-panel paper format (Tables 2 + 3 + A1 + A2 in one file)	Producing the paper-tables block — main + heterogeneity + robustness + placebo as a single document	modelsummary chained with gt::gt_group() for one document with section headers, OR Quarto .qmd rendering multiple modelsummary calls between prose	gt_group(modelsummary(...), modelsummary(...)) · quarto render paper.qmd
3. Full session bundle (the Stata collect / Python Stargazer + pylatex equivalent)	Replication appendix that mixes summary stats + balance + multiple regression tables + headings + prose in one file	Quarto is the modern R-native answer. master.qmd interleaves prose + chunks that emit modelsummary / gtsummary / ggplot2 outputs; one quarto render produces .pdf / .docx / .html	YAML front matter sets format: [pdf, docx, html] for triple-target output

Journal styling — pick the right stars and SE label. The AEA convention is c("*"=.1, "**"=.05, "***"=.01) and notes = "Cluster-robust standard errors in parentheses...". Define a wrapper once at the top of 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
    )
  }
}

For the multi-panel .docx / .xlsx and Quarto cookbook (single-file paper-tables bundle), see references/08-tables-plots.md.

---

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

Before touching the data, write down (a) the population, (b) the design, (c) the minimum detectable effect (MDE) under the planned sample size and α=0.05, β=0.20. Persist the result as pap.json so a referee can verify the design was powered before, not after, the data were seen.

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