ROC曲线与模型评估¶

              预测阳性      预测阴性
实际阳性    TP（真阳性）   FN（假阴性）  ← 漏诊
实际阴性    FP（假阳性）   TN（真阴性）  ← 误诊

灵敏度 = TP / (TP + FN)  → 能抓住多少真正有病的人
特异度 = TN / (FP + TN)  → 能放过多少真正没病的人

ROC曲线就是：遍历所有可能的阈值，画出
  X轴 = 1 - 特异度（FPR，假阳性率）
  Y轴 = 灵敏度（TPR，真阳性率）

AUC = 0.5: 和抛硬币一样（没用的模型）
AUC = 0.7: 随机选一个病人和一个正常人，
           模型有70%的概率给病人打更高的分
AUC = 0.9: 90%的概率能正确区分
AUC = 1.0: 完美区分（现实中几乎不可能）

发文章参考：
  AUC < 0.7:  差
  0.7-0.8:    可以接受
  0.8-0.9:    好
  > 0.9:      优秀

# 安装和加载
install.packages("pROC")  # 安装pROC包
library(pROC)              # 加载

# 准备数据
# actual: 实际标签（0/1）
# predicted: 模型预测的概率（0-1之间的连续值）
actual <- c(1, 1, 1, 0, 0, 0, 1, 0, 1, 0)          # 真实标签
predicted <- c(0.9, 0.8, 0.7, 0.3, 0.2, 0.1, 0.6, 0.4, 0.65, 0.35)  # 预测概率

# 计算ROC
roc_obj <- roc(
  response = actual,      # 真实标签
  predictor = predicted,  # 预测概率
  levels = c(0, 1),       # 标签水平：0=阴性，1=阳性
  direction = "<"         # 预测值越大越可能是阳性
)

# 查看AUC
auc(roc_obj)  # 打印AUC值

# 画ROC曲线
plot(roc_obj,
     print.auc = TRUE,           # 显示AUC值
     print.auc.x = 0.4,          # AUC文字X位置
     print.auc.y = 0.3,          # AUC文字Y位置
     col = "blue",               # 曲线颜色
     lwd = 2,                    # 线宽
     main = "ROC Curve")         # 标题
abline(a = 0, b = 1, lty = 2, col = "gray")  # 对角参考线

# 计算AUC的95%置信区间
ci_auc <- ci.auc(
  roc_obj,
  conf.level = 0.95,   # 95%置信水平
  method = "delong"     # DeLong方法（推荐）
)
print(ci_auc)  # 如 0.85 (0.78-0.92)

# Bootstrap方法（样本小时更稳定）
ci_boot <- ci.auc(roc_obj, method = "bootstrap", boot.n = 2000)
print(ci_boot)

# 方法1：Youden Index（灵敏度+特异度-1 最大）
best_youden <- coords(
  roc_obj,
  x = "best",                 # 最佳切点
  best.method = "youden",     # Youden Index
  ret = c("threshold", "sensitivity", "specificity")  # 返回值
)
print(best_youden)

# 方法2：最近左上角点
best_closest <- coords(
  roc_obj,
  x = "best",
  best.method = "closest.topleft"  # 离(0,1)最近的点
)
print(best_closest)

# 方法3：指定灵敏度要求
coords(roc_obj, x = 0.9, input = "sensitivity",
       ret = c("threshold", "specificity"))  # 灵敏度≥90%时的阈值

# 两个模型的预测结果
pred_model1 <- predict(model1, type = "response")  # 模型1预测
pred_model2 <- predict(model2, type = "response")  # 模型2预测

# 分别画ROC
roc1 <- roc(actual, pred_model1)
roc2 <- roc(actual, pred_model2)

# DeLong检验比较两个AUC
test_result <- roc.test(roc1, roc2, method = "delong")
print(test_result)
# p < 0.05 说明两个模型的AUC有显著差异

# 画在一张图上
plot(roc1, col = "blue", lwd = 2, main = "ROC Comparison")
plot(roc2, col = "red", lwd = 2, add = TRUE)  # add=TRUE叠加画
legend("bottomright",
       legend = c(paste("Model 1: AUC =", round(auc(roc1), 3)),
                  paste("Model 2: AUC =", round(auc(roc2), 3))),
       col = c("blue", "red"), lwd = 2)

