如何从千次随机化回归结果中提取指定变量的系数、标准误与P值

问题背景 #

已通过data.table生成模拟数据，完成以下步骤：

1. 按城市分组，将village_status随机化1000次，生成randomvil1至randomvil1000列；
2. 使用estimatr包的lm_robust执行1000次回归，回归公式为y~ randomvil + randomvil:x1 + x2 + x3 + factor(city)；

需要提取每次回归中randomvil和randomvil:x1的系数（Estimate）、标准误（Std. Error）及P值（Pr(>|t|)），整理为指定格式的data.table。

解决方案 #

使用以下代码提取并整理结果：

# 定义提取函数：从单个回归结果中获取目标变量的统计量
extract_coefs <- function(fit_obj, var_name) {
  # 把系数表转换为data.table，保留变量名列
  coef_table <- as.data.table(fit_obj$coefficients, keep.rownames = "term")
  # 筛选出目标变量行
  target_coefs <- coef_table[term %in% c("randomvil", "randomvil:x1")]
  # 将变量名替换为对应randomvilN格式
  target_coefs[, term := gsub("randomvil", var_name, term)]
  # 保留需要的列并重命名为指定格式
  target_coefs[, .(
    term,
    Estimate = estimate,
    `Std. Error` = std.error,
    `Pr(>|t|)` = p.value
  )]
}

# 批量处理所有回归结果
result_list <- mapply(
  FUN = extract_coefs,
  fit_obj = fit,
  var_name = names(fit),
  SIMPLIFY = FALSE
)

# 合并所有结果为一个data.table
final_result <- rbindlist(result_list)

代码说明 #

• extract_coefs函数：针对单个回归结果，筛选出randomvil和randomvil:x1的行，将变量名修正为randomvilN和randomvilN:x1格式，提取并重命名所需的统计量列；
• mapply：同步遍历回归结果列表fit和对应的随机变量名（如randomvil1），批量执行提取操作；
• rbindlist：将所有提取出的小表合并为一个完整的data.table，结构与需求示例一致。

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