简化版的ROC曲线

测试开头

简化版的ROC曲线今天是生信星球陪你的第863天

大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~

就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~

这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我!

之前的ROC曲线代码，用的是plot基础包画的，这一篇是用pROC包，无缝对接ggplot2，而且代码更加简单了。

1.输入数据

需要tumor表达矩阵exprSet、生存信息表格meta

load("TCGA-KIRC_sur_model.Rdata")
ls()
## [1] "exprSet" "meta"    "proj"
exprSet[1:4,1:4]
##            TCGA-B8-5162-01A TCGA-B4-5834-01A TCGA-CW-6087-01A TCGA-B0-4698-01A
## WASH7P           0.38273762       1.08187834       0.24083089       0.81791121
## MIR6859-1        0.00000000       0.03924379       0.00000000       0.09673435
## AL627309.1       0.02072238       0.03924379       0.00000000       0.41454776
## CICP27           0.00000000       0.03924379       0.01736828       0.09673435
head(meta)
##                                ID event     time  race age gender stage
## TCGA-B8-5162-01A TCGA-B8-5162-01A     0 1.200000 white  62   male    II
## TCGA-B4-5834-01A TCGA-B4-5834-01A     0 1.266667 white  59   male     I
## TCGA-CW-6087-01A TCGA-CW-6087-01A     1 1.366667 white  61   male    IV
## TCGA-B0-4698-01A TCGA-B0-4698-01A     1 1.400000 white  75   male    IV
## TCGA-B0-4690-01A TCGA-B0-4690-01A     1 1.433333 white  65   male    IV
## TCGA-AS-3778-01A TCGA-AS-3778-01A     0 1.433333 white  35   male     I

ROC曲线用于评估某个预测指标的表现，我们这里简单的用年龄和分期画一画，换成某个模型的预测值也是一样的用法。

# 把罗马数字换成阿拉伯数字，预测指标需要是一个连续型变量
meta$stage = as.numeric(as.factor(meta$stage))
head(meta$stage)
## [1] 2 1 4 4 4 
head(meta$age)
## [1] 62 59 61 75 65 35

2.计算并画图

library(pROC)
library(ggplot2)
#第一个指标
m <- roc(meta$event, meta$age)
auc(m)
## Area under the curve: 0.6295

#第二个指标
m2 <- roc(meta$event, meta$stage)
auc(m2)
## Area under the curve: 0.7526
#出图和美化
g <- ggroc(list(m = m,m2 = m2),legacy.axes = T,size = 1)
g + theme_minimal() +
  scale_color_manual(values = c("#2fa1dd", "#f87669"))+
  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), 
               colour = "grey", linetype = "dashed")+
  annotate("text",x = .75, y = .25,
           label = paste("AUC of m = ",format(round(as.numeric(auc(m)),2),nsmall = 2)),color = "#2fa1dd")+
  annotate("text",x = .75, y = .15,
           label = paste("AUC of m2 = ",format(round(as.numeric(auc(m2)),2),nsmall = 2)),color = "#f87669")