最近工作需要繪制ROC曲線，對(duì)該曲線的計(jì)算細(xì)節(jié)進(jìn)行了一番摸索。當(dāng)前搜索ROC曲線一般跟機(jī)器學(xué)習(xí)相關(guān)聯(lián)，導(dǎo)致我對(duì)它的概念有了曲解，理所當(dāng)然地以為它只是一個(gè)用于機(jī)器學(xué)習(xí)的分類(lèi)器評(píng)估標(biāo)準(zhǔn)，所以在繪制曲線前使用邏輯回歸（我的響應(yīng)變量是0-1類(lèi)型）對(duì)數(shù)據(jù)建模分析。實(shí)則不然，ROC曲線適用于任何判斷0-1類(lèi)型（真假、成功失敗等二分類(lèi)）響應(yīng)結(jié)果閾值分割效果的評(píng)估。
關(guān)于ROC曲線的概念，可以參考ROC曲線與AUC值這篇博文學(xué)習(xí)和理解，如果你能看懂下面這條曲線，基本上你就理解了ROC存在的意義。
對(duì)于0，1兩類(lèi)分類(lèi)問(wèn)題，一些分類(lèi)器得到的結(jié)果往往不是0，1這樣的標(biāo)簽，如神經(jīng)網(wǎng)絡(luò)得到諸如0.5，0.8這樣的分類(lèi)結(jié)果。這時(shí)，我們?nèi)藶槿∫粋€(gè)閾值，比如0.4，那么小于0.4的歸為0類(lèi)，大于等于0.4的歸為1類(lèi)，可以得到一個(gè)分類(lèi)結(jié)果。同樣，這個(gè)閾值我們可以取0.1或0.2等等。取不同的閾值，最后得到的分類(lèi)情況也就不同。如下面這幅圖：
藍(lán)色表示原始為負(fù)類(lèi)分類(lèi)得到的統(tǒng)計(jì)圖，紅色表示原始為正類(lèi)得到的統(tǒng)計(jì)圖。那么我們?nèi)∫粭l直線，直線左邊分為負(fù)類(lèi)，直線右邊分為正類(lèi)，這條直線也就是我們所取的閾值。閾值不同，可以得到不同的結(jié)果，但是由分類(lèi)器決定的統(tǒng)計(jì)圖始終是不變的。這時(shí)候就需要一個(gè)獨(dú)立于閾值，只與分類(lèi)器有關(guān)的評(píng)價(jià)指標(biāo)，來(lái)衡量特定分類(lèi)器的好壞。還有在類(lèi)不平衡的情況下，如正樣本有90個(gè)，負(fù)樣本有10個(gè)，直接把所有樣本分類(lèi)為正樣本，得到識(shí)別率為90%，但這顯然是沒(méi)有意義的。如上就是ROC曲線的動(dòng)機(jī)。
-- ROC曲線與AUC值
在R里面，有ROCR與專(zhuān)門(mén)的機(jī)器學(xué)習(xí)包mlr可以進(jìn)行建模和繪制ROC曲線，以及相關(guān)參量的計(jì)算。實(shí)際上，不需要使用任何模型，也可以繪制ROC曲線，因?yàn)镽OC曲線的繪制就是選擇閾值與計(jì)算當(dāng)前閾值下假陽(yáng)性率與真陽(yáng)性率變化的過(guò)程。
上述提到的兩個(gè)包使用有些復(fù)雜，實(shí)際上我要用的也不是它們，關(guān)于ROC的計(jì)算，仔細(xì)思考寫(xiě)個(gè)程序就能搞定。核心是計(jì)算假陽(yáng)性、真陽(yáng)性率，首先要計(jì)算下方混淆矩陣中的各個(gè)參數(shù)。
calcROC <- function(.data, predict_var, target, group_var, positive="success"){
# predic_var here must be a numeric value
require(tidyverse)
predict_var <- enquo(predict_var)
target <- enquo(target)
group_var <- enquo(group_var)
groups <- .data %>% filter(!is.na(!! predict_var)) %>% select(!! group_var) %>%
unlist() %>% table() %>% names()
total_res <- list()
# process groups one by one
j <- 1
for (i in groups){
df <- list()
df <- .data %>% filter(!is.na(!! predict_var), !! group_var == i) %>%
arrange(desc(!! predict_var)) %>%
mutate(isPositive = ifelse(!! target == positive, 1, 0))
# select a threshold, calculate true positive and false positive value
ths <- df %>% select(!! predict_var) %>% unlist
mat <- base::sapply(ths, function(th){
# true positive
tp <- df %>% filter(!! predict_var >= th) %>% filter(isPositive == 1) %>% nrow
# false positive
fp <- df %>% filter(!! predict_var >= th) %>% filter(isPositive == 0) %>% nrow
# true negative
tn <- df %>% filter(!! predict_var < th) %>% filter(isPositive == 0) %>% nrow
# false negative
fn <- df %>% filter(!! predict_var < th) %>% filter(isPositive == 1) %>% nrow
# true positive rate
tpr <- tp / (tp + fn)
# false positive rate
fpr <- fp / (fp + tn)
return(c(tp, fp, tn, fn, tpr, fpr))
})
res <- t(mat)
res <- data.frame(res)
# fake a (0, 0) point
res <- rbind(c(rep(NA, 4), 0, 0), res)
colnames(res) <- c("tp", "fp", "tn", "fn", "tpr", "fpr")
res$Group <- i
total_res[[j]] <- res
j <- j + 1
}
dat <- base::Reduce(rbind, total_res)
return(dat)
}
上面函數(shù)的運(yùn)行需要保證tidyverse包已經(jīng)安裝，寫(xiě)法遵從tidyverse語(yǔ)法，涉及不少管道操作，如果你只想使用，直接拷貝運(yùn)行即可，如果想要理解過(guò)程，需要dplyr使用和編程（列舉一篇筆記）的一些知識(shí)。
略微講解一下使用：
> args(calcROC)
function (.data, predict_var, target, group_var, positive = "success")
函數(shù)有5個(gè)參數(shù)，第一個(gè)是包含數(shù)據(jù)的數(shù)據(jù)框；第二個(gè)是預(yù)測(cè)變量，一個(gè)數(shù)值向量；第三個(gè)是目標(biāo)變量，包含0-1信息（成功或失敗，等等）；第四個(gè)是一個(gè)分組參數(shù)，一般我們會(huì)比較兩組或多組ROC曲線的差異；第五個(gè)是給出成功（或1）是用什么指定的，比如目標(biāo)變量中success指代成功。
結(jié)果會(huì)返回一個(gè)數(shù)據(jù)框：
calcROC(dat, mut, isBenefit, Gender)
tp fp tn fn tpr fpr Group
1 NA NA NA NA 0.0000000 0.0 Female
TMB_NonsynSNP1 1 0 10 7 0.1250000 0.0 Female
TMB_NonsynSNP2 2 0 10 6 0.2500000 0.0 Female
TMB_NonsynSNP3 3 0 10 5 0.3750000 0.0 Female
TMB_NonsynSNP4 4 0 10 4 0.5000000 0.0 Female
TMB_NonsynSNP5 4 1 9 4 0.5000000 0.1 Female
TMB_NonsynSNP6 5 1 9 3 0.6250000 0.1 Female
TMB_NonsynSNP7 6 1 9 2 0.7500000 0.1 Female
TMB_NonsynSNP8 6 2 8 2 0.7500000 0.2 Female
TMB_NonsynSNP9 7 2 8 1 0.8750000 0.2 Female
TMB_NonsynSNP10 7 3 7 1 0.8750000 0.3 Female
TMB_NonsynSNP11 8 3 7 0 1.0000000 0.3 Female
TMB_NonsynSNP12 8 4 6 0 1.0000000 0.4 Female
TMB_NonsynSNP13 8 5 5 0 1.0000000 0.5 Female
TMB_NonsynSNP14 8 6 4 0 1.0000000 0.6 Female
TMB_NonsynSNP15 8 7 3 0 1.0000000 0.7 Female
TMB_NonsynSNP16 8 8 2 0 1.0000000 0.8 Female
TMB_NonsynSNP17 8 9 1 0 1.0000000 0.9 Female
TMB_NonsynSNP18 8 10 0 0 1.0000000 1.0 Female
11 NA NA NA NA 0.0000000 0.0 Male
TMB_NonsynSNP19 0 1 9 6 0.0000000 0.1 Male
TMB_NonsynSNP21 1 1 9 5 0.1666667 0.1 Male
TMB_NonsynSNP31 1 2 8 5 0.1666667 0.2 Male
TMB_NonsynSNP41 2 2 8 4 0.3333333 0.2 Male
TMB_NonsynSNP51 2 3 7 4 0.3333333 0.3 Male
TMB_NonsynSNP61 3 3 7 3 0.5000000 0.3 Male
TMB_NonsynSNP71 4 3 7 2 0.6666667 0.3 Male
TMB_NonsynSNP81 4 4 6 2 0.6666667 0.4 Male
TMB_NonsynSNP91 5 4 6 1 0.8333333 0.4 Male
TMB_NonsynSNP101 5 5 5 1 0.8333333 0.5 Male
TMB_NonsynSNP111 5 6 4 1 0.8333333 0.6 Male
TMB_NonsynSNP121 5 7 3 1 0.8333333 0.7 Male
TMB_NonsynSNP131 5 8 2 1 0.8333333 0.8 Male
TMB_NonsynSNP141 5 9 1 1 0.8333333 0.9 Male
TMB_NonsynSNP151 6 9 1 0 1.0000000 0.9 Male
TMB_NonsynSNP161 6 10 0 0 1.0000000 1.0 Male
注意函數(shù)調(diào)用時(shí)寫(xiě)法符合dplyr，相關(guān)參數(shù)可沒(méi)有打引號(hào)，注意下~
左邊第一列不用管，是一個(gè)無(wú)意義的行名，結(jié)果共有7列，而畫(huà)圖只需要最后的三列。每個(gè)組都有一行是偽造的，用于填充圖像中的（0,0）點(diǎn)，也即閾值取無(wú)窮大時(shí)。
調(diào)用ggpubr包畫(huà)圖：
ggpubr::ggline(data = calcROC(dat, mut, isBenefit, Gender), x = "fpr", y = "tpr", linetype = "Group", shape = "Group")
如果想實(shí)際畫(huà)出美麗的曲線，請(qǐng)參考新文章【r<-ROC|包】分析與可視化ROC——plotROC、pROC
使用 pROC 包的數(shù)據(jù)更新一個(gè)示例
library(pROC) data(aSAH) head(aSAH) calcROC(aSAH, s100b, outcome, gender, positive = "Good")
