如何利用pROC包获取ROC曲线各阈值的似然比及Bootstrap 95%置信区间

API 探路者

2026-7-6

计算诊断阈值对应的似然比及Bootstrap 95%置信区间

需求说明

已通过pROC结合Bootstrap法选择了最优阈值（closest top left方法），并能获取对应90%、95%灵敏度的阈值及其灵敏度、特异度，现需计算这些阈值对应的阳性似然比（LR+）和阴性似然比（LR-），以及它们的Bootstrap 95%置信区间。

实现步骤与代码

1. 准备数据（可复现代码）

set.seed(42)

num_subjects <- 40

continuous_variable <- rnorm(num_subjects, mean = 50, sd = 10)

probability_of_one <- 1 / (1 + exp(-0.1 * (continuous_variable - 50)))
binary_variable <- rbinom(num_subjects, 1, probability_of_one)

df <- data.frame(
  continuous_variable = continuous_variable,
  binary_variable = binary_variable
)

2. 拟合ROC曲线并获取目标阈值

先拟合ROC曲线，再提取三类目标阈值：

最优阈值（closest top left）
灵敏度90%对应的阈值
灵敏度95%对应的阈值

library(pROC)

# 拟合ROC曲线
ROC_1 <- roc(response = df$binary_variable, predictor = df$continuous_variable, ci = TRUE)

# 获取最优阈值（closest top left）
best_thresh <- coords(ROC_1, x = "best", best.method = "closest.topleft")$threshold

# 获取灵敏度90%对应的阈值（取最接近90%的点）
sens90_thresh <- coords(ROC_1, x = 0.9, input = "sensitivity", ret = "threshold")$threshold

# 获取灵敏度95%对应的阈值
sens95_thresh <- coords(ROC_1, x = 0.95, input = "sensitivity", ret = "threshold")$threshold

# 整理所有目标阈值
target_thresholds <- c("最优阈值" = best_thresh, "90%灵敏度阈值" = sens90_thresh, "95%灵敏度阈值" = sens95_thresh)

3. 计算原始数据的似然比

基于原始ROC结果，通过灵敏度、特异度推导似然比：

阳性似然比（LR+）= 灵敏度 / (1 - 特异度)
阴性似然比（LR-）= (1 - 灵敏度) / 特异度

# 定义计算似然比的函数
calc_lr <- function(roc_obj, threshold) {
  coords_res <- coords(roc_obj, x = threshold, ret = c("sensitivity", "specificity"))
  lr_plus <- coords_res$sensitivity / (1 - coords_res$specificity)
  lr_minus <- (1 - coords_res$sensitivity) / coords_res$specificity
  return(data.frame(
    threshold = threshold,
    sensitivity = coords_res$sensitivity,
    specificity = coords_res$specificity,
    lr_plus = lr_plus,
    lr_minus = lr_minus
  ))
}

# 计算所有目标阈值的似然比
original_lr <- do.call(rbind, lapply(target_thresholds, calc_lr, roc_obj = ROC_1))
rownames(original_lr) <- names(target_thresholds)
print("原始数据的似然比结果：")
print(original_lr)

4. Bootstrap法计算似然比的95%置信区间

由于pROC的ci.coords不直接支持似然比的置信区间，自定义Bootstrap抽样流程实现：

# 设置Bootstrap参数
boot_n <- 2000
set.seed(42) # 保证结果可重复

# 自定义Bootstrap函数，每次抽样后计算对应阈值的似然比
boot_lr <- function(data, thresholds) {
  # 有放回抽样
  sample_idx <- sample(nrow(data), replace = TRUE)
  sample_data <- data[sample_idx, ]
  
  # 拟合ROC曲线（处理可能的极端情况，比如样本全为0/1）
  tryCatch({
    roc_sample <- roc(response = sample_data$binary_variable, predictor = sample_data$continuous_variable)
    # 计算每个阈值的似然比
    lr_list <- lapply(thresholds, function(thresh) {
      coords_sample <- coords(roc_sample, x = thresh, ret = c("sensitivity", "specificity"))
      lr_plus <- coords_sample$sensitivity / (1 - coords_sample$specificity)
      lr_minus <- (1 - coords_sample$sensitivity) / coords_sample$specificity
      return(c(lr_plus, lr_minus))
    })
    return(unlist(lr_list))
  }, error = function(e) {
    # 若抽样后无法拟合ROC（如结局全为一类），返回NA
    return(rep(NA, 2 * length(thresholds)))
  })
}

# 执行Bootstrap抽样
boot_results <- replicate(boot_n, boot_lr(df, target_thresholds))

# 整理结果并计算95%置信区间（去掉NA值）
boot_ci <- data.frame()
for (i in seq_along(target_thresholds)) {
  lr_plus_vals <- boot_results[2*i-1, ]
  lr_plus_vals <- lr_plus_vals[!is.na(lr_plus_vals)]
  lr_minus_vals <- boot_results[2*i, ]
  lr_minus_vals <- lr_minus_vals[!is.na(lr_minus_vals)]
  
  boot_ci <- rbind(boot_ci, data.frame(
    threshold_name = names(target_thresholds)[i],
    threshold = target_thresholds[i],
    lr_plus_lower = quantile(lr_plus_vals, 0.025),
    lr_plus_upper = quantile(lr_plus_vals, 0.975),
    lr_minus_lower = quantile(lr_minus_vals, 0.025),
    lr_minus_upper = quantile(lr_minus_vals, 0.975)
  ))
}

# 合并原始似然比和Bootstrap置信区间
final_results <- merge(original_lr, boot_ci, by.x = "row.names", by.y = "threshold_name")
colnames(final_results)[1] <- "阈值类型"
print("最终结果（含Bootstrap 95%置信区间）：")
print(final_results)

关键说明

若Bootstrap抽样中出现极端样本（如结局全为0或1），会返回NA，计算置信区间时需剔除这些无效值
阳性似然比越大，说明该阈值区分阳性病例的能力越强；阴性似然比越小，排除阴性病例的能力越强

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