full-empirical-analysis-skill-r

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

Aller à l'installation

Dépôt: brycewang-stanford/Auto-Empirical-Research-Skills
Dernière activité de la source: 22 juin 2026 à 09:41
Langue détectée de SKILL.md: anglais
Étoiles: 3 291
Forks: 432

Options d'installation

Le prompt qui vérifie d'abord la source est sélectionné par défaut. Vous pouvez passer à une commande directe ou télécharger une copie locale.

Installer avec Codex ou Claude Copiez ce prompt, collez-le dans Codex, Claude ou un autre assistant, puis laissez-le vérifier la page du skill et l'installer pour vous.

Une commande directe contourne le prompt de vérification. Examinez la source avant de l'exécuter.

npx skills add https://github.com/brycewang-stanford/Auto-Empirical-Research-Skills --skill full-empirical-analysis-skill-r

La commande reste sur une seule ligne. Faites défiler horizontalement pour la vérifier avant de la copier.

Vous préférez une copie locale ? Téléchargez les fichiers actuellement disponibles dans SkillsMP.

Vérifiez les fichiers source

Lisez SKILL.md et les fichiers associés affichés par SkillsMP avant de décider de l'installer.

Explorateur de fichiers

10 fichiers

SKILL.md

Instructions source · Aperçu en lecture seule

Retour à la recherche Retour à l’installation

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.

Full Empirical Analysis — Classical R Workflow

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.

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")

