【r<-ROC|包】分析与可视化ROC——plotROC、pROC

在【r<-绘图|ROC】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

这里我们唯一需要理清的是d与m映射是什么，现在我们查看下生成的数据框：

上述画图只使用到了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)

修改它：

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

使用plotROC提供的风格：

styledplot <- basicplot + style_roc()
styledplot

将标签加在曲线上：

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

绘制多条曲线

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

longtest <- melt_roc(test, "D", c("M1", "M2"))
head(longtest)

##     D          M name
## M11 1 1.48117155   M1
## M12 1 0.61994478   M1
## M13 0 0.57613345   M1
## M14 1 0.85433197   M1
## M15 0 0.05258342   M1
## M16 1 0.66703989   M1

画比较图：

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

还有其他一些功能，请查看文档（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在图上，可以修改函数进行自定义

默认会画出两组及其融合的曲线，可以添加选项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)
    }

为：

    if(length(ng) == 7){
        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",
                                        names(auc)[4], " AUC =", round(auc[4], 3), "\n"
                                        names(auc)[5], " AUC =", round(auc[5], 3), "\n"
                                        names(auc)[6], " AUC =", round(auc[6], 3), "\n"
                                        names(auc)[7], " AUC =", round(auc[7], 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次引用。

pROC绘图

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

安装

# 安装
install.packages("pROC")
# 导入
library(pROC)
# 获取帮助
?pROC

使用

不过相对于plotROC，它的图形绘制更为复杂（样例代码参见https://web.expasy.org/pROC/screenshots.html）。

比如

其代码为：

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")

不过我们常用的一般是

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

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种，可以自己选择，有兴趣的朋友不妨再深入看看。

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 206,839评论 6赞 482
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 88,543评论 2赞 382
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 153,116评论 0赞 344
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 55,371评论 1赞 279
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 64,384评论 5赞 374
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 49,111评论 1赞 285
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 38,416评论 3赞 400
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 37,053评论 0赞 259
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 43,558评论 1赞 300
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,007评论 2赞 325
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,117评论 1赞 334
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,756评论 4赞 324
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 39,324评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,315评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,539评论 1赞 262
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 45,578评论 2赞 355
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,877评论 2赞 345

推荐阅读更多精彩内容

  • 面试题
    Swift1> Swift和OC的区别1.1> Swift没有地址/指针的概念1.2> 泛型1.3> 类型严谨 对...
    cosWriter阅读 11,089评论 1赞 32
  • 如梦的生活
    我站在梦的两边. 敲打着生活. 酸菜瘦肉炒饭. 烫青菜. 成了我的主题歌. 一场不大不小的雨. 凝固了整个空气. ...
    几分秋意浓阅读 328评论 2赞 0
  • 2017.6.15 星期四
    第四十天 四十天啦，又多了十天。撒花。 偷拍你的照片一定把你拍胖了，恩，一定是这样。手臂没那么粗的是吧，哈哈，星期...
    漫步的野马阅读 258评论 0赞 0
  • 且歌且行
    且歌且行
    梦之鸭阅读 242评论 0赞 0