「R」ROC三剑客(二)分析与可视化ROC——plotROC、pROC

2020-05-14 26160动手实验室

《使用R语言手撕ROC曲线》这篇文章中我讲了ROC曲线的本质以及如何计算和绘制ROC曲线。注意,我这里谈到的ROC并未曾涉及机器学习模型的拟合与预测,而是指存在一组真实的连续型数值数据设定阈值的不同对响应变量(二分类)的影响(真阳性率、假阳性率)。

这一篇文章我们学习两个跟ROC相关的R包:

  • plotROC - Generate ROC Curve Charts for Print and Interactive Use
  • pROC - display and analyze ROC curves in R and S+

plotROC

plotROC包较为简单与单一,它就是用来绘制ROC曲线的,包中定义的函数基于ggplot2,因此我们可以结合ggplot2使用和修改、美化图形结果。

从GitHub上安装

devtools::install_github("hadley/ggplot2")
devtools::install_github("sachsmc/plotROC") library(plotROC)

从CRAN

install.packages("plotROC")

快速使用

plotROC提供了Shiny应用,只需要键入

shiny_plotROC()

即可通过图形界面使用。
咱们还是来看命令吧,要有点难度不是?

命令行使用

导入包与创建模拟数据:

library(plotROC)
set.seed(2529)
D.ex <- rbinom(200, size = 1, prob = .5)
M1 <- rnorm(200, mean = D.ex, sd = .65)
M2 <- rnorm(200, mean = D.ex, sd = 1.5)
test <- data.frame(D = D.ex, D.str = c("Healthy", "Ill")[D.ex + 1],                    
M1 = M1, M2 = M2, stringsAsFactors = FALSE)

简单绘图:

basicplot <- ggplot(test, aes(d = D, m = M1)) + geom_roc()
basicplot

1-1.png

这里我们唯一需要理清的是d与m映射是什么,现在我们查看下生成的数据框:
1-2.png

上述画图只使用到了D与M1,只关注这两列即可。D是一个0-1列,即表示结果的两分类信息,M1是一个数值型数据。我们可以姑且称d为decision缩写,m为measurement缩写。

一旦我们理解了ggplot中的映射,对这个图的修改和美化其实就是修改geom_roc()函数里面的参数,以及用其他ggplot元素进行优化。

默认曲线上会显示阈值cutoff的数值,我们可以关闭它:

ggplot(test, aes(d = D, m = M1)) + geom_roc(n.cuts = 0)
1-3.png

修改它:

ggplot(test, aes(d = D, m = M1)) + geom_roc(n.cuts = 5, labelsize = 5, labelround = 2)

1-4.png

使用plotROC提供的风格:

styledplot <- basicplot + style_roc() styledplot

1-5.png

将标签加在曲线上:

direct_label(basicplot, labels = "Biomarker", nudge_y = -.1) + style_roc()

1-6.png

绘制多条曲线

plotROC提供的函数melt_roc()可以将多个变量列变为长格式,方便数据的绘制:

1-12.png

画比较图:

ggplot(longtest, aes(d = D, m = M, color = name)) + geom_roc() + style_roc()

1-7.png

