之前的timeROC曲線圖在：(本來(lái)這里應(yīng)該有個(gè)鏈接畸冲，可是已經(jīng)被封了煮仇，又不讓放別的平臺(tái)的鏈接劳跃，所以放不了了。)
剛剛對(duì)它進(jìn)行了一點(diǎn)調(diào)整和升級(jí)浙垫。
1.輸入數(shù)據(jù)
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列刨仑，生存時(shí)間、結(jié)局事件和預(yù)測(cè)值夹姥，下面的代碼里 time 杉武，event 和 fp 就是這三列的列名。
2.完成timeROC分析和畫(huà)圖
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()
做的調(diào)整有：
1.修改顏色辙售、去掉格子轻抱、橫縱坐標(biāo)的expand
2.更改圖例位置、文字旦部，加框
3.修改橫縱坐標(biāo)名字
4.簡(jiǎn)化了生成作圖數(shù)據(jù)的代碼
5.改正了四舍五入末尾0被省略掉的問(wèn)題祈搜，原來(lái)的0.90會(huì)被改成0.9，現(xiàn)在永遠(yuǎn)保留兩位士八。
3.關(guān)于橫縱坐標(biāo)
特異度與靈敏度容燕，是評(píng)價(jià)模型的兩個(gè)重要指標(biāo)
按照疾病診斷，有病是陽(yáng)性婚度，沒(méi)病是陰性來(lái)說(shuō)：
FPR是誤診率蘸秘，把沒(méi)病的人說(shuō)成有病的概率，是1-特異度蝗茁，1-沒(méi)病的人診斷沒(méi)病的概率
TPR是靈敏度醋虏，把有病的人診斷出有病的概率。
我們希望靈敏度高评甜，誤診率低灰粮，然而他倆負(fù)相關(guān)，就像魚(yú)與熊掌不可兼得忍坷。好在不是必須二選一粘舟，而是可以選個(gè)平衡點(diǎn)哦。
FPR是1 - Specificity佩研，TPR是Sensitivity柑肴。
我看到了一個(gè)非常清楚的的解讀：http://www.reibang.com/p/7919ef304b19
4.根據(jù)ROC曲線確定最佳截點(diǎn)
搜了好一圈沒(méi)有看到timeROC包怎樣去找最佳截點(diǎn)，倒是清一色都是survivalROC包的結(jié)果旬薯。似乎并不像timeROC一樣支持同時(shí)計(jì)算好幾個(gè)時(shí)間節(jié)點(diǎn)晰骑，一次只能計(jì)算一個(gè)時(shí)間。我想確定一下兩個(gè)R包算出的結(jié)果是否一致，所以把一年的ROC曲線花在了一起硕舆，可見(jiàn)二者幾乎完全重疊秽荞。
因此我用survivalROC的結(jié)果來(lái)找截點(diǎn)，是完全可以的抚官。
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
不得不說(shuō)一下上面這句代碼扬跋，約登指數(shù)=靈敏度+特異度?1，即=TP?FP凌节，約登指數(shù)最大時(shí)的截?cái)嘀导礊樽罴呀財(cái)嘀怠?/p>
5.獲取最佳截?cái)嘀盗硪环N的方法
我已經(jīng)用過(guò)很多次了钦听，survminer包里的surv_cutpoint函數(shù)，選出讓高低兩組間差異最顯著的截?cái)嘀怠?/p>
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).
用它來(lái)計(jì)算的結(jié)果倍奢，沒(méi)有考慮到時(shí)間因素朴上，因此與上面得到的結(jié)果不相同。
library(survminer)
res.cut <- surv_cutpoint(new_dat, time = "time", event = "event",
variables = "fp")
res.cut
## cutpoint statistic
## fp 11.18487 10.27772
我想了一下卒煞，是因?yàn)閠imeROC的計(jì)算中痪宰，生存時(shí)間超過(guò)1年，最終結(jié)局為死亡的病人畔裕，在1年時(shí)生存狀態(tài)為活著酵镜。當(dāng)我把這部分病人的生存狀態(tài)改為0（活著），這兩個(gè)方法計(jì)算出來(lái)的結(jié)果就相同咯
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
用這個(gè)截?cái)嘀底龇纸M柴钻，看看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)
哈哈，非彻噶福可以贴届。
