ROC曲线

1. ROC曲线是反应一个分类器在某个阈值时对样本的识别能力。
2. 可以借助ROC曲线选择出某一诊断方法最佳的诊断界限值。ROC曲线越是靠近左上角，试验的FPR越高和FPR越低，即灵敏度越高，误判率越低，则诊断方法的性能越好。

1. 下面的这个方法做ROC曲线适用于只有一个因素标签，如有无疗效，是否耐药等，不适合用于生存的ROC制作(两个因素time and censor)。

1. 导入包和数据

library('plotROC')
test=read.table("test.txt",sep="\t",header=T,check.names=F)

test.txt
2.plotROC提供的函数melt_roc()可以将多个变量列变为长格式，方便数据的绘制：

longtest <- melt_roc(test, "survival", c("CCR2", "UCHL1"))
head(longtest)

longtest

3. 画比较图

ggplot(longtest, aes(d = D, m = M, color = name)) + geom_roc() + style_roc()

image.png

4. 自定义函数

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

** 5. 绘图**

plotROC(longtest, predict_col = M, target = D, group = name, positive = 1)
# 参数1：提供数据框
# 参数2：提供预测数值列
# 参数3：提供二分类信息列（尽量为0-1，字符也可以）
# 参数4：提供一个组别
# 参数5：这里1表示成功，如果target是success和failure，可以知道positive="success"
# 注意，这里只有3条曲线绘制时才会给出AUC在图上，可以修改函数进行自定

image.png

参考资料

[roc曲线分析和可视化] https://blog.csdn.net/u013066730/article/details/96476546

2. 生存资料的结果涉及生存状态和生存时间两个变量，所有不能用普通的ROC曲线，必须用时间依赖性ROC曲线，也就是survivalROC，这样才能把两个因素都分析进去。

1. 导入数据

Sys.setlocale('LC_ALL','C')
library(survivalROC)
CCR2=read.table("CCR2.txt",sep="\t",header=T,check.names=F)

CCR2.txt

2. 计算相应的TP和FP

nobsCCR2 <- NROW(CCR2)
cutoff <- 1095 #这里就是时间节点的设置

Mayo4.1= survivalROC(Stime=CCR2$time,##生存时间
                     status=CCR2$censor,## 终止事件    
                     marker = CCR2$CCR2, ## marker value    
                     predict.time = cutoff,## 预测时间截点
                     method = 'KM')##span,NNE法的namda
str(Mayo4.1)## list结构,是每一个marker的cutoff值都计算出相应的TP,FP

Mayo4.2= survivalROC(Stime=CCR2$time,  
                     status=CCR2$censor,      
                     marker = CCR2$UCHL1,     
                     predict.time = cutoff, method="KM")

Mayo4.3= survivalROC(Stime=CCR2$time,  
                     status=CCR2$censor,      
                     marker = CCR2$combine,     
                     predict.time = cutoff, method="KM")

Mayo4.1

3. 画图

require(ggsci)
library("scales")
pal_nejm("default")(8)
show_col(pal_nejm("default")(8))

plot(Mayo4.1$FP, Mayo4.1$TP, ## x=FP,y=TP
     type="l",col="#BC3C29FF", ##线条设置
     xlim=c(0,1), ylim=c(0,1), 
     xlab=("1 - Specificity"), ##连接
     ylab="Sensitivity",
     main="Time dependent ROC")## \n换行符
abline(0,1,col="gray",lty=2)##线条颜色

#lines函数在原有基础上继续绘图 #E18727FF
#legend函数增加legend
lines(Mayo4.2$FP, Mayo4.2$TP, type="l",col="#0072B5FF",xlim=c(0,1), ylim=c(0,1))
lines(Mayo4.3$FP, Mayo4.3$TP, type="l",col="#E18727FF",xlim=c(0,1), ylim=c(0,1))
legend(0.45,0.3,c(paste("AUC of CCR2 =",round(Mayo4.1$AUC,3)),
                 paste("AUC of UCHL1 =",round(Mayo4.2$AUC,3)),
                 paste("AUC of Combine =",round(Mayo4.2$AUC,3))),
       x.intersp=1, y.intersp=0.8, #调整lengend在图中的位置
       lty= 1 ,lwd= 2,col=c("#BC3C29FF","#0072B5FF",'#E18727FF'),
       bty = "n",# bty框的类型
       seg.len=1,cex=0.8)# 

survivalROC