在R中自动化面板回归：变量设置与NA值处理问题

原代码存在变量提取逻辑错误、字符串变量无法直接嵌入公式、基线控制变量未与原数据集匹配等问题，以下是针对性的修正方案，同时解决NA值处理需求：

1. 预处理：提取终末期因变量作为基线控制变量 #

先从面板数据中提取round == 11时的因变量值，匹配到每个个体的所有观测行，同时计算该变量的全局均值：

library(lfe)
library(dplyr)

outcome <- "ps_depression_cesd"

# 提取round=11的因变量值，按个体标识（hhid）匹配到全量数据
df <- df %>%
  group_by(hhid) %>%
  mutate(baseline_outcome = get(outcome)[round == 11]) %>%
  ungroup()

# 计算终末期因变量的均值（排除NA）
mean_baseline <- df %>%
  filter(round == 11) %>%
  pull(outcome) %>%
  mean(na.rm = TRUE)

2. 动态构造回归公式 #

因outcome是字符串变量，需用reformulate动态生成公式，避免直接传入felm时触发报错：

# 定义固定效应与聚类项
fixed_effs <- "level1 + round_cont"
cluster_terms <- "level4 + hhid"

# 回归1：仅核心处理效应
formula1 <- reformulate("treatment", response = outcome)
reg_1 <- felm(formula1 | fixed_effs | 0 | cluster_terms, data = df, na.action = na.exclude)

# 回归2：处理效应+基线控制+均值项
formula2 <- reformulate(
  c("treatment", "baseline_outcome", paste0("I(", mean_baseline, ")")),
  response = outcome
)
reg_2 <- felm(formula2 | fixed_effs | 0 | cluster_terms, data = df, na.action = na.exclude)

# 回归3：分类型处理效应
formula3 <- reformulate(c("treat_uct", "treat_ffa"), response = outcome)
reg_3 <- felm(formula3 | fixed_effs | 0 | cluster_terms, data = df, na.action = na.exclude)

# 回归4：分类型处理效应+基线控制+均值项
formula4 <- reformulate(
  c("treat_uct", "treat_ffa", "baseline_outcome", paste0("I(", mean_baseline, ")")),
  response = outcome
)
reg_4 <- felm(formula4 | fixed_effs | 0 | cluster_terms, data = df, na.action = na.exclude)

3. NA值处理说明 #

• 回归时设置na.action = na.exclude，保留样本结构的同时自动剔除含NA的观测；
• 计算均值时用na.rm = TRUE忽略NA值，避免均值计算失效；
• 若部分个体无round == 11的观测，baseline_outcome会生成NA，可根据研究需求用tidyr::fill填充或在回归中自动过滤。

批量自动化封装（可选） #

如果需要处理多个因变量，可封装为函数实现批量回归：

run_panel_regs <- function(df, outcome_var) {
  # 预处理基线控制变量
  df_processed <- df %>%
    group_by(hhid) %>%
    mutate(baseline_outcome = get(outcome_var)[round == 11]) %>%
    ungroup()
  
  mean_baseline <- df_processed %>%
    filter(round == 11) %>%
    pull(outcome_var) %>%
    mean(na.rm = TRUE)
  
  # 定义公式集合
  formulas <- list(
    reformulate("treatment", response = outcome_var),
    reformulate(c("treatment", "baseline_outcome", paste0("I(", mean_baseline, ")")), response = outcome_var),
    reformulate(c("treat_uct", "treat_ffa"), response = outcome_var),
    reformulate(c("treat_uct", "treat_ffa", "baseline_outcome", paste0("I(", mean_baseline, ")")), response = outcome_var)
  )
  
  # 批量运行回归
  reg_list <- lapply(formulas, function(f) {
    felm(f | level1 + round_cont | 0 | level4 + hhid, data = df_processed, na.action = na.exclude)
  })
  
  names(reg_list) <- paste0("reg_", 1:4)
  return(reg_list)
}

# 使用示例
all_regressions <- run_panel_regs(df, "ps_depression_cesd")

内容的提问来源于stack exchange，提问作者Paula