基于R语言MICE多重插补后的逻辑与有序模型AUC计算问题

一、逻辑回归模型（二分类结局） #

1. 计算单插补数据集的指标 #

借助pROC包计算AUC，rcompanion或pscl包计算伪R²（如Nagelkerke R²），在mice的with()函数中一次性完成模型拟合与指标计算：

library(mice)
library(pROC)
library(rcompanion)
library(pscl)

# 拟合模型并提取每个插补集的AUC、伪R²
imputed_model_results <- with(imputed_Data, {
  model <- glm(Poor ~ age + sex + education + illness + injurycause, family = "binomial")
  pred_prob <- predict(model, type = "response")
  # 计算AUC
  auc_val <- roc(Poor ~ pred_prob, quiet = TRUE)$auc
  # 计算Nagelkerke伪R²（也可替换为pscl包的pR2(model)$McFadden）
  rsq_val <- nagelkerke(model)$Pseudo.R.squared.for.model.vs.null[3]
  list(model = model, auc = auc_val, rsq = rsq_val)
})

2. 合并插补后的指标 #

针对AUC和R²，可采用两种方式合并：

• 简单均值法：插补数量足够时，均值能提供稳定估计

# 提取所有插补集的指标列表
auc_list <- sapply(imputed_model_results$analyses, function(x) x$auc)
rsq_list <- sapply(imputed_model_results$analyses, function(x) x$rsq)

# 计算合并值与标准误
pooled_auc <- mean(auc_list)
auc_se <- sd(auc_list)/sqrt(length(auc_list))

pooled_rsq <- mean(rsq_list)
rsq_se <- sd(rsq_list)/sqrt(length(rsq_list))

# 输出结果
cat("合并后的AUC：", round(pooled_auc, 3), "（SE：", round(auc_se, 3), "）\n")
cat("合并后的Nagelkerke R²：", round(pooled_rsq, 3), "（SE：", round(rsq_se, 3), "）\n")

• Rubin规则合并：严格遵循多重插补的合并逻辑，计算总方差与置信区间

# Rubin规则合并AUC
within_var <- var(auc_list)
between_var <- sum((auc_list - pooled_auc)^2)/(length(auc_list)-1)
total_var <- within_var/length(auc_list) + between_var*(1 + 1/length(auc_list))
auc_ci <- pooled_auc + c(-1,1)*qnorm(0.975)*sqrt(total_var)

cat("Rubin规则合并AUC：", round(pooled_auc,3), "（95%CI：", paste(round(auc_ci,3), collapse="-"), "）\n")

二、有序回归模型（多分类有序结局） #

对于有序逻辑回归（如polr模型），用multiclass.roc计算多分类AUC，伪R²计算逻辑与二分类一致：

library(MASS)

imputed_ordinal_results <- with(imputed_Data, {
  model <- polr(Ordered_Poor ~ age + sex + education + illness + injurycause, Hess = TRUE)
  pred_prob <- predict(model, type = "probs")
  # 采用One-vs-All方式计算多分类AUC（按需调整对应类别概率）
  auc_val <- multiclass.roc(Ordered_Poor, pred_prob[,2])$auc
  rsq_val <- nagelkerke(model)$Pseudo.R.squared.for.model.vs.null[3]
  list(model = model, auc = auc_val, rsq = rsq_val)
})

# 合并指标
auc_ordinal_list <- sapply(imputed_ordinal_results$analyses, function(x) x$auc)
rsq_ordinal_list <- sapply(imputed_ordinal_results$analyses, function(x) x$rsq)

pooled_auc_ordinal <- mean(auc_ordinal_list)
pooled_rsq_ordinal <- mean(rsq_ordinal_list)

三、增量预测效益计算 #

对比基础模型与纳入新变量后的拓展模型，通过AUC差值、R²差值评估增量效益：

# 拟合基础模型（仅含核心变量）
imputed_base_model <- with(imputed_Data, glm(Poor ~ age + sex, family = "binomial"))
# 拟合拓展模型（加入目标变量）
imputed_extended_model <- with(imputed_Data, glm(Poor ~ age + sex + education + illness + injurycause, family = "binomial"))

# 提取两个模型的AUC列表
base_auc_list <- sapply(imputed_base_model$analyses, function(x) roc(x$y, predict(x, type="response"), quiet=TRUE)$auc)
extended_auc_list <- sapply(imputed_extended_model$analyses, function(x) roc(x$y, predict(x, type="response"), quiet=TRUE)$auc)

# 计算每个插补集的AUC差值并合并
auc_diff_list <- extended_auc_list - base_auc_list
pooled_auc_diff <- mean(auc_diff_list)
diff_se <- sd(auc_diff_list)/sqrt(length(auc_diff_list))

cat("拓展模型相对基础模型的AUC增量：", round(pooled_auc_diff,3), "（SE：", round(diff_se,3), "）\n")

注意事项 #

• 确保结局变量无缺失（mice默认处理，可提前用complete.cases()检查）
• 伪R²选择：Nagelkerke R²适合二分类模型解释拟合优度，McFadden R²更侧重模型相对 null 模型的提升
• 多分类有序模型的AUC计算方式可按需调整（One-vs-One或One-vs-All）

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