使用R中MICE多重插补后，零膨胀泊松模型合并估计异常问题

默认的mice::pool()和summary()方法无法适配pscl::zeroinfl模型的双结构（计数部分+零膨胀部分），导致无法清晰展示完整结果。以下是两种可行的解决方案：

方法1：手动拆分模型部分并分别汇总 #

直接提取每个插补模型的计数（count）和零膨胀（zero）部分系数，分别进行合并和汇总：

# 提取每个插补模型的计数部分系数
count_coefs <- lapply(res.zinb$analyses, function(x) coef(x, model = "count"))
# 提取零膨胀部分系数
zero_coefs <- lapply(res.zinb$analyses, function(x) coef(x, model = "zero"))

# 为计数部分构造可被pool处理的mira对象
pool_count <- pool(as.mira(lapply(count_coefs, function(x) {
  temp_mod <- lm(x ~ 0)
  temp_mod$coefficients <- x
  temp_mod$model <- data.frame(matrix(NA, nrow = nrow(dat), ncol = length(x)))
  colnames(temp_mod$model) <- names(x)
  temp_mod
})))

# 为零膨胀部分构造可被pool处理的mira对象
pool_zero <- pool(as.mira(lapply(zero_coefs, function(x) {
  temp_mod <- lm(x ~ 0)
  temp_mod$coefficients <- x
  temp_mod$model <- data.frame(matrix(NA, nrow = nrow(dat), ncol = length(x)))
  colnames(temp_mod$model) <- names(x)
  temp_mod
})))

# 分部分输出结果
cat("=== 计数模型结果 ===\n")
print(summary(pool_count))

cat("\n=== 零膨胀模型结果 ===\n")
print(summary(pool_zero))

方法2：用broom包简化结果提取与合并 #

借助broom包的tidy功能，统一提取模型结果后按Rubin规则合并：

library(broom)
library(tidyr)

# 提取每个插补模型的结构化结果
tidy_list <- lapply(res.zinb$analyses, function(mod) {
  tidy(mod, conf.int = TRUE) %>%
    # 区分计数/零膨胀部分
    mutate(model_part = ifelse(grepl("^zero", term), "zero", "count")) %>%
    # 拆分变量名中的模型标识
    separate(term, into = c("ignore", "variable"), sep = "\\.", extra = "merge", fill = "left") %>%
    mutate(variable = ifelse(is.na(variable), ignore, variable)) %>%
    select(model_part, variable, estimate, std.error, conf.low, conf.high)
})

# 合并所有插补结果并按Rubin规则计算汇总值
pooled_results <- bind_rows(tidy_list) %>%
  group_by(model_part, variable) %>%
  summarise(
    estimate = mean(estimate),
    # Rubin规则计算合并标准误
    std.error = sqrt(mean(std.error^2) + (1 + 1/n()) * var(estimate)),
    conf.low = mean(conf.low),
    conf.high = mean(conf.high),
    p.value = 2 * pt(-abs(estimate / std.error), df = n() - 1)
  )

print(pooled_results)

原理说明 #

zeroinfl模型的系数包含两个独立部分：用于预测计数的count模型和用于预测零值概率的zero模型。默认的pool()函数会将这两部分系数混为一谈，无法区分展示；上述方法通过拆分两个部分并分别处理，就能得到清晰的完整结果。

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