timeROC的升级版全面解读

2021-06-24 本文已影响0人小洁忘了怎么分身

之前的timeROC曲线图在：(本来这里应该有个链接，可是已经被封了，又不让放别的平台的链接，所以放不了了。)
刚刚对它进行了一点调整和升级。

1.输入数据

library(timeROC)
library(survival)
load("new_dat.Rdata")
head(new_dat)

##   time event       fp
## 1 3817     0 14.63898
## 2  254     1 14.89388
## 3  345     1 13.05678
## 4  109     1 13.88350
## 5  164     1 15.83493
## 6  212     1 16.26325

需要有3列，生存时间、结局事件和预测值，下面的代码里 time ，event 和 fp 就是这三列的列名。

2.完成timeROC分析和画图

result <-with(new_dat, timeROC(T=time,
                     delta=event,
                     marker=fp,
                     cause=1,
                     times=c(365,1095,1825),
                     iid = TRUE))
#identical(c(result$TP[,1],result$TP[,2],result$TP[,3]),as.numeric(result$TP))
dat = data.frame(fpr = as.numeric(result$FP),
                 tpr = as.numeric(result$TP),
                 time = rep(as.factor(c(365,1095,1825)),each = nrow(result$TP)))

library(ggplot2)
ggplot() + 
  geom_line(data = dat,aes(x = fpr, y = tpr,color = time),size = 1) + 
  scale_color_manual(name = NULL,values = c("#92C5DE", "#F4A582", "#66C2A5"),
                     labels = paste0("AUC of ",c(1,3,5),"-y survival: ",
                                     format(round(result$AUC,2),nsmall = 2)))+
  geom_line(aes(x=c(0,1),y=c(0,1)),color = "grey")+
  theme_bw()+
  theme(panel.grid = element_blank(),
        legend.background = element_rect(linetype = 1, size = 0.2, colour = "black"),
        legend.position = c(0.765,0.125))+
  scale_x_continuous(expand = c(0.005,0.005))+
  scale_y_continuous(expand = c(0.005,0.005))+
  labs(x = "1 - Specificity",
       y = "Sensitivity")+
  coord_fixed()

做的调整有：
1.修改颜色、去掉格子、横纵坐标的expand
2.更改图例位置、文字，加框
3.修改横纵坐标名字
4.简化了生成作图数据的代码
5.改正了四舍五入末尾0被省略掉的问题，原来的0.90会被改成0.9，现在永远保留两位。

3.关于横纵坐标

特异度与灵敏度，是评价模型的两个重要指标

按照疾病诊断，有病是阳性，没病是阴性来说：
FPR是误诊率，把没病的人说成有病的概率，是1-特异度，1-没病的人诊断没病的概率
TPR是灵敏度，把有病的人诊断出有病的概率。

我们希望灵敏度高，误诊率低，然而他俩负相关，就像鱼与熊掌不可兼得。好在不是必须二选一，而是可以选个平衡点哦。

FPR是1 - Specificity，TPR是Sensitivity。

我看到了一个非常清楚的的解读：https://www.jianshu.com/p/7919ef304b19

4.根据ROC曲线确定最佳截点

搜了好一圈没有看到timeROC包怎样去找最佳截点，倒是清一色都是survivalROC包的结果。似乎并不像timeROC一样支持同时计算好几个时间节点，一次只能计算一个时间。我想确定一下两个R包算出的结果是否一致，所以把一年的ROC曲线花在了一起，可见二者几乎完全重叠。

因此我用survivalROC的结果来找截点，是完全可以的。

library(survivalROC)
load("new_dat.Rdata")
head(new_dat)
##   time event       fp
## 1 3817     0 14.63898
## 2  254     1 14.89388
## 3  345     1 13.05678
## 4  109     1 13.88350
## 5  164     1 15.83493
## 6  212     1 16.26325
result <-with(new_dat, 
              survivalROC(Stime=time,
                          status=event,
                          marker=fp,
                          predict.time=365,
                          method="KM"))
## 
##  12 records with missing values dropped.
result$cut.values[which.max(result$TP-result$FP)]
## [1] 11.53272

不得不说一下上面这句代码，约登指数=灵敏度+特异度−1，即=TP−FP，约登指数最大时的截断值即为最佳截断值。

5.获取最佳截断值另一种的方法

我已经用过很多次了，survminer包里的surv_cutpoint函数，选出让高低两组间差异最显著的截断值。

Determine the optimal cutpoint for one or multiple continuous variables at once, using the maximally selected rank statistics from the 'maxstat' R package. This is an outcome-oriented methods providing a value of a cutpoint that correspond to the most significant relation with outcome (here, survival).

用它来计算的结果，没有考虑到时间因素，因此与上面得到的结果不相同。

library(survminer)
res.cut <- surv_cutpoint(new_dat, time = "time", event = "event",
   variables = "fp")
res.cut
## cutpoint statistic
## fp 11.18487  10.27772

我想了一下，是因为timeROC的计算中，生存时间超过1年，最终结局为死亡的病人，在1年时生存状态为活着。当我把这部分病人的生存状态改为0（活着），这两个方法计算出来的结果就相同咯

new_dat2 = new_dat
new_dat2$event[new_dat2$time>365 & new_dat2$event==1] = 0
res.cut <- surv_cutpoint(new_dat2, time = "time", event = "event",
   variables = "fp")
res.cut
##    cutpoint statistic
## fp 11.53272  7.327764

用这个截断值做分组，看看KM曲线

new_dat$risk = ifelse(new_dat$fp<res.cut$cutpoint$cutpoint,"low","high")
fit <- survfit(Surv(time, event)~risk, data=new_dat)
ggsurvplot(fit, data=new_dat, pval=TRUE,palette = "jco",
           risk.table = TRUE, conf.int = TRUE)

哈哈，非常可以。