ROC及曲线面积汇总学习

ROC基础

生成模拟数据

率的计算

R语言计算测试

ROCR：

pROC

ROC绘制

单个ROC

两个ROC

Logistic回归的ROC曲线

timeROC

ROC基础

ROC曲线的横坐标是假阳性率，纵坐标是真阳性率，需要的结果是这个率表示疾病阳性的率（而不能误算为阴性的率）。

ROC阳性结果还是阴性结果？_roc曲线计算阳性阴性预测值-CSDN博客

生成模拟数据

rm(list = ls())
set.seed(1234)
#模拟数据
ca125 <- c(rnorm(10,80,20),rnorm(20,50,10))
##前10个均值为80，标准差为20的正态分布数据，后20个均值为50，标准差为10的正态分布随机数
group <- c(rep(c("肿瘤","非肿瘤"),c(10,20)))#
data <- data.frame(ca125,group)

head(data)

      ca125 group
1  55.85869  肿瘤
2  85.54858  肿瘤
3 101.68882  肿瘤
4  33.08605  肿瘤
5  88.58249  肿瘤
6  90.12112  肿瘤

率的计算

#设置ca125>60，判断为肿瘤
data$pre <- ifelse(data$ca125 > 60,"猜他是肿瘤","猜他不是肿瘤")
#列联表
table(data$group,data$pre)猜他不是肿瘤 猜他是肿瘤
非肿瘤           19          1
肿瘤              2          8

根据这个四格表我们就能算出目前的真阳性率和假阳性率：

真阳性率：猜他是肿瘤猜对的人数 / 所有肿瘤人数
假阳性率：猜他是肿瘤猜错的人数 / 所有非肿瘤人数

真阳性率 = 8 / (2+8) = 0.8
假阳性率 = 1 / (19+1) = 0.05

如果目标结果是确定非肿瘤患者：

假阳性率：猜他不是肿瘤猜错的人数 / 所有肿瘤人数
真阳性率：猜他不是肿瘤猜对的人数 / 所有非肿瘤人数

假阳性率 = 2 / (2+8) = 0.2
真阳性率 = 19 / (19+1) = 0.95

R语言计算测试

R语言计算AUC（ROC曲线）的注意事项_r语言auc-CSDN博客

`ROCR：`

ROCR默认计算顺序靠后的类别的AUC。如果提供给labels的值是有序因子型变量，则排在前面的默认是阴性结果（negtive），排在后面的默认是阳性结果（positive），默认计算阳性结果（排序靠后）的AUC。如果是无序因子、数值、字符、逻辑型变量，会按照R语言的默认排序，比如按照数字大小、首字母顺序等，也是计算排序靠后的类别的AUC。


