关于多重插补后使用pool.scalar()进行多中介多元分析的方法咨询

嗨，我完全懂你现在的困扰——用mice做完多重插补，想靠bruceR的PROCESS函数跑5个平行中介的分析，结果因为bruceR不兼容broom，用mice::pool直接报错，想自己用pool.scalar()手动合并结果，却找不到多变量场景的指南。别慌，我一步步给你拆解操作流程，亲测可行：

第一步：在每个插补数据集上批量运行PROCESS分析，提取关键统计量 #

首先你需要遍历所有插补后的数据集，在每个数据集上跑PROCESS，并提取每个效应（间接、直接、总效应）的估计值和方差（方差=标准误的平方，pool.scalar()需要这两个值）。

library(bruceR)
library(mice)

# 假设你的多重插补结果对象是imputed_data
# 定义单个插补数据集的分析函数
run_single_process <- function(data) {
  # 替换成你的多中介模型语法，model=4对应平行多重中介
  med_model <- PROCESS(data = data,
                       y = "outcome", x = "predictor",
                       m = c("med1", "med2", "med3", "med4", "med5"),
                       model = 4)
  
  # 提取间接效应的估计值和标准误，计算方差
  indirect_res <- med_model$indirect_effects
  indirect_ests <- indirect_res$Estimate
  indirect_vars <- (indirect_res$SE)^2
  
  # 提取直接、总效应的估计值和方差
  direct_est <- med_model$direct_effect$Estimate
  direct_var <- (med_model$direct_effect$SE)^2
  total_est <- med_model$total_effect$Estimate
  total_var <- (med_model$total_effect$SE)^2
  
  # 返回所有需要的统计量
  list(
    indirect_ests = indirect_ests,
    indirect_vars = indirect_vars,
    direct_est = direct_est,
    direct_var = direct_var,
    total_est = total_est,
    total_var = total_var,
    med_names = indirect_res$Effect  # 保存中介效应名称，方便后续对应
  )
}

# 在所有插补数据集上批量运行分析
all_imp_results <- lapply(complete(imputed_data, "all"), run_single_process)

第二步：整理所有插补数据集的统计量 #

把每个效应的估计值和方差分别整理成可批量处理的格式：

m <- length(all_imp_results)  # 插补数据集的数量
med_names <- all_imp_results[[1]]$med_names  # 中介效应名称

# 整理间接效应的估计值和方差矩阵（每行对应一个插补数据集）
indirect_ests_mat <- do.call(rbind, lapply(all_imp_results, function(x) x$indirect_ests))
indirect_vars_mat <- do.call(rbind, lapply(all_imp_results, function(x) x$indirect_vars))

# 整理直接、总效应的估计值和方差向量
direct_ests <- sapply(all_imp_results, function(x) x$direct_est)
direct_vars <- sapply(all_imp_results, function(x) x$direct_var)
total_ests <- sapply(all_imp_results, function(x) x$total_est)
total_vars <- sapply(all_imp_results, function(x) x$total_var)

第三步：用pool.scalar()逐个合并每个效应的结果 #

pool.scalar()本身是处理单变量的，但我们可以循环处理每个中介的间接效应，以及直接、总效应：

# 定义单个效应的合并函数
pool_single_eff <- function(ests, vars, m, original_n) {
  # 调用pool.scalar合并结果
  pooled_res <- pool.scalar(Q = ests, U = vars, n = original_n, m = m)
  # 计算标准误和95%置信区间（正态近似）
  se <- sqrt(pooled_res$t)
  ci_low <- pooled_res$qbar - 1.96 * se
  ci_high <- pooled_res$qbar + 1.96 * se
  # 计算p值（正态近似）
  p_val <- 2 * pnorm(-abs(pooled_res$qbar / se))
  
  # 返回结构化结果
  data.frame(
    Estimate = pooled_res$qbar,
    SE = se,
    CI_Low = ci_low,
    CI_High = ci_high,
    p_value = p_val
  )
}

# 合并所有间接效应
pooled_indirect <- lapply(1:length(med_names), function(i) {
  pool_single_eff(indirect_ests_mat[,i], indirect_vars_mat[,i], m, nrow(imputed_data$data))
})
names(pooled_indirect) <- med_names

# 合并直接、总效应
pooled_direct <- pool_single_eff(direct_ests, direct_vars, m, nrow(imputed_data$data))
pooled_total <- pool_single_eff(total_ests, total_vars, m, nrow(imputed_data$data))

第四步：整理合并后的结果并解读 #

把所有结果整合成清晰的表格，方便查看：

# 整理间接效应结果
indirect_df <- do.call(rbind, pooled_indirect)
indirect_df$Effect_Type <- "Indirect"
indirect_df$Mediator <- rownames(indirect_df)

# 整理直接、总效应结果
direct_df <- pooled_direct
direct_df$Effect_Type <- "Direct"
direct_df$Mediator <- "None"

total_df <- pooled_total
total_df$Effect_Type <- "Total"
total_df$Mediator <- "None"

# 合并所有结果
final_results <- rbind(indirect_df, direct_df, total_df)
# 调整列顺序，让结果更易读
final_results <- final_results[, c("Effect_Type", "Mediator", "Estimate", "SE", "CI_Low", "CI_High", "p_value")]

# 查看最终结果
print(final_results)

一些注意事项 #

• 先单独跑一次PROCESS，用str(med_model)查看输出结构，确保你提取的Estimate和SE路径正确——不同版本的bruceR可能输出位置略有变化。
• 这里用的是正态近似计算置信区间和p值，如果想要更稳健的结果，可以考虑在每个插补数据集上做Bootstrap再合并样本，但操作会更繁琐，上面的方法是多重插补后中介分析的标准做法。
• pool.scalar()中的n参数必须用原始数据的样本量（也就是nrow(imputed_data$data)），这样自由度调整才正确。

备注：内容来源于stack exchange，提问作者NessD