Full-empirical-analysis-skill-R

by brycewang-stanford4.6k

Communication

Install

The skill file stays in the publisher's repository.

Paste the copied command in a terminal. It installs this SKILL.md into that agent.

npx skills add brycewang-stanford/Auto-Empirical-Research-Skills --skill 00.3-Full-empirical-analysis-skill_R -a cursor -y

Download skill skills/00.3-Full-empirical-analysis-skill_R/SKILL.md on GitHub

Related skills

View all

Classical end-to-end empirical analysis workflow in the modern tidyverse + econometrics R ecosystem — dplyr + tidyr + haven + fixest + sandwich + lmtest + clubSandwich + AER + ivreg + did + bacondecomp + HonestDiD + eventstudyr + rdrobust + rddensity + Synth + gsynth + synthdid + MatchIt + WeightIt + cobalt + ebal + grf + DoubleML + mediation + marginaleffects + modelsummary + kableExtra + gt + ggplot2 + ggpubr + cowplot + binsreg. **Defaults to economics empirical-paper style** (AER / QJE / AEJ) — every run produces a publication-ready output set with a multi-column regression table (M1→M6 progressive controls/FE) as the centerpiece, plus Table 1 (descriptives), mechanism / heterogeneity / robustness tables, and event-study + coefficient + trend figures. Covers the full 8-step R pipeline an applied economist runs on every paper — (1) data import & cleaning (read_dta/read_csv, naniar, janitor, validate-merges), (2) variable construction (mutate/across/winsorize/group_by + lag/lead with dplyr), (3) descriptive statistics & Table 1 (gtsummary, modelsummary::datasummary, tableone), (4) classical diagnostic tests (shapiro/jarque.bera.test/bptest/dwtest/bgtest/vif/adf.test/kpss.test/Hausman), (5) baseline modeling (fixest::feols, ivreg, did::att_gt, eventstudyr, sun_ab, did_imputation, synthdid, rdrobust, MatchIt, WeightIt, grf::causal_forest, DoubleML, mediation), (6) robustness battery (modelsummary stack, clubSandwich CRSE, fwildclusterboot, ri2, robomit Oster, bacondecomp, HonestDiD), (7) further analysis (interactions + marginaleffects, mediation::mediate, gsem via lavaan, dose-response splines, grf CATE), (8) publication-ready tables & figures (modelsummary, kableExtra, gt, stargazer, texreg, flextable to LaTeX/Word/HTML; ggplot2 + ggpubr + cowplot + binsreg + iplot for figures). **Also covers two parallel domain modes that share the same 8-step scaffolding** — **Mode A — Epidemiology / public health** (target-trial emulation, IPTW + g-formula + TMLE doubly-robust triplet via `WeightIt` / `gfoRmula` / `tmle` / `ltmle`, Mendelian randomization via `MendelianRandomization` / `TwoSampleMR` / `MRPRESSO`, KM / Cox / AFT / RMST survival via `survival` / `survminer` / `flexsurv`, E-value sensitivity via `EValue`, principal stratification — STROBE / TRIPOD reporting), and **Mode B — ML causal inference** (DML via `DoubleML`, S/T/X/R/DR meta-learners via `causalweight` / `grf`, causal forest via `grf::causal_forest`, BART/BCF via `bartCause` / `bcf`, matrix completion via `MCPanel`, CATE distribution + policy tree via `policytree`, off-policy evaluation, conformal causal via `conformalInference` / `cfcausal`, fairness audit via `fairmodels`, DAG learning via `pcalg` / `bnlearn` / LLM-assisted). Use when the user asks for a complete R empirical analysis, wants a tidyverse-style reproducible R script / Quarto workflow, prefers fixest over reghdfe, needs the R counterpart to StatsPAI / 00.1 / 00.2, or names a specific R step in isolation ("feols with cluster", "MatchIt nearest neighbor", "bacondecomp in R", "gtsummary table 1", "modelsummary to Word"). Mode A triggers on "target trial emulation R", "tmle ltmle", "MendelianRandomization", "TwoSampleMR", "MRPRESSO", "survival cox AFT", "STROBE R", "EValue R", "公共健康 R", "流行病学 R". Mode B triggers on "DoubleML R", "grf causal forest", "policytree", "bartCause bcf", "conformal causal R", "fairmodels", "pcalg NOTEARS", "因果机器学习 R".

