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 的方差,没有校正乐观偏倚。评价泛化性能时,需要另行设计验证或偏倚校正流程,并让插补、重采样和合并方法相互匹配。
备注:内容仅供参考。
