使用mice包合并多重插补下多分类逻辑回归结果报错求助

可复现问题代码 #

# 加载依赖包
library(mice)
library(nnet)
library(broom)

# 构造带缺失值的多分类数据集
data(iris)
iris$Sepal.Length[sample(1:nrow(iris), 20)] <- NA
iris$Species <- factor(iris$Species, levels = c("setosa", "versicolor", "virginica"))

# 执行多重插补
imp <- mice(iris, m = 5, seed = 123)

# 拟合多分类逻辑回归模型
fit <- with(imp, multinom(Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width))

# 尝试合并结果触发错误
tidy(fit) # 报错：No tidy method recognized for this list

问题原因 #

broom包的tidy函数未适配mice返回的多分类逻辑回归模型列表；pool.scalar需手动处理每个系数的估计值与方差，直接调用无法识别多分类模型的结构。

解决方法 #

方法1：手动提取系数后用pool.scalar合并

# 定义提取单模型系数与方差的函数
extract_coefs <- function(model) {
  coef_mat <- coef(model)
  var_vec <- diag(vcov(model))
  
  # 整理为标准化数据框
  data.frame(
    category = rep(rownames(coef_mat), each = ncol(coef_mat)),
    variable = rep(colnames(coef_mat), nrow(coef_mat)),
    estimate = as.vector(coef_mat),
    std.error = sqrt(var_vec)
  )
}

# 提取所有插补模型的系数
coef_list <- lapply(fit$analyses, extract_coefs)
coef_df <- do.call(rbind, coef_list)
coef_df$.imp <- rep(1:5, each = nrow(coef_list[[1]]))

# 对每个系数分组执行pool.scalar合并
pooled_results <- lapply(split(coef_df, paste(coef_df$category, coef_df$variable)), function(x) {
  pool.scalar(Q = x$estimate, U = x$std.error^2, n = nrow(iris))
})

# 整理最终结果
final_df <- do.call(rbind, lapply(names(pooled_results), function(nm) {
  parts <- strsplit(nm, " ")[[1]]
  res <- pooled_results[[nm]]
  data.frame(
    category = parts[1],
    variable = paste(parts[-1], collapse = " "),
    estimate = res$qbar,
    std.error = sqrt(res$t),
    p.value = 2 * pnorm(-abs(res$qbar / sqrt(res$t)))
  )
}))

print(final_df)

方法2：使用miceadds包简化合并流程

miceadds的pool.mice函数专门适配多分类模型的插补结果合并：

library(miceadds)

# 直接合并插补后的模型结果
pooled_fit <- pool.mice(fit)

# 输出合并后的详细结果
summary(pooled_fit)

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