ROC曲线不用愁，四种R包教你一步搞定！

作图丫

发布于 2022-03-29 09:37:36

6.6K0

发布于 2022-03-29 09:37:36

文章被收录于专栏：作图丫作图丫

导语

GUIDE ╲

前面我们介绍了一个对有害同义突变预测的方法PrDSM，可以发现，在对模型的分析中，大量的使用ROC对模型进行评估，今天我们就来介绍一下ROC的相关内容和两种ROC绘图方法：pROC、plotROC、ggROC和ROCR。

ROC介绍

ROC曲线是受试者工作特征曲线 / 接收器操作特性曲线(receiver operating characteristic curve), 是一个反映二元分类器系统在其识别阈值变化时的诊断能力的图形。

ROC曲线是通过绘制真阳性率(TPR)与假阳性率(FPR)在不同阈值设置下的曲线。在机器学习中，真阳性率也被称为灵敏度、回忆率或检出率。假阳性率也称为误报率，可以计算为(1 -特异度)。ROC曲线也可以被认为是决策规则的Type I Error 的函数（当性能仅从总体的一个样本中计算时，它可以被认为是这些量的估计值）。因此ROC曲线是敏感度或召回率作为降噪的函数。一般情况下，如果真阳性率和假阳性率分布已知，可以通过对y轴上的真阳性率和x轴上的假阳性率绘制的累积分布函数(概率分布下的面积，从-∞到判别阈值)来生成ROC曲线，因此ROC图有时被称为敏感性vs(1−特异性)图。

考虑一个两类预测问题(二元分类)，其中结果被标记为正(p)或负(n)。一个二元分类器有四种可能的结果。①如果预测的结果是p，实际值也是p，则称为真正(true positive, TP)。②如果预测的结果是p，实际值为n，则称为假阳性(FP)。③当预测结果与实际值均为n时，是真阴性(TN)。④当预测结果为n而实际值为p时，是假阴性(FN)。

下图所示各个指标及计算公式:

最好的预测方法是在ROC空间的左上角或坐标(0,1)处找到一个点，表示100%的敏感性(无假阴性)和100%的特异性(无假阳性)。(0,1)点也被称为完美分类。所以ROC曲线越靠近左上角，说明该方法分类效果越好。最靠近左上角的ROC曲线上的点是分类错误最少的最好阈值，其假正例和假反例总数最少。可以对不同的学习器比较性能。将各个学习器的ROC曲线绘制到同一坐标中，直观地鉴别优劣，靠近左上角的ROC曲所代表的学习器准确性最高。

AUC是衡量学习器优劣的一种性能指标，为ROC曲线下与坐标轴围成的面积。其意义是：①因为是在1x1的方格里求面积，AUC必在0~1之间。②假设阈值以上是阳性，以下是阴性；③若随机抽取一个阳性样本和一个阴性样本，分类器正确判断阳性样本的值高于阴性样本的概率 = AUC 。④AUC值越大的分类器，正确率越高。

R包介绍

R包pROC

pROC是一个用于显示、平滑和比较ROC曲线的工具。(部分)曲线下面积AUC(pAUC)可以通过基于U-statistics或bootstrap的统计检验进行比较。可以计算(p)AUC或ROC曲线的置信区间。

install.packages("pROC")
library(pROC)
data(aSAH)
#该数据集总结了113例动脉瘤性蛛网膜下腔出血的临床和实验室变量。

（1）建立ROC对象并计算AUC

roc1 <- roc(aSAH$outcome, aSAH$s100b)
print(roc1)
#或
roc(outcome ~ s100b, aSAH)

#计算部分auc
auc(roc1, partial.auc = c(1, .9))

（2）使ROC曲线平滑

smooth(roc1)

（3）方差计算

roc2 <- roc(aSAH$outcome, aSAH$wfns)
roc3 <- roc(aSAH$outcome, aSAH$ndka)
var(roc1) #计算ROC的方差
cov(roc1, roc3) #计算两组ROC的协方差

（4）绘制ROC

①plot绘制ROC

plot(roc1)

②绘制ROC并计算pROC

