欢迎来到医科研，这里是白介素2的读书笔记，跟我一起聊临床与科研的故事, 生物医学数据挖掘，R语言，TCGA、GEO数据挖掘。
含有删失数据的生存数据，使用survivalROC包
包括Kaplan-Meier (KM) or Nearest Neighbor Estimation (NNE) 两种方法
假设我们有删失的生存数据与基线marker值，我们希望看到marker如何预测数据集中的受试者的存活时间。特别是，假设我们有几天的生存时间，我们想看看标记如何预测一年的存活（predict.time=365）。该功能roc.km.calc()返回感兴趣的时间点的唯一标记值、TP（真阳性）、FP（假阳性）、对应于感兴趣时间点(predict.time)和AUC(ROC)曲线下面积的Kaplan-Meier生存估计。
返回值：
cut.values：由于计算TP和FP的marker值
TP：根据cutoff 判断的TRUE postive 真阳性
FP：根据cutoff判断的假阳性
predict.time：感兴趣的时间截点：可以为5年，3年等等
Survival：kaplan-Meier法的预估生存时间
AUC：Area under ROC，在时间截点的曲线下面积
实际代码演示
1Sys.setlocale('LC_ALL','C')
2library(survivalROC)
3data(mayo)
4head(mayo)
5## time censor mayoscore5 mayoscore4
6## 1 41 1 11.251850 10.629450
7## 2 179 1 10.136070 10.185220
8## 3 334 1 10.095740 9.422995
9## 4 400 1 10.189150 9.567799
10## 5 130 1 9.770148 9.039419
11## 6 223 1 9.226429 9.033388
12nobs <- NROW(mayo)
13cutoff <- 365
marker value可以为表达值，免疫分数，检验参数等任何可以定义为数值的指标
Mayo4.1得到的是列表，其内容是每一个marker的cutoff值都计算出相应的TP,FP
1Mayo4.1= survivalROC(Stime=mayo$time,##生存时间
2 status=mayo$censor,## 终止事件
3 marker = mayo$mayoscore4, ## marker value
4 predict.time = cutoff,## 预测时间截点
5 span = 0.25*nobs^(-0.20))##span,NNE法的namda
6str(Mayo4.1)## list结构
7## List of 6
8## $ cut.values : num [1:313] -Inf 4.58 4.9 4.93 4.93 ...
9## $ TP : num [1:313] 1 0.997 0.995 0.993 0.99 ...
10## $ FP : num [1:313] 1 0.997 0.994 0.99 0.987 ...
11## $ predict.time: num 365
12## $ Survival : num 0.929
13## $ AUC : num 0.931
14## 绘图
15plot(Mayo4.1$FP, Mayo4.1$TP, ## x=FP,y=TP
16 type="l",col="red", ##线条设置
17 xlim=c(0,1), ylim=c(0,1),
18 xlab=paste( "FP", "\n", "AUC = ",round(Mayo4.1$AUC,3)), ##连接
19 ylab="TP",
20 main="Mayoscore 4, Method = NNE \n Year = 1")## \n换行符
21abline(0,1,col="gray",lty=2)##线条颜色
MAYOSCORE 4作为marker, KM法计算
1## MAYOSCORE 4, METHOD = KM
2Mayo4.2= survivalROC(Stime=mayo$time,
3 status=mayo$censor,
4 marker = mayo$mayoscore4,
5 predict.time = cutoff, method="KM")
6
7plot(Mayo4.2$FP, Mayo4.2$TP,
8 type="l",col="red",xlim=c(0,1), ylim=c(0,1),
9 xlab=paste( "FP", "\n", "AUC = ",round(Mayo4.2$AUC,3)),
10 ylab="TP",
11 main="Mayoscore 4, Method = KM \n Year = 1")
12abline(0,1,col="gray",lty=2)
lines函数在原有基础上继续绘图
legend函数增加legend
1## NNE法
2plot(Mayo4.1$FP, Mayo4.1$TP, ## x=FP,y=TP
3 type="l",col="red", ##线条设置
4 xlim=c(0,1), ylim=c(0,1),
5 xlab=("FP"), ##连接
6 ylab="TP",
7 main="Time dependent ROC")## \n换行符
8abline(0,1,col="gray",lty=2)##线条颜色
9
10## KM法
11lines(Mayo4.2$FP, Mayo4.2$TP, type="l",col="green",xlim=c(0,1), ylim=c(0,1))
12legend(0.6,0.2,c(paste("AUC of NNE =",round(Mayo4.1$AUC,3)),
13 paste("AUC of KM =",round(Mayo4.2$AUC,3))),
14 x.intersp=1, y.intersp=0.8,
15 lty= 1 ,lwd= 2,col=c("red","green"),
16 bty = "n",# bty框的类型
17 seg.len=1,cex=0.8)#
1require(ggsci)
2library("scales")
3pal_nejm("default")(8)
4show_col(pal_nejm("default")(8))
5
6## NNE法
7plot(Mayo4.1$FP, Mayo4.1$TP, ## x=FP,y=TP
8 type="l",col="#BC3C29FF", ##线条设置
9 xlim=c(0,1), ylim=c(0,1),
10 xlab=("FP"), ##连接
11 ylab="TP",
12 main="Time dependent ROC")## \n换行符
13abline(0,1,col="gray",lty=2)##线条颜色
14
15## KM法
16lines(Mayo4.2$FP, Mayo4.2$TP, type="l",col="#0072B5FF",xlim=c(0,1), ylim=c(0,1))
17legend(0.6,0.2,c(paste("AUC of NNE =",round(Mayo4.1$AUC,3)),
18 paste("AUC of KM =",round(Mayo4.2$AUC,3))),
19 x.intersp=1, y.intersp=0.8,
20 lty= 1 ,lwd= 2,col=c("#BC3C29FF","#0072B5FF"),
21 bty = "n",# bty框的类型
22 seg.len=1,cex=0.8)#
很显然这样的配色让图上了一个档次。
联系客服