library(modelsummary) library(kableExtra) library(gt) library(fixest) library(ggplot2) # ============================================================ # 8a. ★ TABLE 2 — Main results, multi-column regression M1→M6 # (the centerpiece of every economics paper) # ============================================================ modelsummary( list("(1) Raw" = m1, "(2) +Demog" = m2, "(3) +Tenure" = m3, "(4) +Unit FE" = m4, "(5) +2-way FE" = m5, "(6) +Ind×Yr FE" = m6), stars = c('*' = .1, '**' = .05, '***' = .01), coef_map = c("training" = "Training", "age" = "Age", "edu" = "Education", "tenure" = "Tenure"), gof_map = list( list("raw" = "nobs", "clean" = "N", "fmt" = 0), list("raw" = "r.squared", "clean" = "R²", "fmt" = 3), list("raw" = "adj.r.squared","clean" = "Adj. R²", "fmt" = 3) ), notes = "Cluster-robust SE at worker_id in parentheses. * p<0.10, ** p<0.05, *** p<0.01.", output = "tables/table2_main.tex" ) modelsummary(list("(1)"=m1, "(2)"=m2, "(3)"=m3, "(4)"=m4, "(5)"=m5, "(6)"=m6), stars = TRUE, output = "tables/table2_main.docx") # ============================================================ # 8b. TABLE 1 — Summary statistics & balance # ============================================================ library(gtsummary) tbl1 <- df %>% select(log_wage, age, edu, tenure, female, training) %>% tbl_summary(by = training, missing = "ifany", statistic = all_continuous() ~ "{mean} ({sd})") %>% add_p() %>% add_difference() %>% add_n() %>% bold_labels() tbl1 %>% as_kable_extra(format = "latex", booktabs = TRUE) %>% kableExtra::save_kable("tables/table1_balance.tex") tbl1 %>% as_flex_table() %>% flextable::save_as_docx(path = "tables/table1_balance.docx") # ============================================================ # 8c. TABLE 3 — Mechanism / outcome ladder (3+ outcomes) # ============================================================ ladder <- list() for (y in c("hours_worked", "productivity", "log_wage")) { ladder[[y]] <- feols(as.formula(paste(y, "~ training + age + edu + tenure | worker_id + year")), data = df, cluster = ~ worker_id) } modelsummary(ladder, stars = c('*' = .1, '**' = .05, '***' = .01), coef_map = c("training" = "Training"), notes = "Each column is a separate regression on the labelled outcome. Cluster-robust SE at worker_id.", output = "tables/table3_mechanism.tex") # ============================================================ # 8d. TABLE 4 — Heterogeneity (subgroup × main coef) # ============================================================ het_specs <- list( "All" = df, "Female=0" = df %>% filter(female == 0), "Female=1" = df %>% filter(female == 1), "Age<40" = df %>% filter(age < 40), "Age≥40" = df %>% filter(age >= 40), "Manuf." = df %>% filter(industry == "manufacturing") ) het_models <- imap(het_specs, ~ feols(log_wage ~ training + age + edu + tenure | worker_id + year, data = .x, cluster = ~ worker_id)) modelsummary(het_models, stars = c('*' = .1, '**' = .05, '***' = .01), coef_map = c("training" = "Training"), notes = "Cluster-robust SE at worker_id. Wald p-values for cross-subgroup equality should accompany this table — see references/07.", output = "tables/table4_heterogeneity.tex") # ============================================================ # 8e. TABLE 5 — Robustness battery (alt SE / cluster / sample / placebo) # ============================================================ rob <- list( "Baseline" = feols(log_wage ~ training | worker_id + year, data = df, cluster = ~ worker_id), "Cluster=Firm" = feols(log_wage ~ training | worker_id + year, data = df, cluster = ~ firm_id), "2-way Cluster" = feols(log_wage ~ training | worker_id + year, data = df, cluster = ~ worker_id + firm_id), "Winsor 1/99" = feols(log_wage ~ training | worker_id + year, data = df %>% mutate(log_wage = DescTools::Winsorize(log_wage, probs = c(.01,.99), na.rm = TRUE)), cluster = ~ worker_id), "Drop Manuf." = feols(log_wage ~ training | worker_id + year, data = df %>% filter(industry != "manufacturing"), cluster = ~ worker_id), "Placebo (-3)" = feols(log_wage ~ fake_post | worker_id + year, data = df %>% filter(year < first_treat), cluster = ~ worker_id) ) modelsummary(rob, stars = c('*' = .1, '**' = .05, '***' = .01), output = "tables/table5_robustness.tex") # ============================================================ # 8f. ★ FIGURE 3 — Coefficient plot across M1→M6 # ============================================================ modelplot(list("(1)"=m1, "(2)"=m2, "(3)"=m3, "(4)"=m4, "(5)"=m5, "(6)"=m6), coef_map = c("training" = "Training"), conf_level = 0.95) + geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) + labs(x = "Coefficient on training (95% CI)", y = "Specification", title = "Effect of training across specifications") + theme_classic(base_size = 11) ggsave("figures/fig3_coefplot.pdf", width = 6, height = 4) ggsave("figures/fig3_coefplot.png", width = 6, height = 4, dpi = 300) # ============================================================ # 8g. FIGURE 2 — Event-study plot (dynamic DID, base period = -1) # ============================================================ pdf("figures/fig2_event_study.pdf", width = 7, height = 4) iplot(es, xlab = "Years relative to treatment", ylab = "Coefficient (ATT, 95% CI)", main = "Event study: dynamic effect of training", ref.line = -0.5) dev.off() png("figures/fig2_event_study.png", width = 2100, height = 1200, res = 300) iplot(es, xlab = "Years relative to treatment", ylab = "Coefficient (ATT, 95% CI)",

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 `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".
triggers	["R empirical analysis","tidyverse econometrics workflow","reproducible R script","Quarto empirical pipeline","fixest feols feglm fepois","high-dimensional fixed effects R","clubSandwich cluster-robust","fwildclusterboot wild cluster bootstrap","ivreg AER 2SLS R","did att_gt Callaway SantAnna R","eventstudyr event study R","did_imputation Borusyak R","synthdid R package","bacondecomp R Goodman Bacon","HonestDiD R Rambachan Roth","rdrobust R","rddensity R","Synth gsynth R","MatchIt nearest neighbor R","WeightIt IPW propensity R","cobalt balance check R","ebal entropy balancing R","grf causal forest R","DoubleML R","mediation R Imai","marginaleffects R","gtsummary table 1","modelsummary publication table","kableExtra LaTeX","texreg stargazer","flextable Word","ggplot2 coefplot","iplot fixest","binsreg R","haven read_dta sav","janitor clean_names","naniar missing","epidemiology pipeline R","public health causal inference R","target trial emulation R","g-formula R gfoRmula","IPTW marginal structural model R","WeightIt PSweight","tmle ltmle doubly robust","HAL-TMLE R","Mendelian randomization R","MendelianRandomization package","TwoSampleMR","MRPRESSO","MR-Egger weighted median R","STROBE TRIPOD reporting R","EValue sensitivity R","Kaplan-Meier AFT survival R","survival survminer flexsurv","流行病学 R","公共健康 R","ML causal inference R","DoubleML R","grf causal forest R","meta-learner S T X R DR R","causalweight R","bartCause bcf","Bayesian causal forest BCF R","CATE distribution R","policytree R","off-policy evaluation R","conformalInference cfcausal","conformal causal prediction R","fairmodels fairness audit","causal discovery PC NOTEARS R","pcalg bnlearn","因果机器学习 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 `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".
triggers	["R empirical analysis","tidyverse econometrics workflow","reproducible R script","Quarto empirical pipeline","fixest feols feglm fepois","high-dimensional fixed effects R","clubSandwich cluster-robust","fwildclusterboot wild cluster bootstrap","ivreg AER 2SLS R","did att_gt Callaway SantAnna R","eventstudyr event study R","did_imputation Borusyak R","synthdid R package","bacondecomp R Goodman Bacon","HonestDiD R Rambachan Roth","rdrobust R","rddensity R","Synth gsynth R","MatchIt nearest neighbor R","WeightIt IPW propensity R","cobalt balance check R","ebal entropy balancing R","grf causal forest R","DoubleML R","mediation R Imai","marginaleffects R","gtsummary table 1","modelsummary publication table","kableExtra LaTeX","texreg stargazer","flextable Word","ggplot2 coefplot","iplot fixest","binsreg R","haven read_dta sav","janitor clean_names","naniar missing","epidemiology pipeline R","public health causal inference R","target trial emulation R","g-formula R gfoRmula","IPTW marginal structural model R","WeightIt PSweight","tmle ltmle doubly robust","HAL-TMLE R","Mendelian randomization R","MendelianRandomization package","TwoSampleMR","MRPRESSO","MR-Egger weighted median R","STROBE TRIPOD reporting R","EValue sensitivity R","Kaplan-Meier AFT survival R","survival survminer flexsurv","流行病学 R","公共健康 R","ML causal inference R","DoubleML R","grf causal forest R","meta-learner S T X R DR R","causalweight R","bartCause bcf","Bayesian causal forest BCF R","CATE distribution R","policytree R","off-policy evaluation R","conformalInference cfcausal","conformal causal prediction R","fairmodels fairness audit","causal discovery PC NOTEARS R","pcalg bnlearn","因果机器学习 R"]

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

full-empirical-analysis-skill-r

Informations de source

Options d'installation

Vérifiez les fichiers source

Plus depuis ce dépôt

Plus depuis ce dépôt

Full Empirical Analysis — Classical R Workflow

Philosophy

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)