set.seed(20220840)
ca125_1 <- c(rnorm(10,80,20),rnorm(20,50,10))
ca125_2 <- c(rnorm(10,20,20),rnorm(20,70,10))
class=c(rep(1:0,c(10,20)))tumor <- c(rep(c("癌症","非癌症"),c(10,20)))df <- data.frame(`class`=class,`ca125_1`=ca125_1,`ca125_2`=ca125_2,`tumor`=tumor)
library(ROCR)pred <- prediction(predictions = ca125_1, # 预测指标labels = tumor # 真实结果)performance(pred, "auc")@y.values[[1]]
# [1] 0.075

这里我们想计算癌症的AUC，而不是非癌症的AUC，手动指定顺序！

pred <- prediction(predictions = ca125_1, # 预测指标labels = tumor # 真实结果,label.ordering = c("非癌症","癌症") # 此时就是计算癌症的AUC)performance(pred, "auc")@y.values[[1]]
## [1] 0.925

pROC

pROC包计算AUC也需要用来预测结果的指标以及真实结果。

这个包计算pROC略有不同，它是根据中位数来的，谁的中位数大，就计算谁的AUC，比如我们的这个例子，计算下中位数看看：

ca125_1

# 把ca125_1按照tumor的两个类别进行分组，然后分别计算中位数
tapply(ca125_1, tumor, median)
##     癌症   非癌症 
## 81.34426 49.99926（肿瘤的均值大）library(pROC)
roc(response=tumor, predictor=ca125_1)
#Area under the curve: 0.925（计算结果为癌症的）

ca125_2：是计算非癌症的AUC。

tapply(ca125_2, tumor, median)
##     癌症   非癌症 
## 13.52771 69.69272（非肿瘤的平均数大）roc(response=tumor, predictor=ca125_2)
Data: ca125_2 in 10 controls (tumor 癌症) < 20 cases (tumor 非癌症).
Area under the curve: 0.9

#需要设置levels和direction
# 此时计算的就是癌症的AUC
roc(response=tumor, predictor=ca125_2,levels=c("非癌症", "癌症"), # 这个顺序随便设定，重要的是directiondirection = "<" # 手动设定非癌症 < 癌症
)
#Data: ca125_2 in 20 controls (tumor 非癌症) < 10 cases (tumor 癌症).
#Area under the curve: 0.1

ROC绘制

各章示例代码/Chapter13 临床诊断实验评价.R · 杨敏迪/Analysis of Medical data by R language - 码云 - 开源中国 (gitee.com)

单个ROC

数据为动脉瘤性蛛网膜下腔出血患者的检测数据和预后：s100b是一个血清指标，outcome根据格拉斯哥评分分为good（4-5分：这里最大为5）、poor（1-3分）

rm(list = ls())
library(pROC)
data(aSAH)
roc1 <- roc(outcome ~ s100b, data = aSAH)
attributes(roc1)#查看结果包含内容
roc1$auc#
#Area under the curve: 0.7314

> #求约登指数
> roc.result <- data.frame(threshold = roc1$thresholds,
+                          sensitivity = roc1$sensitivities,
+                          specificity = roc1$specificities)
> View(roc.result)
> roc.result$youden <- roc.result$sensitivity + roc.result$specificity - 1
> head(roc.result)threshold sensitivity specificity      youden
1      -Inf   1.0000000  0.00000000  0.00000000
2     0.035   0.9756098  0.00000000 -0.02439024
3     0.045   0.9756098  0.06944444  0.04505420
4     0.055   0.9756098  0.11111111  0.08672087
5     0.065   0.9756098  0.13888889  0.11449864
6     0.075   0.9024390  0.22222222  0.12466125
> #找出约登指数最大的一行
> which.max(roc.result$youden)
[1] 18
> roc.result[18, ]threshold sensitivity specificity    youden
18     0.205   0.6341463   0.8055556 0.4397019

绘图

?plot.roc#查看参数细节
plot.roc(roc1, print.auc = TRUE, auc.polygon = TRUE,grid = c(0.1,0.2), grid.col = c("green","red"),auc.polygon.col = "lightblue", print.thres = TRUE)
#AUC的置信区间——DeLong法
ci.auc(roc1)

两个ROC

roc1 <- roc(aSAH$outcome, aSAH$s100b)
roc2 <- roc(aSAH$outcome, aSAH$ndka)
#DeLong非参数方法，Venkatraman回归模型法，bootstrap重抽样法
#默认为DeLong法
#默认为两组相关检测结果AUC的比较
#独立：参数paired = FALSE
roc.test(roc1,roc2)
#绘图
plot(roc1)
lines(roc2, col = "red")
test <- roc.test(roc1, roc2)
text(0.5,0.5, labels = paste("p-value = ",round(test$p.value, 3)))
legend("bottomright",legend = c("S100b", "NDKA"),col = c("1","red"), lwd = 2)

Logistic回归的ROC曲线

二分类变量

data = infert
#建立Logistic回归模型
fit <- glm(case ~ induced + spontaneous, family = binomial, data = infert)
library(epiDisplay)
logistic.display(fit)
#ROC绘制
lroc(fit, line.col = "red", lwd = 2)