name: Full-empirical-analysis-skill-R description: Classical end-to-end empirical analysis workflow in the modern tidyverse + econometrics R ecosystem — dplyr + tidyr + haven + fixest + sandwich + lmtest + clubSandwich + AER + ivreg + did + bacondecomp + HonestDiD + eventstudyr + rdrobust + rddensity + Synth + gsynth + synthdid + MatchIt + WeightIt + cobalt + ebal + grf + DoubleML + mediation + marginaleffects + modelsummary + kableExtra + gt + ggplot2 + ggpubr + cowplot + binsreg. Defaults to economics empirical-paper style (AER / QJE / AEJ) — every run produces a publication-ready output set with a multi-column regression table (M1→M6 progressive controls/FE) as the centerpiece, plus Table 1 (descriptives), mechanism / heterogeneity / robustness tables, and event-study + coefficient + trend figures. Covers the full 8-step R pipeline an applied economist runs on every paper — (1) data import & cleaning (read_dta/read_csv, naniar, janitor, validate-merges), (2) variable construction (mutate/across/winsorize/group_by + lag/lead with dplyr), (3) descriptive statistics & Table 1 (gtsummary, modelsummary::datasummary, tableone), (4) classical diagnostic tests (shapiro/jarque.bera.test/bptest/dwtest/bgtest/vif/adf.test/kpss.test/Hausman), (5) baseline modeling (fixest::feols, ivreg, did::att_gt, eventstudyr, sun_ab, did_imputation, synthdid, rdrobust, MatchIt, WeightIt, grf::causal_forest, DoubleML, mediation), (6) robustness battery (modelsummary stack, clubSandwich CRSE, fwildclusterboot, ri2, robomit Oster, bacondecomp, HonestDiD), (7) further analysis (interactions + marginaleffects, mediation::mediate, gsem via lavaan, dose-response splines, grf CATE), (8) publication-ready tables & figures (modelsummary, kableExtra, gt, stargazer, texreg, flextable to LaTeX/Word/HTML; ggplot2 + ggpubr + cowplot + binsreg + iplot for figures). Also covers two parallel domain modes that share the same 8-step scaffolding — Mode A — Epidemiology / public health (target-trial emulation, IPTW + g-formula + TMLE doubly-robust triplet via WeightIt / gfoRmula / tmle / ltmle, Mendelian randomization via / / , KM / Cox / AFT / RMST survival via / / , E-value sensitivity via , principal stratification — STROBE / TRIPOD reporting), and (DML via , S/T/X/R/DR meta-learners via / , causal forest via , BART/BCF via / , matrix completion via , CATE distribution + policy tree via , off-policy evaluation, conformal causal via / , fairness audit via , DAG learning via / / LLM-assisted). Use when the user asks for a complete R empirical analysis, wants a tidyverse-style reproducible R script / Quarto workflow, prefers fixest over reghdfe, needs the R counterpart to StatsPAI / 00.1 / 00.2, or names a specific R step in isolation ("feols with cluster", "MatchIt nearest neighbor", "bacondecomp in R", "gtsummary table 1", "modelsummary to Word"). Mode A triggers on "target trial emulation R", "tmle ltmle", "MendelianRandomization", "TwoSampleMR", "MRPRESSO", "survival cox AFT", "STROBE R", "EValue R", "公共健康 R", "流行病学 R". Mode B triggers on "DoubleML R", "grf causal forest", "policytree", "bartCause bcf", "conformal causal R", "fairmodels", "pcalg NOTEARS", "因果机器学习 R". triggers:

