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

설치로 이동

저장소: caixinyu2017-star/xixi
최근 소스 활동: 2026년 7월 6일 10:00
감지된 SKILL.md 언어: 영어
스타: 1
포크: 0

설치 방법

기본적으로 소스를 먼저 확인하는 Prompt가 선택됩니다. 직접 명령으로 전환하거나 로컬 사본을 다운로드할 수도 있습니다.

Codex 또는 Claude로 설치 이 Prompt를 복사해 Codex, Claude 또는 다른 어시스턴트에 붙여 넣으면 Skill 페이지를 검토하고 설치를 진행할 수 있습니다.

직접 명령은 검토 Prompt를 거치지 않습니다. 실행하기 전에 소스를 확인하세요.

npx skills add https://github.com/caixinyu2017-star/xixi --skill full-empirical-analysis-skill-r

명령은 한 줄로 유지됩니다. 복사하기 전에 가로로 스크롤해 전체 내용을 확인하세요.

로컬 사본을 원하시나요? SkillsMP에서 현재 제공할 수 있는 파일을 다운로드하세요.

소스 파일 검토

설치 여부를 결정하기 전에 SKILL.md와 SkillsMP에 표시된 보조 파일을 읽어 보세요.

파일 탐색기

9 개 파일

SKILL.md

소스 지침 · 읽기 전용 미리보기

검색으로 돌아가기 설치로 돌아가기

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)", main = "Event study: dynamic effect of training", ref.line = -0.5) dev.off() # ============================================================ # 8h. FIGURE 4 — Sensitivity / robustness curve # (HonestDiD / spec curve / forest of robustness battery) # ============================================================ # HonestDiD example (after the event study with stored b/V): library(HonestDiD) honest_out <- createSensitivityResults(betahat = es$coefficients, sigma = vcov(es), numPrePeriods = 5, numPostPeriods = 5, Mbarvec = seq(0, 0.5, by = 0.05)) sens_plot <- createSensitivityPlot(honest_out, originalResults = honest_out$mainResult) ggsave("figures/fig4_sensitivity.pdf", plot = sens_plot, width = 7, height = 4) ggsave("figures/fig4_sensitivity.png", plot = sens_plot, width = 7, height = 4, dpi = 300) # Alternative — robustness forest plot: # rob_summary <- imap_dfr(rob, ~ tibble( # group = .y, # est = coef(.x)[1], # se = se(.x)[1] # )) # ggplot(rob_summary, aes(est, fct_rev(factor(group)))) + # geom_point(size = 3, color = "navy") + # geom_errorbarh(aes(xmin = est - 1.96*se, xmax = est + 1.96*se), # height = 0.2, color = "navy") + # geom_vline(xintercept = 0, linetype = "dashed") + # labs(x = "Coefficient on training (95% CI)", y = NULL, # title = "Robustness forest plot") # ggsave("figures/fig4_sensitivity.pdf", width = 7, height = 4) # ============================================================ # 8i. FIGURE 1 — Trend / motivation (treated vs control over time) # ============================================================ 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", color = "gray40") + scale_color_manual(values = c("0" = "darkred", "1" = "navy"), labels = c("Control", "Treated"), name = "") + labs(x = "Year", y = "Mean log wage", title = "Treated vs control trend") + theme_classic(base_size = 11) + theme(legend.position = "bottom")

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

소스 정보

설치 방법

소스 파일 검토

이 저장소의 다른 Skills

이 저장소의 다른 Skills

Full Empirical Analysis — Classical R Workflow

Philosophy

Full Empirical Analysis — Classical R Workflow

Philosophy

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)