如何正确结合R的MICE多重插补与RMS包绘制限制性立方样条图？

你当前采用的「合并所有插补数据集+拟合单一模型」的做法不符合多重插补的统计原则，会导致标准误估计偏误，无法正确反映缺失数据带来的不确定性。正确的做法是在每个插补数据集上分别拟合模型，再用Rubin法则合并各模型的预测结果。

正确实现步骤与代码 #

1. 加载包并生成测试数据 #

library(rms)
library(ggplot2)
library(mice)
library(data.table)

# 生成测试数据
set.seed(17)
n <- 1000
age <- rnorm(n, 50, 10)
time_at_risk <- rnorm(n, 200, 80)
cholesterol <- rnorm(n, 200, 25)
sex <- factor(sample(c('female','male'), n, TRUE))
L <- .4*(sex=='male') + .045*(age-50) + (log(cholesterol - 10)-5.2)*(-2*(sex=='female') + 2*(sex=='male'))
outcome <- ifelse(runif(n) < plogis(L), 1, 0)
cholesterol[1:3] <- NA
test_data <- data.frame(outcome=outcome, age=age, sex=sex, cholesterol=cholesterol, time_at_risk=time_at_risk)

2. 多重插补与模型拟合 #

用mice::with()在每个插补数据集上拟合Cox模型，并生成限制性立方样条的预测结果：

# 多重插补
test_mice_output <- mice(test_data, m=10, maxit=50, seed=10)

# 在每个插补数据集上拟合模型并生成预测
fit_list <- with(test_mice_output, {
  dd <- datadist(.data); options(datadist='dd')
  cox_model <- cph(Surv(time_at_risk, outcome) ~ rcs(cholesterol, 3) + age + sex, 
                   x=TRUE, y=TRUE, ref.zero=TRUE)
  Predict(cox_model, cholesterol, ref.zero=TRUE, fun=exp)
})

3. 合并插补结果（Rubin法则） #

手动合并多个插补的预测值，计算合并后的点估计、标准误和置信区间：

# 提取所有插补的预测数据
pred_list <- lapply(fit_list$analyses, function(x) as.data.table(x))

# 合并数据并计算均值（点估计）、标准误
combined_pred <- rbindlist(pred_list, idcol='imputation')[, .(
  or_mean = mean(yhat),
  or_se = sqrt(mean(se.yhat^2) + (1 + 1/10) * var(yhat))  # Rubin法则计算合并标准误
), by = .(cholesterol)]

# 计算95%置信区间
combined_pred[, `:=`(
  or_lower = or_mean - 1.96*or_se,
  or_upper = or_mean + 1.96*or_se
)]

4. 绘制限制性立方样条图 #

ggplot(combined_pred, aes(x=cholesterol, y=or_mean)) +
  geom_line(color='blue', linewidth=1) +
  geom_ribbon(aes(ymin=or_lower, ymax=or_upper), alpha=0.2, fill='blue') +
  ylim(c(0.4, 1.8)) +
  xlim(c(140, 260)) +
  ylab("预测风险比（OR）") +
  xlab("胆固醇") +
  theme_bw()

关键说明 #

• mice::with()会自动遍历每个插补数据集，避免手动拆分数据的麻烦
• 合并预测值时，遵循Rubin法则：合并标准误 = 根号(平均内插补方差 + (1+1/m)×插补间方差)，其中m是插补次数
• 最终的图同时展示了点估计和考虑缺失数据不确定性的置信区间，结果更可靠

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