基于多重插补数据的R语言似不相关回归(SUR)结果合并问题

嘿，我来帮你解决用mice生成的多重插补数据结合systemfit做似不相关回归（SUR）的问题！下面是完整的实操步骤和代码示例，一步步来就行：

首先把需要的包和示例数据搞定，用你提到的nhanes2数据集来演示：

library(systemfit)
library(mice)
library(broom) # 辅助提取模型结果，可选但方便
library(miceadds) # 专门支持systemfit的多重插补合并，推荐安装

# 加载内置的缺失数据集
data(nhanes2)

用mice处理缺失值，这里生成5个插补后的完整数据集（你可以根据需求调整m的数量）：

# 生成插补数据，关闭打印日志让输出更干净
imp <- mice(nhanes2, m = 5, printFlag = FALSE)

按照你的需求定义两个回归方程，然后打包成systemfit需要的列表格式：

# 定义两个回归方程
r1 <- bmi ~ hyp
r2 <- bmi ~ age
# 构建SUR模型系统
system_models <- list(r1 = r1, r2 = r2)

循环遍历每个插补后的数据集，分别拟合SUR模型，把结果存起来：

# 创建空列表存储每个插补模型的结果
sur_fits <- list()

# 逐个处理插补数据集
for (i in 1:imp$m) {
  # 提取第i个插补后的完整数据
  complete_dat <- complete(imp, i)
  # 拟合SUR模型，method="SUR"指定用似不相关回归
  fit <- systemfit(system_models, data = complete_dat, method = "SUR")
  # 把结果存入列表
  sur_fits[[i]] <- fit
}

这一步是关键，必须按照Rubin规则合并多个插补模型的结果，这里给你两种方法：

方法一：用miceadds一键合并（推荐） #

miceadds包专门做了pool.sysfit函数，直接帮你搞定合并，省心又准确：

# 合并所有插补的SUR结果
pooled_result <- pool.sysfit(sur_fits)

# 查看最终合并后的系数、标准误等信息
summary(pooled_result)

方法二：手动合并（无需额外包） #

如果不想装新包，也可以手动按照Rubin规则计算合并结果：

# 提取每个模型的系数
coef_list <- lapply(sur_fits, coef)
coef_df <- do.call(rbind, coef_list)

# 提取每个模型的方差协方差矩阵
vcov_list <- lapply(sur_fits, vcov)

# 按照Rubin规则计算合并指标
m <- imp$m # 插补数量
pooled_coef <- colMeans(coef_df) # 合并系数
W <- apply(simplify2array(vcov_list), 1:2, mean) # 组内方差
B <- var(coef_df) # 组间方差
total_var <- W + (1 + 1/m)*B # 总方差
pooled_se <- sqrt(diag(total_var)) # 合并标准误

# 整理成直观的结果数据框
final_results <- data.frame(
  系数名称 = names(pooled_coef),
  估计值 = pooled_coef,
  标准误 = pooled_se,
  stringsAsFactors = FALSE
)

# 查看结果
print(final_results)

• 绝对不能直接用带缺失值的原始数据跑SUR！必须对每个插补后的完整数据集单独拟合，再合并结果，这是多重插补的核心要求。
• 如果你的SUR模型有特殊设定（比如工具变量、加权等），只需要调整systemfit的参数就行，循环和合并的逻辑完全不变。

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