使用multipleROC快速绘制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)

也可以使用箱示图和密度图进行展示，见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)

Figure 1: group1和group2两组在a变量中的差别

同法可以显示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)

Figure 2: group1和group2两组在三变量中的差别

-- 虽然探索性分析可以判断两组的差异，但是无法确定最佳截断值，也无妨评估预测效能。

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)

Figure 3: a变量在两组的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(aauc直接看到完整的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)

Figure 4: pROC和multipleROC的结果对比

多个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 #不显示各种率
         )

Figure 5: 三条曲线同时绘制

当然，如果你不想写那么多代码的话，也可以直接使用plot_ROC2()函数直接绘制，是不是简单的多。

plot_ROC2(yvar="group",xvars=c("a","b","c"),dataname="ROC")

Figure 6: 三条曲线同时绘制简单函数

分面显示

将三张图放在一起，可以看到数值重叠影响颜值，我们可以用ggplot2的facet进行分面绘制。

plot_ROC(list(a,b,c),facet=TRUE)

Figure 7: 三条曲线的分面显示

可以发现分面的标签默认是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添加标题，还可以更换主题等等。。。