PAYATHON 2026

R 中使用 mice 多重插补后,如何计算并合并 glmer() 的 AUC?

技术老齐

可以先在每个插补数据集上用 predict() 提取 glmer() 的预测概率,再用 pROC::roc() 计算 AUC。如果还需要合并置信区间,就要估计每个 AUC 的插补内方差,再按 Rubin 规则合并。

pROC::roc() 接收结局和预测分数,因此也适用于 glmer() 的预测结果。bobyqa 影响模型拟合过程,不改变 AUC 的计算方式。计算前,需要确定预测是否包含随机效应,以及如何处理组内相关性。

先确定要评价哪一种预测

随机截距模型有两种预测方式:

# 包含估计的随机效应
predict(fit, newdata = dat, type = "response", re.form = NULL)

# 将随机效应设为 0,仅使用固定效应
predict(fit, newdata = dat, type = "response", re.form = NA)

前者利用已有组的随机效应信息;后者不依赖已估计的组效应。两种方式算出的 AUC 可能不同,所有插补数据集必须采用同一种预测方式。

re.form = NA 将随机效应设为零,并未对随机效应分布积分,因此得到的不是边际概率。

下面的示例包含随机效应,并在拟合数据上计算 AUC。这样得到的是表观 AUC,通常比模型在新数据上的表现乐观。

在各个插补数据集上计算 AUC

假设 outcome 已编码为数值型 0/1,其中 1 表示事件;模型变量已完成插补,randomeffect 是正确的分组标识。

library(mice)
library(lme4)

model <- with(
  mids_data,
  glmer(
    outcome ~ var1 + var2 + (1 | randomeffect),
    family = binomial,
    control = glmerControl(optimizer = "bobyqa"),
    na.action = na.fail
  )
)

datasets <- lapply(
  seq_len(mids_data$m),
  function(i) mice::complete(mids_data, action = i)
)

fits <- model$analyses

stopifnot(length(fits) == length(datasets), length(fits) >= 2L)

auc_value <- function(fit, dat) {
  y <- dat$outcome

  stopifnot(
    is.numeric(y),
    !anyNA(y),
    all(y %in% c(0, 1)),
    length(unique(y)) == 2L
  )

  prob <- predict(
    fit,
    newdata = dat,
    type = "response",
    re.form = NULL
  )

  stopifnot(
    length(prob) == nrow(dat),
    all(is.finite(prob))
  )

  roc_obj <- pROC::roc(
    response = y,
    predictor = prob,
    levels = c(0, 1),
    direction = "<",
    quiet = TRUE
  )

  as.numeric(pROC::auc(roc_obj))
}

A <- vapply(
  seq_along(fits),
  function(i) auc_value(fits[[i]], datasets[[i]]),
  numeric(1)
)

data.frame(
  imputation = seq_along(A),
  AUC = A
)

指定 levels 和 direction 后,较大的预测概率对应事件 1,避免每次插补自动选择 ROC 方向。

如果结局是因子,应先核对事件水平,再显式转换,不要直接使用 as.numeric(factor)。

用按组 bootstrap 估计插补内方差

计算合并标准误时,除了各次插补的 AUC,Rubin 规则还需要两部分不确定性:

  • 每个插补数据集内部的估计方差;
  • 不同插补数据集之间的估计差异。

普通 DeLong 方差依赖观测独立等条件。混合效应模型中的同组观测通常相关,直接使用 pROC::var(..., method = "delong") 可能低估不确定性。

对于这里的表观 AUC,可以考虑用按组 bootstrap 近似估计方差:在每个完成插补的数据集中,按独立组重采样,每次重新拟合模型并计算 AUC。以下代码假设组之间独立,且 randomeffect 是唯一的聚类层级。