library(validate); library(assertr) data_contract <- function(df, y, treatment, id = NULL, time = NULL, covariates = c()) { keys <- c(y, treatment, id, time, covariates) contract <- list( n_obs = nrow(df), # 1. shape dtypes = sapply(df[keys], function(x) class(x)[1]), # 2. dtypes n_missing = sapply(df[keys], function(x) sum(is.na(x))), # 3. missingness n_dupes_on_keys = if (!is.null(id) && !is.null(time)) sum(duplicated(df[, c(id, time)])) else 0, # 4. duplicates panel_balanced = NULL, cohort_sizes = NULL ) if (!is.null(id) && !is.null(time)) { bal <- df %>% count(.data[[id]]) contract$panel_balanced <- all(bal$n == max(bal$n)) # 5. balance contract$n_dropped_by_balance <- sum(bal$n != max(bal$n)) if ("first_treat" %in% names(df)) { contract$cohort_sizes <- df %>% distinct(.data[[id]], .keep_all = TRUE) %>% count(first_treat) %>% deframe() } } contract$y_range <- range(df[[y]], na.rm = TRUE) contract$treatment_share <- mean(df[[treatment]], na.rm = TRUE) # MCAR sniff test (Rubin) — if missing(y) is associated with covariates, # listwise deletion biases the estimate. Use mice / IPW instead. miss_y <- is.na(df[[y]]) contract$mcar_hint <- "likely MCAR (listwise OK)" if (any(miss_y) && any(!miss_y)) { for (cov in covariates) { if (is.numeric(df[[cov]])) { p <- t.test(df[[cov]][miss_y], df[[cov]][!miss_y])$p.value if (p < 0.05) { contract$mcar_hint <- sprintf("NOT MCAR (y-miss differs on %s, p=%.3f) → use mice / IPW", cov, p) break } } } } contract } contract <- data_contract(df, y = "wage", treatment = "training", id = "worker_id", time = "year", covariates = c("age", "edu", "tenure")) stopifnot(contract$n_dupes_on_keys == 0) stopifnot(all(contract$n_missing == 0)) write_json(contract, "artifacts/data_contract.json", pretty = TRUE, auto_unbox = TRUE)

library(gtsummary) library(modelsummary) # 3a. Full-sample summary — one line, publication ready df %>% select(log_wage, age, edu, tenure, training) %>% datasummary_skim() # Or df %>% select(log_wage, age, edu, tenure, training) %>% tbl_summary( type = list(all_continuous() ~ "continuous2"), statistic = all_continuous() ~ c("{N_nonmiss}", "{mean} ({sd})", "{min} – {median} – {max}") ) %>% bold_labels() %>% as_kable_extra() %>% kableExtra::save_kable("tables/table1_full.tex") # 3b. Stratified Table 1 (treated vs control, with SMDs + p-values) df %>% select(log_wage, age, edu, tenure, female, training) %>% tbl_summary(by = training, missing = "ifany") %>% add_p() %>% add_difference() %>% add_n() %>% modify_header(label = "**Variable**") %>% bold_labels() %>% as_gt() %>% gt::gtsave("tables/table1_balance.html") # Or via modelsummary (writes LaTeX/Word/HTML) datasummary_balance(~ training, data = df %>% select(training, age, edu, tenure, female), output = "tables/table1_balance.tex") # 3c. Correlation matrix with stars library(corrplot); library(psych) corr_obj <- corr.test(df %>% select(log_wage, age, edu, tenure, training), method = "pearson") corrplot(corr_obj$r, method = "color", type = "upper", p.mat = corr_obj$p, sig.level = 0.05, insig = "blank", addCoef.col = "black", number.cex = 0.7, tl.col = "black", tl.srt = 45, col = colorRampPalette(c("#B2182B","white","#2166AC"))(200)) # 3d. Distribution plots library(ggplot2) p1 <- ggplot(df, aes(log_wage, fill = factor(training))) + geom_density(alpha = 0.5) + scale_fill_manual(values = c("0" = "darkred", "1" = "navy"), labels = c("Control", "Treated"), name = "") + labs(x = "Log wage", y = "Density", title = "Log-wage density by treatment") + theme_classic() p2 <- ggplot(df, aes(sample = log_wage)) + stat_qq() + stat_qq_line() + labs(title = "Normal Q-Q") + theme_classic() cowplot::plot_grid(p1, p2, labels = "auto") %>% ggsave("figures/distributions.pdf", plot = ., width = 10, height = 4) # 3e. Time-trend (DID motivation) df %>% group_by(year, training) %>% summarise(mean_log_wage = mean(log_wage, na.rm = TRUE), .groups = "drop") %>% ggplot(aes(year, mean_log_wage, color = factor(training))) + geom_line(linewidth = 1) + geom_point(size = 2) + geom_vline(xintercept = policy_year, linetype = "dashed") + scale_color_manual(values = c("0" = "darkred", "1" = "navy"), labels = c("Control","Treated"), name = "") + labs(x = "Year", y = "Mean log wage") + theme_classic() ggsave("figures/trend_did.pdf", width = 7, height = 4)

library(fixest) # 5a. OLS with cluster-robust SEs — feols is the modern primary ols <- feols(log_wage ~ training + age + edu + tenure, data = df, cluster = ~ firm_id) summary(ols) # 5b. Two-way FE — single line fe <- feols(log_wage ~ training + age + edu + tenure | worker_id + year, data = df, cluster = ~ worker_id) # Multi-way clustering fe_mw <- feols(log_wage ~ training | worker_id + year, data = df, cluster = ~ worker_id + firm_id) # High-dim interaction FE fe_hd <- feols(log_wage ~ training | worker_id + industry^year, data = df, cluster = ~ firm_id) # 5c. 2×2 DID did22 <- feols(log_wage ~ i(treated, post, ref = 0) + age + edu, data = df, cluster = ~ worker_id) # Or with absorbed FE: did22 <- feols(log_wage ~ i(treated, post, ref = 0) | worker_id + year, data = df, cluster = ~ worker_id) # 5d. Event study — base period at -1 es <- feols(log_wage ~ i(rel_time, ref = -1) | worker_id + year, data = df %>% filter(!is.na(first_treat)), cluster = ~ worker_id) iplot(es, xlab = "Years relative to treatment", main = "Event study") # 5e. Staggered DID — modern estimators (see references/05-modeling.md §5.4) library(did) cs <- att_gt(yname = "log_wage", tname = "year", idname = "worker_id", gname = "first_treat", data = df, control_group = "nevertreated", est_method = "dr", clustervars = "firm_id") ggdid(cs) # event-study plot # Sun & Abraham via fixest::sunab sa <- feols(log_wage ~ sunab(first_treat, year) | worker_id + year, data = df, cluster = ~ worker_id) iplot(sa, sub.title = "Sun-Abraham (2021)") # Borusyak–Jaravel–Spiess (didimputation) library(didimputation) bjs <- did_imputation(data = df, yname = "log_wage", gname = "first_treat", tname = "year", idname = "worker_id", horizon = 0:5, pretrends = -5:-1, cluster_var = "worker_id") # Synthetic DID library(synthdid) sdid_setup <- synthdid::panel.matrices(df, unit = "worker_id", time = "year", outcome = "log_wage", treatment = "training") sdid_fit <- synthdid_estimate(sdid_setup$Y, sdid_setup$N0, sdid_setup$T0) # 5f. IV / 2SLS iv <- feols(log_wage ~ age + edu | training ~ draft_lottery + z2, data = df, cluster = ~ firm_id) summary(iv, stage = 1) fitstat(iv, ~ ivf + ivwald + sargan) # first-stage F + Wald + overid # Or via AER: library(AER) iv_aer <- ivreg(log_wage ~ training + age + edu | draft_lottery + z2 + age + edu, data = df) summary(iv_aer, vcov. = sandwich, diagnostics = TRUE) # 5g. Sharp RD library(rdrobust); library(rddensity) rd <- rdrobust(y = df$outcome, x = df$running_var, c = 0, kernel = "triangular", bwselect = "mserd") summary(rd) rdplot(y = df$outcome, x = df$running_var, c = 0) rddensity(X = df$running_var, c = 0) # manipulation test # 5h. Binary outcome logit <- feglm(employed ~ training + age + edu | firm_id + year, data = df, family = binomial(link = "logit"), cluster = ~ firm_id) library(marginaleffects) avg_slopes(logit, variables = "training") # AME # 5i. Count w/ HD FE pois <- fepois(citations ~ training + age | firm_id + year, data = df, cluster = ~ firm_id)

library(modelsummary) library(fixest) # 6a. Progressive specs (M1 → M6) m1 <- feols(log_wage ~ training, data = df, cluster = ~ firm_id) m2 <- feols(log_wage ~ training + age + edu, data = df, cluster = ~ firm_id) m3 <- feols(log_wage ~ training + age + edu + tenure | worker_id, data = df, cluster = ~ worker_id) m4 <- feols(log_wage ~ training + age + edu + tenure | worker_id + year, data = df, cluster = ~ worker_id) m5 <- feols(log_wage ~ training + age + edu + tenure | worker_id + year + region, data = df, cluster = ~ worker_id) m6 <- feols(log_wage ~ training + age + edu + tenure | worker_id + year + industry^year, data = df, cluster = ~ worker_id) modelsummary(list("(1)" = m1, "(2)" = m2, "(3)" = m3, "(4)" = m4, "(5)" = m5, "(6)" = m6), stars = c('*' = .1, '**' = .05, '***' = .01), gof_omit = "BIC|AIC|F|Log", coef_map = c("training" = "Training", "age" = "Age", "edu" = "Education", "tenure" = "Tenure"), output = "tables/table_main.tex") # 6b. Alternative cluster levels for (cl in c("worker_id", "firm_id", "industry", "state")) { fit <- feols(log_wage ~ training | worker_id + year, data = df, cluster = as.formula(paste0("~", cl))) cat(cl, ": b=", coef(fit)["training"], " se=", se(fit)["training"], "\n") } # 6c. Wild cluster bootstrap (when few clusters) library(fwildclusterboot) boot <- boottest(m4, param = "training", clustid = "state", B = 9999, seed = 42) summary(boot) # 6d. Subsample splits splits <- list( "Female=0" = df %>% filter(female == 0), "Female=1" = df %>% filter(female == 1), "Young (<40)" = df %>% filter(age < 40), "Old (>=40)" = df %>% filter(age >= 40) ) sub_fits <- imap(splits, ~ feols(log_wage ~ training | worker_id + year, data = .x, cluster = ~ worker_id)) modelsummary(sub_fits, stars = TRUE) # 6e. Placebo — fake timing df_placebo <- df %>% mutate(fake_first = first_treat - 3, fake_post = year >= fake_first) %>% filter(year < first_treat) feols(log_wage ~ fake_post | worker_id + year, data = df_placebo, cluster = ~ worker_id) # 6f. Randomization inference library(ri2) ri_out <- conduct_ri(formula = log_wage ~ training + age + edu, declaration = randomizr::declare_ra(N = nrow(df), prob = mean(df$training)), assignment = "training", sharp_hypothesis = 0, data = df, sims = 1000) summary(ri_out); plot(ri_out) # 6g. TWFE bias diagnosis library(bacondecomp) bacon_out <- bacon(log_wage ~ training, data = df, id_var = "worker_id", time_var = "year") ggplot(bacon_out, aes(weight, estimate, color = type)) + geom_point() ggsave("figures/bacon.pdf") # 6h. Parallel-trends sensitivity library(HonestDiD) honest_out <- createSensitivityResults(betahat = es$coefficients, sigma = vcov(es), numPrePeriods = 5, numPostPeriods = 5, Mbarvec = seq(0, 0.5, by = 0.05)) createSensitivityPlot(honest_out, originalResults = honest_out$mainResult) ggsave("figures/honestdid.pdf") # 6i. Oster (2019) δ* library(robomit) o_test(y = "log_wage", x = "training", con = "age + edu + tenure | worker_id + year", id = "worker_id", time = "year", data = df, R2max = 1.3 * fitstat(m6, "r2"), beta = 0) # ============================================================ # 6j. Pattern H — Robustness master table (Table A1, one column per check) # ============================================================ library(modelsummary); library(MatchIt); library(WeightIt) base <- feols(log_wage ~ training + age + edu + tenure | industry + year, data = df, cluster = ~ firm_id) no99 <- feols(log_wage ~ training + age + edu + tenure | industry + year, data = df %>% filter(wage < quantile(wage, 0.99, na.rm = TRUE)), cluster = ~ firm_id) balpan <- feols(log_wage ~ training + age + edu + tenure | industry + year, data = df %>% group_by(worker_id) %>% filter(n_distinct(year) == max(n_distinct(year))) %>% ungroup(), cluster = ~ firm_id) dropearly <- feols(log_wage ~ training + age + edu + tenure | industry + year, data = df %>% filter(first_treat > 2008), cluster = ~ firm_id) wfe <- feols(log_wage ~ training + age + edu + tenure | worker_id + year, data = df, cluster = ~ firm_id) cl2way <- feols(log_wage ~ training + age + edu + tenure | industry + year, data = df, cluster = ~ firm_id + year) logy <- feols(log(wage + 1) ~ training + age + edu + tenure | industry + year, data = df, cluster = ~ firm_id) ihsy <- feols(asinh(wage) ~ training + age + edu + tenure | industry + year, data = df, cluster = ~ firm_id) m_psm <- matchit(training ~ age + edu + tenure + firm_size, data = df, method = "nearest") psm_lm <- lm(log_wage ~ training + age + edu + tenure, data = match.data(m_psm), weights = weights) ebal_w <- weightit(training ~ age + edu + tenure + firm_size, data = df, method = "ebal") ebal_lm <- lm(log_wage ~ training + age + edu + tenure, data = df, weights = ebal_w$weights) modelsummary( list("(1) Baseline" = base, "(2) Drop top 1%" = no99, "(3) Balanced" = balpan, "(4) Drop early" = dropearly, "(5) Worker FE" = wfe, "(6) 2-way cluster" = cl2way, "(7) log Y" = logy, "(8) IHS Y" = ihsy, "(9) PSM" = psm_lm, "(10) Entropy bal." = ebal_lm), output = "tables/tableA1_robustness.tex", stars = c("*" = 0.1, "**" = 0.05, "***" = 0.01), coef_map = c("training" = "Training (β̂)"), gof_omit = "BIC|AIC|F|Log|Adj", notes = "Each column is one robustness check. β̂ on training is the focal coefficient." ) # ============================================================ # 6k. Specification curve (Simonsohn-Simmons-Nelson 2020) via `specr` # ============================================================ library(specr); library(ggplot2) specs <- setup(data = df, y = c("log_wage", "ihs_wage"), x = "training", model = c("feols"), controls = c("age", "edu", "tenure", "firm_size"), subsets = list(industry = c("manuf", "construction", "transport"))) results <- specr(specs) plot(results, choices = c("x", "y", "controls", "subsets")) ggsave("figures/fig5_spec_curve.pdf", width = 10, height = 6) ggsave("figures/fig5_spec_curve.png", width = 10, height = 6, dpi = 300) # Hand-rolled alternative when `specr` doesn't fit (with custom FE / SE): # spec_grid <- expand.grid(controls = list(c("age"), c("age","edu"), c("age","edu","tenure")), # ytrans = c("log_wage", "ihs_wage"), # sample = c("all", "manuf", "no99"), # cluster = c("firm_id", "firm_id+year")) # Loop, run feols, collect b/se, ggplot::geom_pointrange. # ============================================================ # 6l. Sensitivity dashboard — HonestDiD + Oster + E-value # ============================================================ # (a) HonestDiD — Rambachan-Roth (2023): bound on β̂ under bounded PT violation library(HonestDiD) es_pre <- coef(es)[grep("year::-", names(coef(es)))] es_post <- coef(es)[grep("year::[0-9]", names(coef(es)))] honest_out <- createSensitivityResults(betahat = c(es_pre, es_post), sigma = vcov(es)[c(names(es_pre), names(es_post)), c(names(es_pre), names(es_post))], numPrePeriods = length(es_pre), numPostPeriods = length(es_post), Mbarvec = seq(0, 0.5, by = 0.05)) createSensitivityPlot(honest_out, originalResults = honest_out$mainResult) ggsave("figures/fig6_honestdid.pdf", width = 7, height = 4) # (b) Oster δ — `robomit::o_test` (already shown in 6i) # (c) E-value (VanderWeele-Ding 2017) — for risk-ratio outcomes library(EValue) evalue(RR(1.45), lo = 1.10, hi = 1.91) # → reports the minimum strength of unmeasured confounding to nullify the result.

library(marginaleffects) library(fixest) # 7a. Heterogeneity via interaction het <- feols(log_wage ~ i(female, training, ref = 0) + age + edu | worker_id + year, data = df, cluster = ~ worker_id) summary(het) iplot(het) # visualize interaction # Continuous moderator + marginsplot het_c <- feols(log_wage ~ training * tenure + age + edu | worker_id + year, data = df, cluster = ~ worker_id) plot_slopes(het_c, variables = "training", condition = list(tenure = seq(0, 20, by = 2))) + geom_hline(yintercept = 0, linetype = "dashed") + labs(x = "Tenure", y = "Marginal effect of training") ggsave("figures/het_tenure.pdf", width = 6, height = 4) # 7b. Triple difference ddd <- feols(log_wage ~ treated * post * high_exposure | worker_id + year, data = df, cluster = ~ firm_id) # 7c. Outcome ladder out_ladder <- list() for (y in c("hours_worked", "productivity", "log_wage")) { out_ladder[[y]] <- feols(as.formula(paste(y, "~ training | worker_id + year")), data = df, cluster = ~ worker_id) } modelsummary(out_ladder, stars = TRUE, coef_map = c("training" = "Training"), output = "tables/outcome_ladder.tex") # 7d. Mediation — Imai et al. (2010) library(mediation) med_M <- lm(hours_worked ~ training + age + edu, data = df) med_Y <- lm(log_wage ~ training + hours_worked + age + edu, data = df) med <- mediate(med_M, med_Y, treat = "training", mediator = "hours_worked", boot = TRUE, sims = 1000) summary(med); plot(med) # Sensitivity to unobserved M-Y confounding medsens <- medsens(med, rho.by = 0.05, effect.type = "indirect") plot(medsens) # 7e. CATE via causal forest library(grf) cf <- causal_forest(X = as.matrix(df %>% select(age, edu, tenure, firm_size)), Y = df$log_wage, W = df$training, num.trees = 2000, min.node.size = 5) df$tau_hat <- predict(cf)$predictions variable_importance(cf) average_treatment_effect(cf, target.sample = "all") # Plot CATE by a moderator ggplot(df, aes(tenure, tau_hat)) + geom_smooth(method = "loess", se = TRUE) + labs(x = "Tenure", y = "Estimated CATE") ggsave("figures/cate_tenure.pdf") # 7f. Dose-response — splines library(splines) dr <- feols(log_wage ~ ns(training_hours, df = 4) + age + edu | worker_id + year, data = df, cluster = ~ worker_id) plot_predictions(dr, condition = "training_hours")

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

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

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

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

§	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`

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

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`

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

Install

Related skills

Questions

Mode A — Epidemiology / public health

Mode B — ML causal inference

Full Empirical Analysis — Classical R Workflow

Philosophy

SkillOpt-style execution gate

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

Default Output Spec — Economics Empirical Paper

Required tables (always produced)

Required figures (always produced)

Output file layout (default)

When to deviate

Required packages

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

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

Export cookbook — LaTeX / Word / Excel in one block

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

Step 0 — Sample-construction log & 5-check data contract

0.1 Sample-construction log (footnote 4)

0.2 Five-check data contract (go / no-go gate)

Step 1 — Data import & cleaning

Step 2 — Variable construction & transformation

Step 2.5 — Empirical strategy (write the equation + identifying assumption)

Equation × identifying assumption × R estimator (decision table)

Design picker (when the user is unsure)

Pre-registration strategy.md template

Step 3 — Descriptive statistics & Table 1

Step 3.5 — Identification graphics (Section "Identification, graphical evidence")

3.5.1 Event-study figure + numerical pre-trends test (DID identification)

3.5.2 First-stage F-statistic + scatter (IV identification)

3.5.3 RD: McCrary density + canonical RD plot

3.5.4 Matching: love plot (standardized differences pre vs post)

3.5.5 SCM: synthetic-control trajectory + gap plot

Step 4 — Diagnostic statistical tests

Step 5 — Baseline empirical modeling (Section 4: Main Results)

5.A Pattern A — Progressive controls (the canonical Table 2)

5.B Pattern B — Design horse race (Table 2-bis)

5.C Pattern C — Multi-outcome table (same X, several Y's)

5.D Pattern D — Stacked Panel A / Panel B table

5.E Pattern E — IV reporting triplet (first-stage / reduced-form / 2SLS)

5.F Pattern F — Causal-orchestrator main via did::att_gt / synthdid / grf::causal_forest

5.G Pattern G — Subgroup modelsummary (Table 3, see Step 7)

5.H Pattern H — Robustness master (Table A1, see Step 6)

Canonical estimator commands (the underlying primitives)

Step 6 — Robustness battery

Step 7 — Further analysis

Step 8 — Publication tables & figures

Pre-registration `strategy.md` template

5.F Pattern F — Causal-orchestrator main via `did::att_gt` / `synthdid` / `grf::causal_forest`

5.G Pattern G — Subgroup `modelsummary` (Table 3, see Step 7)