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
- Tidyverse + fixest, the modern R idioms.
feols(... | unit + year, cluster = ~unit), not Frankenstein-y lm(y ~ x + factor(unit) + factor(year)).
- Reproducible scripts / Quarto. Every example below is paste-runnable.
renv for package locking; Quarto (.qmd) for combined narrative + code + tables/figures.
- 8 steps, first-class. R users historically over-invest in Step 5; this skill treats Steps 1–4 and 6–8 as core.
- Rich outputs. Every step yields at least one table or figure — tex/docx/png/pdf.
- 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>
- Route card: record the mode (
econ, epi, or ml-causal), estimand, identification design, focal outcome/treatment, R package family, and required artifacts.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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:
- M1 raw bivariate (
feols(y ~ treat, data))
- M2 + demographics (
+ age + edu)
- M3 + sector controls (
+ tenure / firm_size)
- M4 + unit FE (
| worker_id)
- M5 + 2-way FE (
| worker_id + year)
- 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
install.packages(c(
"tidyverse", "haven", "readxl", "data.table", "janitor",
"naniar", "VIM", "mice", "validate",
"gtsummary", "tableone", "modelsummary", "kableExtra", "gt",
"stargazer", "texreg", "flextable", "psych", "summarytools",
"lmtest", "sandwich", "car", "tseries", "urca", "plm",
"clubSandwich"
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)+ |
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:
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
gof_omit
coef_map coef_map
notes paste
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)
library(WebPower)
library(jsonlite)
pwr.t.test(d = 0.20, power = 0.80, sig.level = 0.05,
type = "two.sample", alternative = "two.sided")
pwr.t.test(n = 2000, power = 0.80, sig.level = 0.05,
type = "two.sample")$d
WebPower::wp.crt2arm(f = 0.20, J = NULL n icc power
alpha alternative
pap
population
treatment
outcome
estimand
design
alpha
power_target
mde_d
n_planned
frozen_at
git_sha
write_jsonpap pretty auto_unbox
For richer DAG-aware power analysis (write down the DAG, declare estimands, simulate the design), use DeclareDesign — it is the R-native equivalent of EGAP's pre-analysis flow.
Commit artifacts/pap.json in the repo before Step 1. AEA RCT Registry / OSF preregistration tools accept it as the analysis-plan exhibit.
Step 0 — Sample-construction log & 5-check data contract
An AER §1 Data section has three jobs: (a) describe sources, (b) document every sample restriction (the "footnote 4" sample log), (c) lock the panel structure.
0.1 Sample-construction log (footnote 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_rowstep n nrowdf2
df3 df2 filterindustry
sample_log sample_log add_rowstep n nrowdf3
df df3
printsample_log
write_jsonsample_log pretty
Paste the printed tibble verbatim as footnote 4 of the paper.
0.2 Five-check data contract (go / no-go gate)
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),
dtypes = sapply(df[keys], function(x) class(x)[1]),
n_missing = sapplydfkeys x x
n_dupes_on_keys id time
duplicateddf id time
panel_balanced
cohort_sizes
id time
bal df count.dataid
contractpanel_balanced baln baln
contractn_dropped_by_balance baln baln
df
contractcohort_sizes df distinct.dataid .keep_all
countfirst_treat deframe
contracty_range dfy na.rm
contracttreatment_share meandftreatment na.rm
miss_y dfy
contractmcar_hint
miss_y miss_y
cov covariates
dfcov
p t.testdfcovmiss_y dfcovmiss_yp.value
p
contractmcar_hint sprintf
cov p
contract
contract data_contractdf y treatment
id time
covariates
stopifnotcontractn_dupes_on_keys
stopifnotcontractn_missing
write_jsoncontract
pretty auto_unbox
If any stopifnot fires, stop and fix it in dplyr first. R estimators silently drop NA rows downstream — this contract is the cheapest insurance against "why did N drop from 12,000 to 9,800 between Table 1 and Table 2?" referee questions.
Step 1 — Data import & cleaning
Deeper patterns: references/01-data-cleaning.md — every format (haven/readxl/data.table::fread/arrow::read_parquet/DBI), janitor::clean_names, naniar missingness viz, MCAR/MAR/MNAR triage with mice, validation with validate/assertr, panel structure checks.
library(tidyverse)
library(haven)
library(janitor)
library(naniar)
library(skimr)
df <- read_dta("raw/panel.dta") %>%
clean_names()
skim(df)
naniar::miss_var_summary(df)
naniar::vis_miss(df)
df <- df %>%
mutate(
year = as.integer(year),
wage = as.numeric(wage),
gender = as.factor(gender),
date as.Datedate
key_vars
df df
drop_naall_ofkey_vars
cat nrowdf
df df
mutate
tenure_missing tenure
tenure if_elsetenure mediantenure na.rm tenure
union fct_explicit_naas.factorunion na_level
df df
mutatewage_z scalewage
outlier_z4 wage_z
cat dfoutlier_z4 na.rm
stopifnotnrowdf distinctworker_id year nrowdf
firm_chars read_dta
n_before nrowdf
df df
left_joinfirm_chars by relationship
stopifnotnrowdf n_before
df countyear
df countworker_id summary
Key principle: dplyr + explicit stopifnot() assertions. No silent row drops downstream.
Step 2 — Variable construction & transformation
Deeper patterns: references/02-data-transformation.md — log/IHS/Box–Cox via MASS::boxcox, group winsorization with dplyr, scale() and bestNormalize, factor handling, lag/lead with dplyr::lag, panel timing.
library(DescTools)
df <- df %>%
mutate(
log_wage = log(pmax(wage, 1)),
ihs_assets = asinh(assets),
wage_w1 = DescTools::Winsorize(wage, probs = c(0.01, 0.99), na.rm = TRUE),
age_std = as.numeric(scale(age)),
age_sq = age^2,
trt_x_edu = training * edu
) %>%
group_byindustry year
mutatewage_w1_iy DescToolsWinsorizewage probs
na.rm
ungroup
arrangeworker_id year
group_byworker_id
mutate
log_wage_l1 laglog_wage
log_wage_f1 leadlog_wage
d_log_wage log_wage laglog_wage
wage_mean_i meanlog_wage na.rm
log_wage_dm log_wage wage_mean_i
ungroup
group_byworker_id
mutatefirst_treat ifelsetraining
yeartraining
ungroup
mutaterel_time year first_treat
never_treated first_treat
cpi read_csv
df df
left_joincpi by
mutatecpi_base cpiyear
wage_real wage cpi_base cpi
log_wage_real pmaxwage_real
Step 2.5 — Empirical strategy (write the equation + identifying assumption)
This is the heart of an AER paper. Before any code, write down the equation explicitly and state the identifying assumption. Vague identification language is the single most common reason a referee rejects an applied paper. Persist the strategy as strategy.md so it is a dated, version-controlled artifact — not a post-hoc rationalization written after seeing the coefficient.
Equation × identifying assumption × R estimator (decision table)
| Design | Estimating equation | Identifying assumption | R estimator |
|---|
| 2×2 DID | Y_it = α_i + λ_t + β·D_it + X'γ + ε_it | parallel trends conditional on X | `feols(y ~ i(treated, post, ref=0) |
| Event-study (CS / SA) | Y_it = α_i + λ_t + Σ_{e≠-1} β_e · 1{t-G_i = e} + ε_it | no anticipation + group-time PT | `feols(y ~ sunab(G, t) |
| 2SLS | Y_i = α + β·D_i + X'γ + ε_i; D_i = π·Z_i + X'δ + u_i | exclusion + relevance + monotonicity | `feols(y ~ X |
| Sharp RD | Y_i = α + β·1{X_i ≥ c} + f(X_i) + ε_i (local poly) | continuity of E[Y(0)|X] at c, no manipulation | rdrobust::rdrobust(y, x, c=0) (+ rddensity) |
| SCM | Ŷ_1t(0) = Σ_j ŵ_j Y_jt, τ_t = Y_1t − Ŷ_1t(0) for t≥T_0 | pre-period fit + interpolation validity | Synth::synth · gsynth::gsynth · synthdid::synthdid_estimate · tidysynth |
| Selection-on-observables (matching/IPW/DML) | Y_i = m(X_i) + β·D_i + ε_i (Robinson partialling-out) | unconfoundedness + overlap | MatchIt::matchit + lm · WeightIt · DoubleML::DoubleMLPLR · grf::causal_forest |
Design picker (when the user is unsure)
┌─ 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)
Pre-registration strategy.md template
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")
Commit artifacts/strategy.md in the repo before running Step 5 / Step 6. The git log of this file is the analysis plan.
Step 3 — Descriptive statistics & Table 1
Deeper patterns: references/03-descriptive-stats.md — gtsummary::tbl_summary (the modern Table 1 standard), modelsummary::datasummary_balance with SMDs, tableone::CreateTableOne, correlation matrices with significance via corrplot / psych::corr.test, distribution plots via ggplot2.
library(gtsummary)
library(modelsummary)
df %>%
select(log_wage, age, edu, tenure, training) %>%
datasummary_skim()
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()
kableExtrasave_kable
df
selectlog_wage age edu tenure female training
tbl_summaryby training
add_p
add_difference
add_n
modify_headerlabel
bold_labels
as_gt
gtgtsave
datasummary_balance training
data df selecttraining age edu tenure female
output
librarycorrplot; librarypsych
corr_obj corr.testdf selectlog_wage age edu tenure training
method
corrplotcorr_objr method type
p.mat corr_objp sig.level insig
addCoef.col number.cex
tl.col tl.srt
col colorRampPalette
libraryggplot2
p1 ggplotdf aeslog_wage fill factortraining
geom_densityalpha
scale_fill_manualvalues
labels name
labsx y
title
theme_classic
p2 ggplotdf aessample log_wage
stat_qq stat_qq_line
labstitle theme_classic
cowplotplot_gridp1 p2 labels
ggsave plot . width height
df
group_byyear training
summarisemean_log_wage meanlog_wage na.rm .groups
ggplotaesyear mean_log_wage color factortraining
geom_linelinewidth geom_pointsize
geom_vlinexintercept policy_year linetype
scale_color_manualvalues
labels name
labsx y theme_classic
ggsave width height
Step 3.5 — Identification graphics (Section "Identification, graphical evidence")
AER convention: the identification figure precedes the regression table. The reader should see graphical evidence that PT holds / first stage is strong / RD jumps cleanly before you ask them to trust your point estimate.
3.5.1 Event-study figure + numerical pre-trends test (DID identification)
Pre-period coefficients ≈ 0 (with the −1 reference period normalized to zero) is the visual evidence for parallel trends. Pair the figure with a numerical pre-trends test so reviewers don't have to eyeball it.
library(fixest); library(ggplot2)
es <- feols(log_wage ~ sunab(first_treat, year) | worker_id + year,
data = df, cluster = ~ worker_id)
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)
pre_idx grep coefesgrepl coefes
W waldes coefespre_idx
catsprintf Wstat Wp
librarybacondecomp
bd baconlog_wage training data df
id_var time_var
ggplotbd aesweight estimate color type shape type
geom_pointsize
labstitle
x y
ggsave width height
librarydid
cs att_gtyname tname idname
gname data df
control_group est_method
clustervars
ggdidaggtecs type
labstitle
ggsave width height
3.5.2 First-stage F-statistic + scatter (IV identification)
Rule of thumb: first-stage F ≥ 10 for OLS-style inference; F ≥ 23 for AR-equivalent inference (Stock–Yogo / Lee 2022). fixest::feols reports F automatically; AER::ivreg requires 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)
library(binsreg)
binsreg(y = df$training, x = df$Z1, w = df[, c("age","edu")],
nbins = 20, polyreg = 2, ci = c(
ggsave width height
3.5.3 RD: McCrary density + canonical RD plot
The signature RD figure is rdplot (CCT-style binned scatter with local-polynomial fit on each side), paired with the McCrary manipulation test.
library(rdrobust); library(rddensity)
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)
rdd <- rddensity(X = df$running_var, c = 0)
print(summary(rdd))
rdplotdensity(rdd, X = df$running_var
title
ggsave width height
3.5.4 Matching: love plot (standardized differences pre vs post)
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: synthetic-control trajectory + gap plot
For synthetic-control designs the canonical Figure 2 is the treated-vs-synthetic time series with treatment time annotated.
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_differencessc; ggsave width height
librarysynthdid
sdid_setup panel.matricesdf unit time
outcome treatment
sdid_fit synthdid_estimatesdid_setupY sdid_setupN0 sdid_setupT0
plotsdid_fit control.name
ggsave width height
Identification-specific checks (PT for DID, weak-IV F, density for RD, common support for matching) are also auto-run inside the Step-5 estimators — don't duplicate the numerics here, but DO produce the figures: a referee scans the figures first.
Step 4 — Diagnostic statistical tests
Deeper patterns: references/04-statistical-tests.md — every classical test. lmtest/sandwich/car/tseries/urca/plm.
library(lmtest)
library(sandwich)
library(car)
library(tseries)
library(urca)
ols <- lm(log_wage ~ training + age + edu + tenure, data = df)
shapiro.test(sample(residuals(ols), min(5000, length(residuals(ols)))))
tseries::jarque.bera.test(residuals(ols))
bptest(ols)
bptest(ols, ~ I(fitted(ols)^ . data df
dwtestols
bgtestols order
Box.testresidualsols lag type
libraryplm
pdata pdata.framedf index
plm_fe plmlog_wage training age edu data pdata model
pbgtestplm_fe
pcdtestplm_fe test
vifols
kappamodel.matrixols exact
adf.testdflog_wage k
kpss.testdflog_wage null
plm_re plmlog_wage training age edu data pdata model
phtestplm_fe plm_re
resettestols power type
Decision table:
| Test | Null | Action if rejected |
|---|
shapiro.test / jarque.bera.test | residuals Normal | bootstrap CIs if N small |
bptest | homoskedastic | use HC3 via coeftest(ols, vcov = vcovHC(ols, "HC3")) or cluster |
dwtest / bgtest | no autocorr | HAC SEs (vcovHAC) or cluster by unit |
pbgtest (panel) | no panel autocorr | cluster by entity |
pcdtest | no CSD | Driscoll–Kraay (vcovDC) |
vif > 10 | — | drop / combine |
| ADF rejects + KPSS doesn't | stationary | levels |
| ADF doesn't reject | unit root | first-difference |
phtest | RE consistent | use FE |
Step 5 — Baseline empirical modeling (Section 4: Main Results)
Deeper patterns: references/05-modeling.md — every estimator. fixest is the workhorse.
This is the densest section of an applied paper. A modern AER §4 typically contains 2–3 multi-regression tables and one coefficient plot:
- Table 2 (main): progressive controls, 4–6 columns — Pattern A below
- Table 2-bis (design horse race): same coefficient under OLS / IV / DID / DML — Pattern B
- Table 2-ter (multi-outcome): same treatment, several outcomes side-by-side — Pattern C
- Figure 3 (coefplot): visual summary of β̂ and 95% CI across specs
Estimator routing (memorize this — getting it wrong silently produces nonsense):
- No FE / single low-card FE →
feols(y ~ X, data, cluster = ~i)
- High-dim FE →
feols(y ~ X | fe1 + fe2, data, cluster = ~i)
- Two-way cluster →
feols(..., cluster = ~ firm_id + year)
- 2SLS / IV →
feols(y ~ X | D ~ Z, data, cluster = ~ firm_id) (or AER::ivreg for diagnostics)
- DID / event-study →
feols(y ~ sunab(G, t) | i + t, data) (SA) · did::att_gt (CS) · didimputation::did_imputation (BJS)
Pick by identification strategy:
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
Canonical calls (the eight patterns A–H below are the AER table cookbook — modelsummary(...) and fixest::etable(...) are the two workhorses, equivalent to Stata outreg2/esttab and Python pf.etable/Stargazer).
5.A Pattern A — Progressive controls (the canonical Table 2)
Stable β̂ across columns ⇒ less concern that selection on observables is driving the estimate (Oster 2019 selection-stability logic; quantified in 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 <- feolslog_wage training age edu tenure firm_size worker_id year
data df cluster firm_id
m6 feolslog_wage training age edu tenure firm_size worker_id year industryyear
data df cluster firm_id
modelsummary
m1
m2
m3
m4
m5
m6
output
stars
gof_omit
coef_map
notes
modelsummarym1m2m3m4m5m6
output
AER convention: show ALL controls (and the intercept). Pass NEITHER keep = NOR coef_omit = so every parameter is visible. Use coef_map = c("training" = "Training") (single mapping) only when a focal-coefficient-only table is intentional (interaction-form heterogeneity, IV first-stage triplet); use coef_omit = "Intercept" only when you want to suppress the constant for paper aesthetics.
5.B Pattern B — Design horse race (Table 2-bis)
Show the same coefficient of interest under multiple identification strategies. This is the AER credibility move: convergent evidence across designs each making different identifying assumptions.
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 est_method
clustervars
psm matchittraining age edu tenure data df
method ratio
psm_lm lmlog_wage training age edu tenure
data match.datapsm weights weights
ebal weightittraining age edu tenure data df method
ebal_lm lmlog_wage training age edu tenure
data df weights ebalweights
modelsummary
ols
iv
aggtecs type
psm_lm
ebal_lm
output
stars
coef_map
gof_omit
notes
5.C Pattern C — Multi-outcome table (same X, several Y's)
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
notes
5.D Pattern D — Stacked Panel A / Panel B table
Same model family, two horizons (short-run / long-run) or two samples. Use gt::gt_group() to stack two modelsummary blocks with panel headers.
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 modelsummarypanelA output
tab_headertitle
ms_B modelsummarypanelB output
tab_headertitle
gt_groupms_A ms_B
gtsave
gt_groupms_A ms_B
gtsave
5.E Pattern E — IV reporting triplet (first-stage / reduced-form / 2SLS)
The textbook AER IV table presents the first stage, the reduced form, and the 2SLS in three columns so the reader can verify Wald-ratio = RF / FS.
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 =
coef_map
gof_map raw clean fmt
notes
IV triplet is intentionally focal: show only Z + endogenous regressor so the reader can eyeball the Wald ratio. Drop coef_map= only if a referee asks for the full coefficient list.
5.F Pattern F — Causal-orchestrator main via did::att_gt / synthdid / grf::causal_forest
For DID / SCM / matching / forest mains, the modern R estimator returns a self-contained estimate + automatic placebos / pre-trends / overlap diagnostics. Pipe into modelsummary via the auto-tidiers.
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"))
print(aggte(cs, type = "dynamic", min_e = -4, max_e = 4))
library(synthdid)
sdid_setup <- panel.matrices(df, unit time
outcome treatment
sdid_fit synthdid_estimatesdid_setupY sdid_setupN0 sdid_setupT0
printsummarysdid_fit
librarygrf
cf causal_forestX as.matrixdf
Y dflog_wage W dftraining num.trees
average_treatment_effectcf target.sample
test_calibrationcf
variable_importancecf
5.G Pattern G — Subgroup modelsummary (Table 3, see Step 7)
One column per subgroup. Detailed code in §Step 7 — Heterogeneity.
5.H Pattern H — Robustness master (Table A1, see Step 6)
Stack every robustness specification next to the baseline. Detailed code in §Step 6.
Canonical estimator commands (the underlying primitives)
library(fixest)
ols <- feols(log_wage ~ training + age + edu + tenure,
data = df, cluster = ~ firm_id)
summary(ols)
fe <- feols(log_wage ~ training + age + edu + tenure | worker_id + year,
data = df, cluster = ~ worker_id)
fe_mw <- feols(log_wage ~ training | worker_id + year,
data = df, cluster = ~ worker_id + firm_id)
fe_hd <- feols(log_wage ~ training | worker_id + industry^year
data df cluster firm_id
did22 feolslog_wage itreated post ref age edu
data df cluster worker_id
did22 feolslog_wage itreated post ref worker_id year
data df cluster worker_id
es feolslog_wage irel_time ref worker_id year
data df filterfirst_treat
cluster worker_id
iplotes
xlab
main
librarydid
cs att_gtyname tname idname
gname data df
control_group
est_method
clustervars
ggdidcs
sa feolslog_wage sunabfirst_treat year worker_id year
data df cluster worker_id
iplotsa sub.title
librarydidimputation
bjs did_imputationdata df yname gname
tname idname
horizon pretrends
cluster_var
librarysynthdid
sdid_setup synthdidpanel.matricesdf unit time
outcome treatment
sdid_fit synthdid_estimatesdid_setupY sdid_setupN0 sdid_setupT0
iv feolslog_wage age edu training draft_lottery z2
data df cluster firm_id
summaryiv stage
fitstativ ivf ivwald sargan
libraryAER
iv_aer ivreglog_wage training age edu
draft_lottery z2 age edu data df
summaryiv_aer vcov. sandwich diagnostics
libraryrdrobust; libraryrddensity
rd rdrobusty dfoutcome x dfrunning_var
kernel bwselect
summaryrd
rdploty dfoutcome x dfrunning_var
rddensityX dfrunning_var
logit feglmemployed training age edu firm_id year
data df family binomiallink
cluster firm_id
librarymarginaleffects
avg_slopeslogit variables
pois fepoiscitations training age firm_id year
data df cluster firm_id
Step 6 — Robustness battery
Deeper patterns: references/06-robustness.md — modelsummary for M1–M6; clubSandwich/fwildclusterboot; bacondecomp/HonestDiD/robomit; ri2 randomization inference.
library(modelsummary)
library(fixest)
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 feolslog_wage training age edu tenure worker_id year industryyear
data df cluster worker_id
modelsummary m1 m2 m3
m4 m5 m6
stars
gof_omit
coef_map
output
cl
fit feolslog_wage training worker_id year data df
cluster as.formulapaste0 cl
catcl coeffit sefit
libraryfwildclusterboot
boot boottestm4 param clustid
B seed
summaryboot
splits
df filterfemale
df filterfemale
df filterage
df filterage
sub_fits imapsplits feolslog_wage training worker_id year
data .x cluster worker_id
modelsummarysub_fits stars
df_placebo df
mutatefake_first first_treat
fake_post year fake_first
filteryear first_treat
feolslog_wage fake_post worker_id year
data df_placebo cluster worker_id
libraryri2
ri_out conduct_riformula log_wage training age edu
declaration randomizrdeclare_raN nrowdf
prob meandftraining
assignment
sharp_hypothesis
data df
sims
summaryri_out; plotri_out
librarybacondecomp
bacon_out baconlog_wage training
data df id_var time_var
ggplotbacon_out aesweight estimate color type geom_point
ggsave
libraryHonestDiD
honest_out createSensitivityResultsbetahat escoefficients
sigma vcoves
numPrePeriods numPostPeriods
Mbarvec seq by
createSensitivityPlothonest_out originalResults honest_outmainResult
ggsave
libraryrobomit
o_testy x
con
id time
data df R2max fitstatm6 beta
librarymodelsummary; libraryMatchIt; libraryWeightIt
base feolslog_wage training age edu tenure industry year
data df cluster firm_id
no99 feolslog_wage training age edu tenure industry year
data df filterwage quantilewage na.rm
cluster firm_id
balpan feolslog_wage training age edu tenure industry year
data df group_byworker_id
filtern_distinctyear n_distinctyear ungroup
cluster firm_id
dropearly feolslog_wage training age edu tenure industry year
data df filterfirst_treat cluster firm_id
wfe feolslog_wage training age edu tenure worker_id year
data df cluster firm_id
cl2way feolslog_wage training age edu tenure industry year
data df cluster firm_id year
logy feolswage training age edu tenure industry year
data df cluster firm_id
ihsy feolswage training age edu tenure industry year
data df cluster firm_id
m_psm matchittraining age edu tenure firm_size data df method
psm_lm lmlog_wage training age edu tenure data match.datam_psm weights weights
ebal_w weightittraining age edu tenure firm_size data df method
ebal_lm lmlog_wage training age edu tenure data df weights ebal_wweights
modelsummary
base
no99
balpan
dropearly
wfe
cl2way
logy
ihsy
psm_lm
ebal_lm
output
stars
coef_map
gof_omit
notes
libraryspecr; libraryggplot2
specs setupdata df
y
x
model
controls
subsets industry
results specrspecs
plotresults choices
ggsave width height
ggsave width height dpi
libraryHonestDiD
es_pre coefesgrep coefes
es_post coefesgrep coefes
honest_out createSensitivityResultsbetahat es_pre es_post
sigma vcoveses_pre es_post
es_pre es_post
numPrePeriods es_pre
numPostPeriods es_post
Mbarvec seq by
createSensitivityPlothonest_out originalResults honest_outmainResult
ggsave width height
libraryEValue
evalueRR lo hi
Step 7 — Further analysis
Deeper patterns: references/07-further-analysis.md — marginaleffects is the post-estimation workhorse; mediation::mediate for Imai mediation; lavaan for SEM; grf::causal_forest for CATE.
library(marginaleffects)
library(fixest)
het <- feols(log_wage ~ i(female, training, ref = 0) + age + edu | worker_id + year,
data = df, cluster = ~ worker_id)
summary(het)
iplot(het)
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 by
geom_hlineyintercept linetype
labsx y
ggsave width height
ddd feolslog_wage treated post high_exposure worker_id year
data df cluster firm_id
out_ladder
y
out_laddery feolsas.formulapastey
data df cluster worker_id
modelsummaryout_ladder stars
coef_map
output
librarymediation
med_M lmhours_worked training age edu data df
med_Y lmlog_wage training hours_worked age edu data df
med mediatemed_M med_Y treat mediator
boot sims
summarymed; plotmed
medsens medsensmed rho.by effect.type
plotmedsens
librarygrf
cf causal_forestX as.matrixdf selectage edu tenure firm_size
Y dflog_wage W dftraining
num.trees min.node.size
dftau_hat predictcfpredictions
variable_importancecf
average_treatment_effectcf target.sample
ggplotdf aestenure tau_hat
geom_smoothmethod se
labsx y
ggsave
librarysplines
dr feolslog_wage nstraining_hours df age edu worker_id year
data df cluster worker_id
plot_predictionsdr condition
Step 8 — Publication tables & figures
This step is mandatory — every analysis run produces all 5 required tables (T1–T5) and all 4 required figures (F1–F4) defined in the Default Output Spec at the top of this skill. Do not skip Step 8 because "the regression already ran". A coefficient without a table and a figure is not how applied economics communicates a result.
Deeper patterns: references/08-tables-plots.md — modelsummary is the modern default (LaTeX/Word/HTML/Excel from one call); kableExtra for further LaTeX styling; gt for HTML/Word; ggplot2 + iplot + ggpubr + cowplot + binsreg for figures.
library(modelsummary)
library(kableExtra)
library(gt)
library(fixest)
library(ggplot2)
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
gof_map
notes
output
modelsummarym1 m2 m3 m4 m5 m6
stars output
librarygtsummary
tbl1 df
selectlog_wage age edu tenure female training
tbl_summaryby training
statistic all_continuous
add_p add_difference add_n bold_labels
tbl1 as_kable_extraformat booktabs
kableExtrasave_kable
tbl1 as_flex_table
flextablesave_as_docxpath
ladder
y
laddery feolsas.formulapastey
data df cluster worker_id
modelsummaryladder
stars
coef_map
notes
output
het_specs
df
df filterfemale
df filterfemale
df filterage
df filterage
df filterindustry
het_models imaphet_specs
feolslog_wage training age edu tenure worker_id year
data .x cluster worker_id
modelsummaryhet_models
stars
coef_map
notes
output
rob
feolslog_wage training worker_id year data df
cluster worker_id
feolslog_wage training worker_id year data df
cluster firm_id
feolslog_wage training worker_id year data df
cluster worker_id firm_id
feolslog_wage training worker_id year
data df mutatelog_wage DescToolsWinsorizelog_wage
probs
na.rm
cluster worker_id
feolslog_wage training worker_id year
data df filterindustry
cluster worker_id
feolslog_wage fake_post worker_id year
data df filteryear first_treat
cluster worker_id
modelsummaryrob
stars
output
modelplotm1 m2 m3 m4 m5 m6
coef_map
conf_level
geom_vlinexintercept linetype alpha
labsx y
title
theme_classicbase_size
ggsave width height
ggsave width height dpi
pdf width height
iplotes
xlab
ylab
main
ref.line
dev.off
png width height res
iplotes
xlab
ylab