resample_clusters <- function(dat) {
  rows <- split(
    seq_len(nrow(dat)),
    dat$randomeffect,
    drop = TRUE
  )

  G <- length(rows)
  if (G < 2L) stop("至少需要两个组。")

  selected <- sample.int(G, size = G, replace = TRUE)

  pieces <- lapply(seq_along(selected), function(j) {
    block <- dat[rows[[selected[j]]], , drop = FALSE]

    # 同一个原始组可能被抽中多次。
    # 每份副本需要独立的 bootstrap 组标识。
    block$randomeffect <- paste0("boot_", j)
    block
  })

  out <- do.call(rbind, pieces)
  out$randomeffect <- factor(out$randomeffect)
  out
}

logit_auc <- function(a) {
  if (any(!is.finite(a) | a <= 0 | a >= 1)) {
    stop("AUC 必须严格位于 0 和 1 之间,才能使用 logit 变换。")
  }
  qlogis(a)
}

Q <- logit_auc(A)

set.seed(2026)
B_boot <- 500L  # 示例次数,实际分析应检查结果的稳定性

U <- vapply(seq_along(fits), function(i) {
  z_boot <- replicate(B_boot, {
    dat_b <- resample_clusters(datasets[[i]])
    fit_b <- update(fits[[i]], data = dat_b)

    logit_auc(auc_value(fit_b, dat_b))
  })

  stats::var(z_boot)
}, numeric(1))

方差直接在 logit 尺度上估计,因此 Q 和 U 的尺度一致。每次重采样后都重新拟合模型,也就纳入了重新估计模型带来的变化。

代码遇到错误会停止。如果某次重采样只有一种结局、模型拟合失败,或 AUC 恰好等于 0 或 1,应检查数据和方法。不能直接删除这些重采样结果,也不应随意截断边界值后继续计算。收敛警告和奇异拟合同样需要检查。

按 Rubin 规则合并,再变换回 AUC

对第 (i) 个插补数据集,记:

[ Q_i=\operatorname{logit}(AUC_i),\qquad U_i=\widehat{\operatorname{Var}}(Q_i). ]

Rubin 总方差为:

[ T=\bar U+\left(1+\frac{1}{m}\right)B, ]

其中 (B) 是各个 (Q_i) 的样本方差。可以用 mice::pool.scalar() 合并:

pooled <- mice::pool.scalar(
  Q = Q,
  U = U,
  n = Inf
)

se <- sqrt(pooled$t)
critical <- stats::qt(0.975, df = pooled$df)

result <- data.frame(
  pooled_AUC = plogis(pooled$qbar),
  lower_95 = plogis(pooled$qbar - critical * se),
  upper_95 = plogis(pooled$qbar + critical * se)
)

result

n = Inf 将完整数据自由度近似为无限大,插补带来的有限自由度仍由合并步骤处理。这一设置并不代表聚类数量足够多;组数很少时,bootstrap 和区间估计都需要更谨慎。

这种合并要求插补模型合理、各次插补估计同一个目标,并满足 logit 尺度上的近似正态性。反变换后的合并 AUC 一般不等于 mean(A)。

结果应如何解释

上述计算评价模型在拟合数据上的总体区分能力,包含估计的随机效应,ROC 中也包含跨组的病例与非病例比较。如果研究目标是判断同一组内能否区分个体,需要另行定义组内 AUC。

解释和使用结果时,还需注意:

  • mean(A) 可以作为描述性汇总,但各次插补 AUC 的标准差不能直接当作合并标准误。
  • 不要先合并回归系数再计算一次 AUC,也不要直接平均各次插补的置信区间端点。
  • 如果要评价对新组的预测能力,应按组划分训练集和验证集,并采用与实际应用一致的预测方式。
  • 交叉验证中的插补和模型拟合都应限制在相应训练折内,避免验证信息泄漏。如果结局也做了插补,不能把插补出的验证结局直接当作真实标签。
  • 上述 bootstrap 估计的是表观 AUC 的方差,没有校正乐观偏倚。评价泛化性能时,需要另行设计验证或偏倚校正流程,并让插补、重采样和合并方法相互匹配。

备注:内容仅供参考。