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

2022-03-29 09:37:36 浏览数 (1)

导语

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包介绍

01

R包pROC

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

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

1.

(1)建立ROC对象并计算AUC

代码语言:javascript复制
roc1 <- roc(aSAH$outcome, aSAH$s100b)
print(roc1)
#或
roc(outcome ~ s100b, aSAH)
代码语言:javascript复制
#计算部分auc
auc(roc1, partial.auc = c(1, .9))

(2)使ROC曲线平滑

代码语言:javascript复制
smooth(roc1)

(3)方差计算

代码语言:javascript复制
roc2 <- roc(aSAH$outcome, aSAH$wfns)
roc3 <- roc(aSAH$outcome, aSAH$ndka)
var(roc1) #计算ROC的方差
cov(roc1, roc3) #计算两组ROC的协方差

(4)绘制ROC

①plot绘制ROC

代码语言:javascript复制
plot(roc1)

②绘制ROC并计算pROC

代码语言:javascript复制
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曲线

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

2. 绘制置信区间

(1)计算置信区间

代码语言:javascript复制
#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)绘制置信区间

代码语言:javascript复制
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曲线

代码语言:javascript复制
plot(roc2, add=TRUE) #加入roc2
plot(ci.thresholds(roc2)) #计算特异性和敏感性CI阈值

3. 比较ROC曲线

(1)Bootstrap比较ROC曲线

代码语言:javascript复制
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)修改部分参数

代码语言:javascript复制
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)限定参数

代码语言:javascript复制
power.roc.test(ncases=41, ncontrols=72, sig.level=0.05, power=0.95)
#ncases,观察样本数
#ncontrols,对照样本数
#sig.level,期望的显著性水平(第一类错误的probability)
#power,测试的期望power(第二类错误的1 -probability)

02

R包plotROC

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

代码语言:javascript复制
install.packages("plotROC")
library(plotROC)

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

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

代码语言:javascript复制
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曲线

代码语言:javascript复制
ggroc2 <- ggplot(rocdata, aes(m = M, d = D, color = Z))   geom_roc()

3. calc_auc计算ROC曲线下面积

代码语言:javascript复制
calc_auc(ggroc2)

4. direct_label添加标签

代码语言:javascript复制
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

代码语言:javascript复制
ggplot(rocdata, aes(m = M, d = D))   geom_roc(n.cuts = 20)
代码语言:javascript复制
ggplot(rocdata, aes(m = M, d = D))   geom_roc(cutoffs.at = c(1.5, 1, .5, 0, -.5))

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

(1)置信区间

代码语言:javascript复制
ggplot(rocdata, aes(m = M, d = D))   geom_roc()   geom_rocci()

(2)置信区间的显著性水平

代码语言:javascript复制
ggplot(rocdata, aes(m = M, d = D, color = Z))   geom_roc()   geom_rocci(sig.level = .01)
#sig.level置信区间的显著性水平

(3)biomarker处展示置信区间

代码语言:javascript复制
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()展示置信区间

代码语言:javascript复制
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曲线

代码语言:javascript复制
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曲线

代码语言:javascript复制
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图添加指南和注释

代码语言:javascript复制
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)

03

R包ggROC

代码语言:javascript复制
install.packages("ggROC")
library(ggROC)
data(roc)

直接输出PDF文件

代码语言:javascript复制
ggroc(roc, 0.01,"green",sp =19,output="F:/roc.pdf")

04

R包ROCR

代码语言:javascript复制
install.packages("ROCR")
library(ROCR)
data(ROCR.simple)
pred <- prediction( ROCR.simple$predictions, ROCR.simple$labels)
perf <- performance(pred,"tpr","fpr")
plot(perf)

小编总结

ROC曲线在生信分析中是一个非常常见的分析工具,一般用在我们对自己构建的方法模型的进行验证分析的时候,有一点要注意的是,前提需要有金标准做对照。今天介绍了四种ROC曲线的绘制方法,大家可以试试哦!

0 人点赞