之前我寫過一篇使用pROC包畫好看的ROC曲線的教程,那篇教程中使用的是pROC撩轰,這個(gè)包可以快速擬合ROC曲線胯甩,然而這個(gè)包需要提前進(jìn)行運(yùn)算結(jié)果堪嫂,并且不能直接顯示AUC值等偎箫,今天推薦一個(gè)另一個(gè)繪制ROC的包multipleROC皆串,顧名思義淹办,這個(gè)包是可以一次性繪制多條ROC曲線的恶复,并且也是基于ggplot2怜森。
目前這個(gè)包作者沒有上傳CRAN或BiocManager谤牡,只能通過Github安裝副硅,地址為https://github.com/cardiomoon/multipleROC
安裝multipleROC
remotes::install_github("cardiomoon/multipleROC")
如果無法訪問GitHub翅萤,也可以導(dǎo)入到Gitee后進(jìn)行安裝
remotes::install_git("https://gitee.com/swcyo/multipleROC/")
數(shù)據(jù)演示
我們使用仙桃學(xué)術(shù)上的一個(gè)診斷性ROC示例數(shù)據(jù)為例進(jìn)行演示(下載請(qǐng)點(diǎn)擊xlxs鏈接)恐疲。
library(readxl)
ROC <- read_excel("~/Desktop/ROC曲線.xlsx")
探索性分析
我們可以事先看一下group1和group2兩組在a變量中的差別套么,使用webr包培己,先看看結(jié)果如何
library(webr)
library(ggplot2)
library(dplyr)
library(tidyr)
ROC %>%
group_by(outcome) %>%
numSummaryTable(a)
outcome | n | mean | sd | median | trimmed | mad | min | max | range |
---|---|---|---|---|---|---|---|---|---|
group1 | 40.00 | 1.51 | 0.55 | 1.45 | 1.51 | 0.63 | 0.59 | 2.44 | 1.86 |
group2 | 32.00 | 1.00 | 0.55 | 1.00 | 0.99 | 0.63 | 0.13 | 1.99 | 1.86 |
也可以使用箱示圖和密度圖進(jìn)行展示胚泌,見Figure 1所示省咨。
p1<- ggplot(data=ROC)+geom_density(aes(x=a,fill=outcome),alpha=0.5)
p2<-ggplot(data=ROC)+geom_boxplot(aes(x=outcome,y=a,fill=outcome),alpha=0.5)
cowplot::plot_grid(p1,p2)
同法可以顯示b和c變量的結(jié)果诸迟,我們暫時(shí)以boxplot展示
p3<-ggplot(data=ROC)+geom_boxplot(aes(x=outcome,y=b,fill=outcome),alpha=0.5)
p4<-ggplot(data=ROC)+geom_boxplot(aes(x=outcome,y=c,fill=outcome),alpha=0.5)
cowplot::plot_grid(p2,p3,p4,labels = "AUTO",nrow = 1)
-- 雖然探索性分析可以判斷兩組的差異茸炒,但是無法確定最佳截?cái)嘀嫡笪矡o妨評(píng)估預(yù)測效能壁公。
ROC曲線的繪制
繪制ROC曲線是確定最佳截?cái)嘀档挠杏梅椒ㄖ簧鹣睢D梢允褂靡韵翿命令執(zhí)行ROC分析紊册。下面的R命令使一個(gè)類為multipleROC的對(duì)象,并進(jìn)行繪圖囊陡。
由于默認(rèn)的函數(shù)中分組需為0和1芳绩,因此需要將group1和group2進(jìn)行賦值撞反,我們將group1定義為0妥色,group2定義為1,我們繪制a變量在兩組中的ROC圖片嘹害,我們可以使用multipleROC()
語句一步計(jì)算,可以看到最佳截?cái)嘀邓北悖珹UC值笔呀,另外敏感度、特異度都是可以直接顯示的髓需,見Figure 3所示。僚匆。
ROC$group<-ifelse(ROC$outcome=='group1',0,1) # 將group1定義為0微渠,否則為1
library(multipleROC)
a=multipleROC(group~a,data=ROC)
如果不想顯示那么多結(jié)果的話,也可以plot_ROC()
函數(shù)一個(gè)個(gè)設(shè)置是否顯示
plot_ROC(a,
show.points = TRUE,
show.eta = TRUE,
show.sens = TRUE,
show.AUC = TRUE,
facet = FALSE )
AUC和p值
在Figure 3的右下角敛助,您可以看到曲線下面積(AUC)和Wilcoxon秩和檢驗(yàn)的p值。p值來自以下計(jì)算結(jié)果屋确。
wilcox.test(ROC$a,ROC$group)
##
## Wilcoxon rank sum test with continuity correction
##
## data: ROC$a and ROC$group
## W = 4416, p-value = 1.294e-13
## alternative hypothesis: true location shift is not equal to 0
AUC值則通過multipleROC包的simpleAUC()
函數(shù)進(jìn)行運(yùn)算纳击,函數(shù)如下:
simpleAUC <- function(df){
df=df[order(df$x,decreasing=TRUE),]
TPR=df$sens
FPR=df$fpr
dFPR <- c(diff(FPR), 0)
dTPR <- c(diff(TPR), 0)
sum(TPR * dFPR) + sum(dTPR * dFPR)/2
}
那么,我們直接直接只有simpleAUC(aauc直接看到完整的AUC值
simpleAUC(a$df) ## 函數(shù)法
## [1] 0.7328125
a$auc # 直接提取法
## [1] 0.7328125
同樣的,我們直接提取截?cái)帱c(diǎn)(cutpoint)和最佳截?cái)嘀担∣ptimal Cutoff value)
a$cutpoint
## [1] 0.5136663
a$cutoff
## a ## 54 1.082828
將結(jié)果轉(zhuǎn)換為pROC對(duì)象
如果你更習(xí)慣pROC的結(jié)果刨啸,使用multipleROC2roc()
函數(shù)堡赔,可以直接將結(jié)果轉(zhuǎn)換為 pROC的roc 對(duì)象
a2<-multipleROC2roc(a)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
class(a) ##a的類型為multipleROC
## [1] "multipleROC"
class(a2) ##a2已經(jīng)轉(zhuǎn)換為roc的類型了
## [1] "roc"
pROC::auc(a2) ## 我們用pROC看auc的結(jié)果
## Area under the curve: 0.7328
我們可以使用pROC的繪圖函數(shù)對(duì)a2進(jìn)行繪圖了设联,我們比較以下兩種結(jié)果吧善已,見Figure 4所示。
library(pROC)
p5<-ggroc(a2, legacy.axes = TRUE)+
geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1),
color="darkgrey",linetype=4)+
theme_bw()+ggtitle("pROC")
p6<-plot(a)+ggtitle("multipleROC")
cowplot::plot_grid(p5,p6)
多個(gè)ROC曲線的繪制
可以用多個(gè)函數(shù)進(jìn)行多個(gè)ROC的曲線换团,可以使用plot_ROC(list())
一個(gè)個(gè)繪制曲線,見Figure 5所示宫蛆。
a=multipleROC(group~a,data=ROC,plot=FALSE)
b=multipleROC(group~b,data=ROC,plot=FALSE)
c=multipleROC(group~c,data=ROC,plot=FALSE)
plot_ROC(list(a,b,c),
show.eta=FALSE,#不顯示截點(diǎn)
show.sens=FALSE #不顯示各種率
)
當(dāng)然艘包,如果你不想寫那么多代碼的話,也可以直接使用plot_ROC2()
函數(shù)直接繪制想虎,是不是簡單的多卦尊。
plot_ROC2(yvar="group",xvars=c("a","b","c"),dataname="ROC")
分面顯示
將三張圖放在一起舌厨,可以看到數(shù)值重疊影響顏值岂却,我們可以用ggplot2的facet進(jìn)行分面繪制邓线。
plot_ROC(list(a,b,c),facet=TRUE)
可以發(fā)現(xiàn)分面的標(biāo)簽?zāi)J(rèn)是1淌友,2骇陈,3,我們可以使用Y叔團(tuán)隊(duì)開發(fā)的ggfun這個(gè)包的facet_set()
函數(shù)進(jìn)行快速的修改
library(ggfun) ## 必須要先安裝好
plot_ROC(list(a,b,c),facet=TRUE)+
facet_set(label=c(`1`="a", `2`="b", `3`="c"))
換一種分面顯示
使用ggplot2包的facet_grid
可以換一個(gè)分面顯示方式
plot_ROC(list(a,b,c))+facet_grid(no~.)+
facet_set(label=c(`1`="a", `2`="b", `3`="c"))
由于是基于ggplot2語句你雌,所以我們可以使用ggtitle
添加標(biāo)題,還可以更換主題等等二汛。婿崭。。