R中调查加权多重插补数据如何计算合并平均边际效应

当前处理带抽样权重的调查数据，已通过mice()完成缺失值多重插补，后续需要拟合含复杂交互项的模型，计算模型的平均边际效应（AME）。
该操作在Stata中流程简便，但需要在R环境完成分析。已知可以对每个单独插补数据集分别计算AME后对估计值取平均，但要得到正确标准误，必须调用mice包的pool()函数完成结果合并。

可复现问题的示例代码 #

library(tidyverse)
library(survey)
library(mice)
library(margins)

df <- tibble(y = c(0, 5, 0, 4, 0, 1, 2, 3, 1, 12), region = c(1, 1, 1, 1, 1, 3, 3, 3, 3, 3), 
             weight = c(7213, 2142, 1331, 4342, 9843, 1231, 1235, 2131, 7548, 2348), 
             x1 = c(1.14, 2.42, -0.34, 0.12, -0.9, -1.2, 0.67, 1.24, 0.25, -0.3),
             x2 = c(12, NA, 10, NA, NA, 12, 11, 8, 9, 9))

已验证可正常运行的场景 #

• 对单一非多重插补数据拟合的svyglm对象调用margins()可正常运行
• 通过with()基于每个插补数据集拟合svyglm后，用pool()合并回归结果可正常执行，对应代码：

m <- with(surv_obj, svyglm(y ~ x1 * x2))
pool(m)

报错/警告场景 #

1. 在with()中嵌套调用margins()，返回错误：Error in .svycheck(design) : argument "design" is missing, with no default，对应代码：

with(surv_obj, margins(svyglm(y ~ x1 * x2), design = surv_obj))

2. 在svyglm调用中显式指定design参数，返回错误：Error in UseMethod("svyglm", design) : no applicable method for 'svyglm' applied to an object of class "svyimputationList"，对应代码：

with(surv_obj, margins(svyglm(y ~ x1 * x2, design = surv_obj), design = surv_obj))

3. 移除调查设计层，直接对每个插补数据集拟合lm模型、调用margins()计算边际效应后用pool()合并，出现警告：Warning in get.dfcom(object, dfcom) : Infinite sample size assumed.，对应代码：

m1 <- with(imputed_df, margins(lm(y ~ x1 * x2)))
pool(m1)

pool()计算过程需要调用样本量参数，该警告会导致结果可靠性存疑。

方案a：组合函数实现AME自动合并 #

with()处理多重插补对象时，会在每个插补子集的独立环境中运行表达式，全局环境中的调查设计对象无法被正确识别，因此嵌套调用会报错。正确流程为逐个处理每个插补数据集，将结果整理为mice支持的mira类对象后传入pool()即可。
完整可运行代码：

# 完成多重插补，可自行调整插补次数
imp <- mice(df, m = 5, printFlag = F)
# 构建多重插补调查设计对象
surv_obj <- svydesign(ids = ~1, weights = ~weight, data = imp)

# 逐个处理每个插补数据集
ame_list <- lapply(1:imp$m, function(i) {
  imp_data <- complete(imp, i)
  design_i <- svydesign(ids = ~1, weights = ~weight, data = imp_data)
  fit_i <- svyglm(y ~ x1 * x2, design = design_i)
  ame_i <- margins(fit_i, design = design_i)
  list(model = fit_i, ame = summary(ame_i))
})

# 整理为pool可识别的mira类结构
ame_mira <- list(
  call = call("margins"),
  nmis = imp$nmis,
  analyses = lapply(ame_list, function(x) {
    res <- list(
      coefficients = setNames(x$ame$AME, x$ame$factor),
      vcov = vcov(x$ame),
      df.residual = df.residual(x$model),
      call = x$model$call,
      terms = x$model$terms
    )
    class(res) <- c("svyglm", "glm", "lm")
    res
  }),
  m = imp$m
)
class(ame_mira) <- "mira"

# 合并结果，无无限样本量警告
pooled_ame <- pool(ame_mira)
summary(pooled_ame)

方案b：手动提取要素合并结果 #

如果不构造适配pool()的对象结构，可直接提取margins()输出的核心统计量，传入pool.scalar()按Rubin规则合并。需要提取的要素仅两类：

• 每个边际效应在各插补数据集的点估计值
• 每个边际效应在各插补数据集的抽样方差（即标准误的平方）

手动合并代码如下：

# 提取各插补集的AME估计值和方差
ame_est <- lapply(ame_list, function(x) setNames(x$ame$AME, x$ame$factor))
ame_var <- lapply(ame_list, function(x) {
  ses <- x$ame$SE
  setNames(ses^2, x$ame$factor)
})

# 逐变量调用pool.scalar合并
var_names <- names(ame_est[[1]])
pooled_res <- lapply(var_names, function(v) {
  q <- sapply(ame_est, function(x) x[v])
  u <- sapply(ame_var, function(x) x[v])
  ps <- pool.scalar(Q = q, U = u, n = nrow(df))
  data.frame(
    term = v,
    estimate = ps$qbar,
    std.error = sqrt(ps$t),
    statistic = ps$qbar / sqrt(ps$t),
    p.value = 2 * pnorm(abs(ps$qbar / sqrt(ps$t)), lower.tail = F),
    df = ps$df
  )
})
pooled_res <- do.call(rbind, pooled_res)

两种方法得到的结果完全一致，传入样本量参数后不会再触发无限样本量警告，标准误计算完全符合多重插补Rubin合并规则，同时正确纳入了调查权重的影响。

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