之前的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)
哈哈，非常可以。
