利用多重 ROC 快速绘制 ROC 曲线
之前我写过一篇使用pROC包画好看的ROC曲线的教程，那篇教程中使用的是pROC，这个包可以快速拟合ROC曲线，然而这个包需要提前进行运算结果，并且不能直接显示AUC值等，今天推荐一个另一个绘制ROC的包multipleROC，顾名思义，这个包是可以一次性绘制多条ROC曲线的，并且也是基于ggplot2。
目前这个包作者没有上传CRAN或BiocManager，只能通过Github安装，地址为https://github.com/cardiomoon/multipleROC
安装multipleROC
remotes::install_github("cardiomoon/multipleROC")
如果无法访问GitHub，也可以导入到Gitee后进行安装
remotes::install_git("https://gitee.com/swcyo/multipleROC/")
数据演示
我们使用仙桃学术上的一个诊断性ROC示例数据为例进行演示（下载请点击xlxs链接）。
library(readxl)
ROC <- read_excel("~/Desktop/ROC曲线.xlsx")
探索性分析
我们可以事先看一下group1和group2两组在a变量中的差别，使用webr包，先看看结果如何
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 |
也可以使用箱示图和密度图进行展示，见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变量的结果，我们暂时以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)
-- 虽然探索性分析可以判断两组的差异，但是无法确定最佳截断值，也无妨评估预测效能。
ROC曲线的绘制
绘制ROC曲线是确定最佳截断值的有用方法之一。您可以使用以下R命令执行ROC分析。下面的R命令使一个类为multipleROC的对象，并进行绘图。
由于默认的函数中分组需为0和1，因此需要将group1和group2进行赋值，我们将group1定义为0，group2定义为1，我们绘制a变量在两组中的ROC图片，我们可以使用multipleROC()语句一步计算，可以看到最佳截断值，AUC值，另外敏感度、特异度都是可以直接显示的，见Figure 3所示。。
ROC$group<-ifelse(ROC$outcome=='group1',0,1) # 将group1定义为0，否则为1
library(multipleROC)
a=multipleROC(group~a,data=ROC)
如果不想显示那么多结果的话，也可以plot_ROC()函数一个个设置是否显示
plot_ROC(a,
show.points = TRUE,
show.eta = TRUE,
show.sens = TRUE,
show.AUC = TRUE,
facet = FALSE )
AUC和p值
在Figure 3的右下角，您可以看到曲线下面积（AUC）和Wilcoxon秩和检验的p值。p值来自以下计算结果。
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()函数进行运算，函数如下:
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(a
auc直接看到完整的AUC值
simpleAUC(a$df) ## 函数法
## [1] 0.7328125
a$auc # 直接提取法
## [1] 0.7328125
同样的，我们直接提取截断点(cutpoint)和最佳截断值（Optimal Cutoff value）
a$cutpoint
## [1] 0.5136663
a$cutoff
## a ## 54 1.082828
将结果转换为pROC对象
如果你更习惯pROC的结果，使用multipleROC2roc()函数，可以直接将结果转换为 pROC的roc 对象
a2<-multipleROC2roc(a)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
class(a) ##a的类型为multipleROC
## [1] "multipleROC"
class(a2) ##a2已经转换为roc的类型了
## [1] "roc"
pROC::auc(a2) ## 我们用pROC看auc的结果
## Area under the curve: 0.7328
我们可以使用pROC的绘图函数对a2进行绘图了，我们比较以下两种结果吧，见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)
多个ROC曲线的绘制
可以用多个函数进行多个ROC的曲线，可以使用plot_ROC(list())一个个绘制曲线，见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,#不显示截点
show.sens=FALSE #不显示各种率
)
当然，如果你不想写那么多代码的话，也可以直接使用plot_ROC2()函数直接绘制，是不是简单的多。
plot_ROC2(yvar="group",xvars=c("a","b","c"),dataname="ROC")
分面显示
将三张图放在一起，可以看到数值重叠影响颜值，我们可以用ggplot2的facet进行分面绘制。
plot_ROC(list(a,b,c),facet=TRUE)
可以发现分面的标签默认是1，2，3，我们可以使用Y叔团队开发的ggfun这个包的facet_set()函数进行快速的修改
library(ggfun) ## 必须要先安装好
plot_ROC(list(a,b,c),facet=TRUE)+
facet_set(label=c(`1`="a", `2`="b", `3`="c"))
换一种分面显示
使用ggplot2包的facet_grid可以换一个分面显示方式
plot_ROC(list(a,b,c))+facet_grid(no~.)+
facet_set(label=c(`1`="a", `2`="b", `3`="c"))
由于是基于ggplot2语句，所以我们可以使用ggtitle添加标题，还可以更换主题等等。。。