在【r<-繪圖|ROC】ROC的計(jì)算與繪制這篇文章中我講了ROC曲線的本質(zhì)以及如何計(jì)算和繪制ROC曲線。注意，我這里談到的ROC并未曾涉及機(jī)器學(xué)習(xí)模型的擬合與預(yù)測(cè)，而是指存在一組真實(shí)的連續(xù)型數(shù)值數(shù)據(jù)設(shè)定閾值的不同對(duì)響應(yīng)變量（二分類）的影響（真陽(yáng)性率、假陽(yáng)性率）。

這一篇文章我們學(xué)習(xí)兩個(gè)跟ROC相關(guān)的R包：

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

plotROC

plotROC包較為簡(jiǎn)單與單一，它就是用來(lái)繪制ROC曲線的，包中定義的函數(shù)基于ggplot2，因此我們可以結(jié)合ggplot2使用和修改、美化圖形結(jié)果。

# 從GitHub上安裝
devtools::install_github("hadley/ggplot2")
devtools::install_github("sachsmc/plotROC")
library(plotROC)
# 從CRAN
install.packages("plotROC")

快速使用

plotROC提供了Shiny應(yīng)用，只需要鍵入

shiny_plotROC()

即可通過(guò)圖形界面使用。

咱們還是來(lái)看命令吧，要有點(diǎn)難度不是？

命令行使用

導(dǎo)入包與創(chuàng)建模擬數(shù)據(jù)：

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)

簡(jiǎn)單繪圖：

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

這里我們唯一需要理清的是d與m映射是什么，現(xiàn)在我們查看下生成的數(shù)據(jù)框：

上述畫圖只使用到了D與M1，只關(guān)注這兩列即可。D是一個(gè)0-1列，即表示結(jié)果的兩分類信息，M1是一個(gè)數(shù)值型數(shù)據(jù)。我們可以姑且稱d為decision縮寫，m為measurement縮寫。

一旦我們理解了ggplot中的映射，對(duì)這個(gè)圖的修改和美化其實(shí)就是修改geom_roc()函數(shù)里面的參數(shù)，以及用其他ggplot元素進(jìn)行優(yōu)化。

默認(rèn)曲線上會(huì)顯示閾值cutoff的數(shù)值，我們可以關(guān)閉它：

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提供的風(fēng)格：

styledplot <- basicplot + style_roc()
styledplot

將標(biāo)簽加在曲線上：

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

繪制多條曲線

plotROC提供的函數(shù)melt_roc()可以將多個(gè)變量列變?yōu)殚L(zhǎng)格式，方便數(shù)據(jù)的繪制：

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

還有其他一些功能，請(qǐng)查看文檔（http://sachsmc.github.io/plotROC/）學(xué)習(xí)，這里最后介紹一下我封裝的一個(gè)函數(shù)，便于兩組ROC比較的使用，感興趣的朋友可以自定義再修改和優(yōu)化。

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)
# 參數(shù)1：提供數(shù)據(jù)框
# 參數(shù)2：提供預(yù)測(cè)數(shù)值列
# 參數(shù)3：提供二分類信息列（盡量為0-1，字符也可以）
# 參數(shù)4：提供一個(gè)組別
# 參數(shù)5：這里1表示成功，如果target是success和failure，可以知道positive="success"
# 注意，這里只有3條曲線繪制時(shí)才會(huì)給出AUC在圖上，可以修改函數(shù)進(jìn)行自定義

默認(rèn)會(huì)畫出兩組及其融合的曲線，可以添加選項(xiàng)all=FALSE去掉ALL曲線。

有讀者談到如何修改，之前之所以沒(méi)寫多條曲線添加AUC，是因?yàn)樯婕耙恍┪谋緢D像的微調(diào)，實(shí)際使用時(shí)需要自定義一下
如果想要添加6條曲線，在加上ALL，就是7條，請(qǐng)補(bǔ)充函數(shù)中的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)
    }

曲線太多時(shí)可能文本注釋添加需要注意下位置和大小。注意上述更改未測(cè)試，請(qǐng)根據(jù)實(shí)際情況調(diào)整。

pROC

pROC是一個(gè)相對(duì)plotROC更強(qiáng)大的R包，不同于plotROC基于ggplot2的創(chuàng)建，pROC自身構(gòu)建了比較完整的ROC分析和繪圖體系。

該包發(fā)表文章為：

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.

目前谷歌搜索已經(jīng)有超過(guò)2000次引用。

pROC繪圖

該包創(chuàng)建的圖像似乎更加圓潤(rùn)。

安裝

# 安裝
install.packages("pROC")
# 導(dǎo)入
library(pROC)
# 獲取幫助
?pROC

使用

不過(guò)相對(duì)于plotROC，它的圖形繪制更為復(fù)雜（樣例代碼參見(jiàn)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")

不過(guò)我們常用的一般是

這樣的圖形，我們參考代碼修改和自定義即可：

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)

這個(gè)圖顯示了pROC包最重要幾個(gè)函數(shù)的使用，第一個(gè)是plot.roc()，它可以繪制ROC曲線，并返回一個(gè)ROC對(duì)象，里面包含該曲線的眾多有用信息，并為后續(xù)的分析做基礎(chǔ)，lines.roc()為當(dāng)前ROC曲線上增添新的ROC曲線。不僅如此，roc.test()函數(shù)提供了對(duì)曲線進(jìn)行檢驗(yàn)，檢驗(yàn)的方法分為3種，可以自己選擇，有興趣的朋友不妨再深入看看。