# 方法1：用rms包的校准（推荐）
library(rms)
dd <- datadist(data)
options(datadist = "dd")
fit <- lrm(outcome ~ var1 + var2, data = data, x = TRUE, y = TRUE)

cal <- calibrate(fit, B = 1000, method = "boot")  # Bootstrap校准
plot(cal,
     xlab = "Predicted Probability",
     ylab = "Observed Probability")
# 45度对角线 = 完美校准
# 曲线偏上 = 模型低估了风险
# 曲线偏下 = 模型高估了风险

# 方法2：手动画校准曲线
library(ggplot2)

# 将预测概率分成10组（十分位数）
pred <- predict(model, type = "response")
cal_data <- data.frame(pred = pred, actual = actual)
cal_data$group <- cut(pred, breaks = quantile(pred, probs = seq(0, 1, 0.1)),
                      include.lowest = TRUE)

# 每组的平均预测值和实际事件率
cal_summary <- aggregate(cbind(pred, actual) ~ group, data = cal_data, mean)

ggplot(cal_summary, aes(x = pred, y = actual)) +
  geom_point(size = 3) +                         # 散点
  geom_abline(slope = 1, intercept = 0, lty = 2) +  # 对角线
  geom_smooth(method = "loess", se = FALSE) +     # 平滑曲线
  labs(x = "Predicted", y = "Observed",
       title = "Calibration Plot") +
  theme_minimal()

# 定量评估校准度
library(ResourceSelection)
hl_test <- hoslem.test(actual, pred, g = 10)  # 分10组
print(hl_test)
# p > 0.05 表示校准良好（预测和实际无显著差异）
# p < 0.05 表示校准不良

# 生存分析中的ROC：不同时间点有不同的AUC
library(timeROC)

# 计算时间依赖AUC
time_roc <- timeROC(
  T = clinical$time,        # 生存时间
  delta = clinical$status,  # 事件状态（0/1）
  marker = risk_score,      # 预测标记（如风险评分）
  cause = 1,                # 事件类别
  times = c(365, 365*3, 365*5),  # 评估时间点（1年、3年、5年）
  iid = TRUE                # 计算置信区间
)

# 查看各时间点AUC
time_roc$AUC  # AUC值
confint(time_roc)  # 置信区间

# 画时间依赖ROC
plot(time_roc, time = 365*3,   # 画3年的ROC
     title = FALSE, lwd = 2)
legend("bottomright",
       paste("3-Year AUC =", round(time_roc$AUC[2], 3)))

# 一张图展示多时间点ROC
library(ggplot2)

# 提取数据画图
roc_data <- data.frame(
  FPR_1yr = 1 - time_roc$SP[, 1],
  TPR_1yr = time_roc$TP[, 1],
  FPR_3yr = 1 - time_roc$SP[, 2],
  TPR_3yr = time_roc$TP[, 2],
  FPR_5yr = 1 - time_roc$SP[, 3],
  TPR_5yr = time_roc$TP[, 3]
)

ggplot() +
  geom_line(data = roc_data, aes(FPR_1yr, TPR_1yr, color = "1-Year")) +
  geom_line(data = roc_data, aes(FPR_3yr, TPR_3yr, color = "3-Year")) +
  geom_line(data = roc_data, aes(FPR_5yr, TPR_5yr, color = "5-Year")) +
  geom_abline(slope = 1, intercept = 0, lty = 2) +
  labs(x = "1 - Specificity", y = "Sensitivity",
       title = "Time-Dependent ROC Curves") +
  scale_color_manual(values = c("red", "blue", "green")) +
  theme_minimal()

# ROC分析流程
roc() → auc() → ci.auc() → coords("best") → plot()

# 模型比较
roc.test(roc1, roc2, method="delong")

# 校准分析
calibrate() 或 hoslem.test()

# 时间依赖ROC
timeROC() → plot() → confint()

# 关键R包
pROC      — ROC曲线和AUC
timeROC   — 时间依赖ROC
rms       — 校准曲线
ResourceSelection — Hosmer-Lemeshow检验
ggDCA     — DCA决策曲线

# AUC评价标准
< 0.7:   差（不适合临床使用）
0.7-0.8: 可接受
0.8-0.9: 好
> 0.9:   优秀

# 最佳切点方法
Youden Index:        max(灵敏度+特异度-1)
Closest to topleft:  离(0,1)最近的点
Cost-based:          考虑漏诊/误诊成本