R中Weibull参数生存模型测试集校准方法咨询

方法一：使用rms包的cph模型配合calibrate函数（推荐） #

rms包的calibrate函数直接支持测试集校准，但需要先使用rms的cph函数拟合Weibull参数化生存模型（替代基础包的survreg）：

library(rms)
library(survival)

# 拆分数据集
train <- lung[1:100,]
test <- lung[101:nrow(lung),]

# 用cph拟合Weibull分布的生存模型，开启surv参数以支持生存概率预测
cph_model <- cph(Surv(time, status) ~ sex, data = train, dist = "weibull", surv = TRUE)

# 选择校准的时间点（这里用训练集的事件中位时间作为参考）
median_event_time <- median(train$time[train$status == 1])

# 计算测试集校准结果，B为bootstrap重复次数（控制稳定性）
test_cal <- calibrate(cph_model, newdata = test, u = median_event_time, B = 100)

# 绘制校准曲线（对角线为完美校准，点越贴近对角线说明校准效果越好）
plot(test_cal)

方法二：手动计算校准曲线（兼容survreg模型） #

如果坚持使用基础包的survreg拟合模型，可以手动计算预测生存概率与实际观察概率的差异，绘制校准曲线：

library(survival)

train <- lung[1:100,]
test <- lung[101:nrow(lung),]

# 用survreg拟合Weibull模型
s <- survreg(Surv(time, status) ~ sex, data = train, dist = "weibull")
weibull_shape <- 1 / s$scale  # Weibull分布形状参数

# 选择校准时间点
cal_time <- median(train$time[train$status == 1])

# 预测测试集每个样本的生存概率
lp <- predict(s, newdata = test)
pred_surv_prob <- exp(-(cal_time / exp(lp))^weibull_shape)

# 按预测概率分组（这里分5组），计算每组平均预测概率与实际生存概率
pred_groups <- cut(pred_surv_prob, breaks = 5)
mean_pred <- tapply(pred_surv_prob, pred_groups, mean)

# 计算每组在cal_time时刻的实际生存概率
actual_surv_groups <- tapply(seq(nrow(test)), pred_groups, function(idx) {
  sf <- survfit(Surv(time, status) ~ 1, data = test[idx,])
  summary(sf, times = cal_time)$surv
})

# 绘制校准曲线
plot(mean_pred, unlist(actual_surv_groups), 
     xlab = "组平均预测生存概率", 
     ylab = "组实际生存概率",
     pch = 16, cex = 1.2)
abline(0, 1, col = "red", lwd = 2)  # 完美校准参考线

校准结果解读 #

• 校准曲线越贴近对角线（y=x），说明模型预测的生存概率与实际观察结果越一致，校准效果越好。
• 若曲线在对角线上方，说明模型低估了实际生存概率；若在下方，则高估了实际生存概率。

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