--- name: Full Empirical Analysis R slug: full-empirical-analysis-r category: Data description: Full Empirical Analysis R runs a complete 8-step empirical research pipeline in R using tidyverse and econometrics packages, generating publication-ready regression tables, descriptive statistics, robustness checks, and figures for applied economics, epidemiology, and ML causal inference. github: "https://github.com/brycewang-stanford/Auto-Empirical-Research-Skills/tree/main/skills/00.3-Full-empirical-analysis-skill_R" language: Stata stars: 3419 forks: 445 install: "npx degit https://github.com/brycewang-stanford/Auto-Empirical-Research-Skills/tree/main/skills/00.3-Full-empirical-analysis-skill_R ~/.claude/skills/00.3-Full-empirical-analysis-skill_R" installs_to: ~/.claude/skills/00.3-Full-empirical-analysis-skill_R source_path: skills/00.3-Full-empirical-analysis-skill_R/SKILL.md collection_size: 25 category_size: 668 collection_url: "https://dirskills.com/collections/brycewang-stanford/Auto-Empirical-Research-Skills" added: 2026-08-16T07:02:46.950Z last_synced: 2026-08-16T07:02:46.950Z canonical_url: "https://dirskills.com/skills/full-empirical-analysis-r" --- # Full Empirical Analysis R Full Empirical Analysis R runs a complete 8-step empirical research pipeline in R using tidyverse and econometrics packages, generating publication-ready regression tables, descriptive statistics, robustness checks, and figures for applied economics, epidemiology, and ML causal inference. **Install:** ```bash npx degit https://github.com/brycewang-stanford/Auto-Empirical-Research-Skills/tree/main/skills/00.3-Full-empirical-analysis-skill_R ~/.claude/skills/00.3-Full-empirical-analysis-skill_R ``` ## README # 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/`](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: ```text best_skill: train_signal: selection_split: heldout_gate: accepted_patterns: rejected_patterns: patch_scope: reject_if: ``` 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 ```r # 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-.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`](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