roc4 <- roc(aSAH$outcome,aSAH$s100b, 
            percent=TRUE,
            #percent敏感性、特异性和AUC是否必须用百分数或分数表示
            partial.auc=c(100, 90),
            #计算部分auc
            partial.auc.correct=TRUE,
            #是否对AUC校正，以获得最大AUC为1.0和非歧视性AUC 0.5，
            partial.auc.focus="sens",
            #计算partial.auc根据"sens"，即敏感性
            ci=TRUE, boot.n=100,
            #ci置信区间，boot.n是bootstrap重复数
            ci.alpha=0.9, stratified=FALSE,
            #stratified表示bootstrap是否分层
            # arguments for plot
            plot=TRUE, auc.polygon=TRUE,
            #auc.polygon是否将area显示为多边形
            max.auc.polygon=TRUE, grid=TRUE,
            #max.auc.polygon是否将最大可能的区域显示为多边形
            #grid是否添加背景网格
            print.auc=TRUE, show.thres=TRUE
            #print.auc，AUC的数值是否应该打印在图上
            )

③添加ROC曲线

roc5 <- roc(aSAH$outcome, aSAH$wfns,
            plot=TRUE, add=TRUE, percent=roc4$percent)
#在上述ROC绘图基础上再绘制
#add是否将其他ROC曲线将被添加到现有的plot中

2. 绘制置信区间

（1）计算置信区间

#ROC曲线的坐标系
coords(roc1, "best", ret=c("threshold", "specificity", "1-npv"))
coords(roc2, "local maximas", ret=c("threshold", "sens", "spec", "ppv", "npv"))
#第二个参数为要找的坐标。“all”:ROC曲线上的所有点。
#“local maximas”:ROC曲线的局部极大值。“best”:用一个方法确定最佳阈值
#ret返回坐标
ci(roc2)  #计算置信区间

（2）绘制置信区间

sens.ci <- ci.se(roc1, specificities=seq(0, 1, .1))
#ci.se,在特定情况下计算灵敏度的置信区间
plot(sens.ci, type="shape", col="lightblue")
plot(sens.ci, type="bars")
#type，plot的形状

（3）加入ROC曲线

plot(roc2, add=TRUE) #加入roc2
plot(ci.thresholds(roc2)) #计算特异性和敏感性CI阈值

3. 比较ROC曲线

（1）Bootstrap比较ROC曲线

roc.test(roc1, roc2, reuse.auc=FALSE, partial.auc=c(1, .8),
         #reuse.auc=FALSE(默认值)，“roc”对象包含一个“auc”字段
         #请在测试中重用这些specifications。
         allow.invalid.partial.auc.correct = TRUE,
         #当试图校正对角线以下的pAUC时，指示校正是否必须返回NA(带有警告)。
         #设置为TRUE则返回一个(可能无效)已更正的AUC。
         #这对于在bootstrap操作中避免对低pAUCs造成偏倚尤其有用
         partial.auc.focus="se", partial.auc.correct=TRUE)

（2）修改部分参数

roc.test(roc1, roc2, reuse.auc=FALSE, partial.auc=c(1, .8),
         allow.invalid.partial.auc.correct = TRUE,
         partial.auc.correct=TRUE,
         boot.n=1000, boot.stratified=FALSE)

4. ROC曲线的样本量power计算

计算ROC曲线的样本量、power、显著性水平或最小曲线下面积

（1）一条曲线

（2）两条曲线

（3）限定参数

power.roc.test(ncases=41, ncontrols=72, sig.level=0.05, power=0.95)
#ncases，观察样本数
#ncontrols,对照样本数
#sig.level，期望的显著性水平（第一类错误的probability）
#power，测试的期望power(第二类错误的1 -probability)

R包plotROC

大多数ROC曲线绘图模糊了cutoff 值，限制了多条曲线的解释和比较。plotROC试图通过提供绘图和交互式工具来解决这些缺点。提供可以生成用于web使用的交互式ROC曲线图，以及打印版本的功能。plotROC是基于ggplot2绘图的。

install.packages("plotROC")
library(plotROC)

1.生成展示绘图的ROC数据

以下一段语句是用于生成展示绘图的ROC数据（在我们自己绘图的时候可用自己的数据）