还有其他一些功能,请查看文档(http://sachsmc.github.io/plotROC/)学习,这里最后介绍一下我封装的一个函数,便于两组ROC比较的使用,感兴趣的朋友可以自定义再修改和优化。

plotROC <- function(.data, predict_col, target, group, positive=1, all=TRUE){    
if(!(require(tidyverse) & require(plotROC))){         stop("--> tidyverse and plotROC packages are required..")    }     
predict_col <- enquo(predict_col)    
target <- enquo(target)     group  <- enquo(group)    
predictN <- quo_name(predict_col)    
groupN   <- quo_name(group)    
df <- .data %>% dplyr::select(!! predict_col, !! target, !! group) %>%        
mutate(targetN = ifelse(!! target == positive, 1, 0)) %>% as.data.frame()    
if (all){        
df2 <- df         
df2[, groupN] <- "ALL"        
df <- rbind(df, df2)   
 }    
 p  <- df %>%  ggplot(aes_string(m = predictN,                                     d = "targetN",                                    color = groupN)) + geom_roc(show.legend = TRUE, labels=FALSE)    
 p <- p + ggpubr::theme_classic2()   
  ng <- levels(factor(df[, groupN]))   
   if(length(ng) == 3){        
   auc <- calc_auc(p)$AUC        
   names(auc) <- ng        
   auc <- base::sort(auc, decreasing = TRUE)       
   p <- p + annotate("text", x = .75, y = .25,                           label = paste(names(auc[1], " AUC =", round(auc[1], 3), "\n",                                        
   names(auc[2], " AUC =", round(auc[2], 3), "\n",                                        
   names(auc[3], " AUC =", round(auc[3], 3), "\n"),                          size = 6)   
    }   
     p + xlab("1 - Specificity") + ylab("Sensitivity") +  
            scale_x_continuous(expand = c(0, 0)) + scale_y_continuous(expand = c(0, 0))
            }

使用:

plotROC(longtest, predict_col = M, target = D, group = name, positive = 1)

参数1:提供数据框

参数2:提供预测数值列

参数3:提供二分类信息列(尽量为0-1,字符也可以)

参数4:提供一个组别

参数5:这里1表示成功,如果target是success和failure,可以知道positive="success"

注意,这里只有3条曲线绘制时才会给出AUC在图上,可以修改函数进行自定义

1-8.png

默认会画出两组及其融合的曲线,可以添加选项all=FALSE去掉ALL曲线。

有读者谈到如何修改,之前之所以没写多条曲线添加AUC,是因为涉及一些文本图像的微调,实际使用时需要自定义一下如果想要添加6条曲线,在加上ALL,就是7条,请补充函数中的if代码块

if(length(ng) == 3){        
auc <- calc_auc(p)$AUC        
names(auc) <- ng       
 auc <- base::sort(auc, decreasing = TRUE)        
 p <- p + annotate("text", x = .75, y = .25,                           label = paste(names(auc[1], " AUC =", round(auc[1], 3), "\n",                                       
  names(auc[2], " AUC =", round(auc[2], 3), "\n",                                        
  names(auc[3], " AUC =", round(auc[3], 3), "\n"),                          size = 6)    
  }

曲线太多时可能文本注释添加需要注意下位置和大小。注意上述更改未测试,请根据实际情况调整。

pROC

pROC是一个相对plotROC更强大的R包,不同于plotROC基于ggplot2的创建,pROC自身构建了比较完整的ROC分析和绘图体系。

该包发表文章为:

Xavier Robin, Natacha Turck, Alexandre Hainard, Natalia Tiberti, Frédérique Lisacek, Jean-Charles Sanchez and Markus Müller (2011). pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics, 12, p. 77. DOI: 10.1186/1471-2105-12-77.

目前谷歌搜索已经有超过2000次引用。

1-9.png

该包创建的图像似乎更加圆润。

安装

安装

install.packages("pROC")

导入

library(pROC)

获取帮助

?pROC

使用

不过相对于plotROC,它的图形绘制更为复杂(样例代码参见https://web.expasy.org/pROC/screenshots.html)。比如
1-10.png

其代码为:

library(pROC)
data(aSAH)
plot.roc(aSAH$outcome, aSAH$s100b, # data
percent=TRUE, # show all values in percent
partial.auc=c(100, 90), partial.auc.correct=TRUE, # define a partial AUC (pAUC)
print.auc=TRUE, #display pAUC value on the plot with following options:
print.auc.pattern="Corrected pAUC (100-90%% SP):\n%.1f%%", print.auc.col="#1c61b6",
auc.polygon=TRUE, auc.polygon.col="#1c61b6", # show pAUC as a polygon
max.auc.polygon=TRUE, max.auc.polygon.col="#1c61b622", # also show the 100% polygon
main="Partial AUC (pAUC)")
plot.roc(aSAH$outcome, aSAH$s100b,
percent=TRUE, add=TRUE, type="n", # add to plot, but don't re-add the ROC itself (useless)
partial.auc=c(100, 90), partial.auc.correct=TRUE,
partial.auc.focus="se", # focus pAUC on the sensitivity
print.auc=TRUE, print.auc.pattern="Corrected pAUC (100-90%% SE):\n%.1f%%", print.auc.col="#008600",
print.auc.y=40, # do not print auc over the previous one
auc.polygon=TRUE, auc.polygon.col="#008600",
max.auc.polygon=TRUE, max.auc.polygon.col="#00860022")

不过我们常用的一般是

1-11.png

这样的图形,我们参考代码修改和自定义即可:

library(pROC)
data(aSAH)
rocobj1 <- plot.roc(aSAH$outcome, aSAH$s100,
main="Statistical comparison", percent=TRUE, col="#1c61b6")
rocobj2 <- lines.roc(aSAH$outcome, aSAH$ndka, percent=TRUE, col="#008600")
testobj <- roc.test(rocobj1, rocobj2)
text(50, 50, labels=paste("p-value =", format.pval(testobj$p.value)), adj=c(0, .5))legend("bottomright", 
legend=c("S100B", "NDKA"), col=c("#1c61b6", "#008600"), lwd=2)

这个图显示了pROC包最重要几个函数的使用,第一个是plot.roc(),它可以绘制ROC曲线,并返回一个ROC对象,里面包含该曲线的众多有用信息,并为后续的分析做基础,lines.roc()为当前ROC曲线上增添新的ROC曲线。不仅如此,roc.test()函数提供了对曲线进行检验,检验的方法分为3种,可以自己选择,有兴趣的朋友不妨再深入看看。

▍本文版权(包括图片及文字)属于“优雅R”(微信公众号:elegant-r),禁止二次转载,如需转载请联系w_shixiang@163.com

上一篇下一篇