D.ex <- rbinom(50, 1, .5)
#rbinom生成一组二项分布随机数，参数为（生成的随机数数量，进行随机试验的次数，二项分布概率）
rocdata <- data.frame(D = c(D.ex, D.ex),
    #D是分类标签，必须为0和1。
           M = c(rnorm(50, mean = D.ex, sd = .4),
           #rnorm是生成50个服从正态分布的随机数，均值是D.ex，方差是0.4
           rnorm(50, mean = D.ex, sd = 1)),
    #M是连续marker值或预测的类标签
           Z = c(rep("A", 50), rep("B", 50))
#z是分组参数，指的是A、B两种分类
    )
rocdata

2. ggplot绘制ROC曲线

ggroc2 <- ggplot(rocdata, aes(m = M, d = D, color = Z)) + geom_roc()

3. calc_auc计算ROC曲线下面积

calc_auc(ggroc2)

4. direct_label添加标签

direct_label(ggroc2, labels = "Biomarker",
             label.angle = 45,
             nudge_x=0.05,nudge_y = -0.1)
#labels标签向量直接添加到图中
#label.angle调整标签角度
#nudge_x, nudge_y水平和垂直的调整，以推动标签。
#Size标签大小

5. 更改cutoffs

ggplot(rocdata, aes(m = M, d = D)) + geom_roc(n.cuts = 20)

ggplot(rocdata, aes(m = M, d = D)) + geom_roc(cutoffs.at = c(1.5, 1, .5, 0, -.5))

6. geom_rocci添加ROC曲线的置信区间

（1）置信区间

ggplot(rocdata, aes(m = M, d = D)) + geom_roc() + geom_rocci()

（2）置信区间的显著性水平

ggplot(rocdata, aes(m = M, d = D, color = Z)) + geom_roc() + geom_rocci(sig.level = .01)
#sig.level置信区间的显著性水平

（3）biomarker处展示置信区间

ggplot(rocdata, aes(m = M, d = D)) + geom_roc(n.cuts = 0) +
  geom_rocci(ci.at = quantile(rocdata$M, c(.1, .25, .5, .75, .9)),linetype = 1)
#在biomarker处展示置信区间

（4）stat_rocci()展示置信区间

ggplot(rocdata, aes(m = M, d = D)) + geom_roc() +
  stat_rocci(ci.at = quantile(rocdata$M, c(.1, .3, .5, .7, .9)))

7. plot_interactive_roc生成交互式的ROC曲线

rocdata <- data.frame(D = c(D.ex, D.ex),
                      M = c(rnorm(50, mean = D.ex, sd = .4), rnorm(50, mean = D.ex, sd = 1)),
                      Z = c(rep("A", 50), rep("B", 50)))
rocplot <- ggplot(rocdata, aes(m = M, d = D)) + geom_roc()
plot_interactive_roc(rocplot + aes(color = Z))

8. stat_roc()计算经验ROC曲线

D.ex <- rbinom(50, 1, .5)
rocdata <- data.frame(D = c(D.ex, D.ex),
                      M = c(rnorm(50, mean = D.ex, sd = .4), rnorm(50, mean = D.ex, sd = 1)),
                      Z = c(rep("A", 50), rep("B", 50)))
ggplot(rocdata, aes(m = M, d = D)) + stat_roc()

9. style_roc()在ROC图添加指南和注释

D.ex <- rbinom(50, 1, .5)
fakedata <- data.frame(M1 = rnorm(50, mean = D.ex),
                       D = D.ex)
ggplot(fakedata, aes(m = M1, d = D)) + geom_roc() +
  style_roc(xlab = "1 - Specificity",theme = theme_grey)

R包ggROC

install.packages("ggROC")
library(ggROC)
data(roc)

直接输出PDF文件

ggroc(roc, 0.01,"green",sp =19,output="F:/roc.pdf")

R包ROCR

install.packages("ROCR")
library(ROCR)
data(ROCR.simple)
pred <- prediction( ROCR.simple$predictions, ROCR.simple$labels)
perf <- performance(pred,"tpr","fpr")
plot(perf)