打开APP
userphoto
未登录

开通VIP,畅享免费电子书等14项超值服

开通VIP
生存分析模型的时依ROC曲线绘制
userphoto

2023.08.26 上海

关注
受试者工作特征曲线(Receiver Operating Characteristic Curve, ROC曲线)是将数据进行二分类时,根据不同的界值,以真阳性率(TPR灵敏度)为纵坐标,假阳性率(FPR1-特异度)为横坐标绘制的曲线ROC曲线下的面积(Area Under ROC CurveAUC常用于诊断试验的评估,在模型评估中也常作为模型效果(区分度)的一个指标,比如众所周知的在二分类的logistic回归中应用。

标准的ROC分析中结局和预测变量都是固定不变的,操作可参见《ROC分析》。但在生存分析中,受试者的生存结局是与时间相关的(time-to-event outcomes),结局事件的发生会随着时间的变化而变化,随着随访时间的延长越来越多的受试者会发生失效事件,也有一些受试者开始失访(截尾),在随访开始时测量的预测变量与结局的相关性也可能会变的越来越低,甚至有些解释变量的取值也在长期的随访过程中发生了变化。时间依赖性(时依)ROCtime-dependent ROC)可以更全面地描述这种生存数据分析模型的评估情况。

删失时间数据的时依灵敏度和特异度有多种定义,常见的有:

1Cumulative sensitivity and dynamic specificity (C/D),采用累积病例(Cumulative case)来计算敏感性、动态对照(Dynamic control)来计算特异性。对于任意时间t任何一个个体都可以按照其在时间t的状态被划分为病例组或对照组。C/D中的病例组定义为在基线t=0和时间t之间发生失效事件的个体,而对照组被定义为在时间t时仍未发生失效事件的个体。根据定义,受试者的状态(病例/对照)会随着时间的变化而变化。相比I/DI/SC/D具有更多的临床相关性,在临床中应用更广泛。

2Incident sensitivity and dynamic specificity (I/D)。采用事件病例(Incident case)来计算敏感性、动态对照(Dynamic control)来计算特异性。这里的事件组定义为在时间t时刻发生失效事件的个体,对照组被定义同C/D

3Incident sensitivity and static specificity (I/S)。采用事件病例(Incident case)来计算敏感性、静态对照(Static control)来计算特异性。事件组定义同I/D,对照组是经过固定的随访期间没有发生失效事件个体。

Closed circles indicate individuals who had an event, open circles indicate individuals who had censored event-times.C/D: A, B and E are cases and C, D and F are controls; I/D: Only A is the case and C, D and F are controls; I/S: Only A is the case and D and F are controls.

BMC Medical Research Methodology (2017):17-53.

R中可实现时依ROC的包有很多,比如timeROCsurvivalROCtdROCrisksetROC等。

示例:299名心力衰竭患者的医疗记录。
数据下载地址:
https://archive.ics.uci.edu/dataset/519/heart+failure+clinical+records
我们只选择AgeHBPAnaemiaCPKEF几个解释变量进行演示。

1】数据导入及因子设置

library(bruceR)

HF<- import("D:/Temp/HF.xlsx")

HF$HBP<-factor(HF$HBP,levels=c(0,1),labels=c("No","Yes"))

HF$Anaemia<-factor(HF$Anaemia,levels=c(0,1),labels=c("No","Yes"))

2】计算多因素效应的指标
#将多因素的效应合并为一个值。本例用Cox回归估计的预后指数作为多因素的合并效应,对应后续函数中的marker值。

library(survival)

Cox.hf<-coxph(Surv(Time,DEATH==1)~Age+HBP+Anaemia+CPK+EF,data=HF,ties="breslow")

PI<-Cox.hf$linear.predictor
#PI<-predict(Cox.hf, newdata=HF,type="lp") #在预测模型估计验证集和/或测试集的ROC时可使用predict函数来计算marker
3.1timeROC包计算C/D,并绘制相应的时依ROC
#采用逆概率删失加权法(Inverse Probability of Censoring WeightingIPCW)估计累积/动态(C/D)时依ROC曲线,权重计算方法可选择Kaplan-Meier法、Cox模型和Aalen加性风险模型。

library(timeROC)

tROC.hf<-timeROC(T=HF$Time,delta=HF$DEATH,marker=PI,cause=1,weighting="marginal",times=c(30.44,91.31,182.63),iid=TRUE)
#T:生存时间;
#delta:失效事件,需要注意删失必须用0表示;
#marker:感兴趣的预测/解释变量。如果只是考察单个因素对结局的预测效果,直接使用原始值即可;多因素可以用上一步计算的PI值。该值越大与较高的事件风险相关,如果值越小与较高的事件风险相关,则需要将该值加负号;
#other_markers:矩阵形式,对计算逆概率删失权重计算有影响的变量,类似需要控制的变量;
#cause:失效事件的指示值,默认为1。在竞争风险模型中需要在此正确指定感兴趣的失效事件;
#weighting:计算权重的方法。默认weighting="marginal" 采用Kaplan-Meier法估计删失值分布,而weighting="cox"weighting="aalen"分别采用Cox模型和Aalen加性风险模型;
#times:感兴趣的时间点,函数将估计该时间点上的时依ROC。本例一年按365.25天计算,每月平均365.25/12=30.44天;
#ROC:是否保存敏感度和特异度值的估计值;
#iid:时依ROC曲线下面积估计时的一些涉及独立同分布的指标,默认FALSE。如需要用到置信区间、AUC比较时,比如后续可能会用plotAUCcurve绘制带有置信区间的时依AUC变化图,用compare函数比较两条时依ROCAUC,则需要设置为TRUE。样本量较大(超过2000)时比较耗时。说得挺复杂,个人理解其实就是是否需要计算标准误。

#绘制多条时依ROC于一图

plot(tROC.hf, time=30.44, col = "red", title=FALSE,add = FALSE)
plot(tROC.hf, time=91.31, col = "blue", add =TRUE)
plot(tROC.hf, time=182.63, col = "green", add =TRUE)

#添加标题和图例

title(main="Time Dependent ROC")
legend("bottomright",paste0("AUC of ",c(1,3,6),"month survival:",format(round(tROC.hf$AUC,4),nsmall=2)), col=c("red","blue","green"),lty=1,lwd=2,bty = "n")

#绘制AUC随时间的变化趋势
plotAUCcurve(tROC.hf, conf.int =TRUE,col="blue")

3.2survivalROC包绘制相应的时依ROC

#感兴趣时间点处删失分布的估计方法可选KM法(KaplanMeier)或最邻近估计法( Nearest Neighbor EstimationNNE),采用NNE时需要指定平滑参数lanbda或者span
library(survivalROC)

tROC.hf2.km<-survivalROC(Stime=HF$Time,status=HF$DEATH,marker=PI,predict.time=30.44,method = "KM")

tROC.hf2.nne1<-survivalROC(Stime=HF$Time,status=HF$DEATH,marker=PI,predict.time=91.31,lambda=0.001)
tROC.hf2.nne2<-survivalROC(Stime=HF$Time,status=HF$DEATH,marker=PI,predict.time=182.63,span = 0.25*(NROW(HF))^(-0.20))
plot(tROC.hf2.km$FP, tROC.hf2.km$TP, type="l", col="green",xlim=c(0,1), ylim=c(0,1), xlab="FP:1-Specificity",ylab="TP:sensitivity",main="HF Death Cox Model, 1/3/6Month \n Method = KM/NNE")
abline(0,1,col="red",lty=2)

#绘制多条时依ROC于一图

lines(tROC.hf2.nne1$FP, tROC.hf2.nne1$TP, type="l",col="blue",xlim=c(0,1), ylim=c(0,1))
lines(tROC.hf2.nne2$FP, tROC.hf2.nne2$TP, type="l",col="yellow",xlim=c(0,1), ylim=c(0,1))

#添加图例

legend(0.55,0.15,c(paste("AUC of 1 month: ",round(tROC.hf2.km$AUC,3)), paste("AUC of 3 month: ",round(tROC.hf2.nne1$AUC,3)), paste("AUC of 6 month : ",round(tROC.hf2.nne2$AUC,3))), lty= 1 ,lwd= 2,col=c("green","blue",'yellow'),bty = "n", x.intersp=1, y.intersp=1,cex=0.8,seg.len=1) 

3.3tdROC包:tdROC函数用非参数加权校正计算灵敏度、特异度及ROC曲线下面积,plot.tdROC绘制时依ROC
library(tdROC)

tROC.hf3<-tdROC(X=PI,Y=HF$Time,delta=HF$DEATH,tau=30.44,nboot=100)

plot.tdROC(tROC.hf3, lwd = 2, xlab = "1-specificity", ylab = "sensitivity", xlim = c(0, 1), ylim = c(0, 1), col = "black", main = "Time Dependent ROC", abline = T)
如想把其他感兴趣的时间点的ROC绘制在同一图上,可使用lineslegend函数来完成,可参考survivalROC包。

3.4risksetROC包计算I/D,并绘制相应的时依ROC

library(risksetROC)

tROC.hf4.6<-risksetROC(Stime=HF$Time,status=HF$DEATH,marker=PI,predict.time=182.63,method="Cox",col="green",lty=1,lwd=1, xlab="FP:1-Specificity",ylab="TP:sensitivity",main="Time Dependent ROC")

#method:默认"cox"假设模型为比例风险模型,非比例风险模型可选"LocalCox""Schoenfeld" ,两者的平滑方法不一样
tROC.hf4.3<-risksetROC(Stime=HF$Time,status=HF$DEATH,marker=PI,predict.time=91.31,method="Cox",plot=FALSE)
tROC.hf4.1<-risksetROC(Stime=HF$Time,status=HF$DEATH,marker=PI,predict.time=30.44,method="Cox",plot=FALSE)

#绘制多条时依ROC于一图

lines(tROC.hf4.3$FP, tROC.hf4.3$TP, type="l",col="blue",xlim=c(0,1), ylim=c(0,1))
lines(tROC.hf4.1$FP, tROC.hf4.1$TP, type="l",col="red",xlim=c(0,1), ylim=c(0,1))
#添加图例
legend("bottomright",c(paste("AUC of 1 month: ",round(tROC.hf4.1$AUC,3)), paste("AUC of 3 month: ",round(tROC.hf4.3$AUC,3)), paste("AUC of 6 month : ",round(tROC.hf4.6$AUC,3))),col=c("red","blue","green"),lty=1,lwd=2,bty = "n") 
在计算C/D时病例组是感兴趣时间点之前所有发生结局事件的患者,早期发生结局事件的高风险患者一直为后续时间上的C/D计算做贡献,随着感兴趣的的时间点的延长,C/D时依ROCAUC会趋于稳定。与此不同,计算I/D时的病例组只是感兴趣的时间点发生的结局事件的患者,也就是只有在每个时间点处于风险集中的个人才能贡献数据,随着感兴趣的的时间点的延长,I/D时依ROCAUC会趋于变小,也许是因为生存时间足够长的患者与基线变量关系不大,如果没有结局事件发生,ROC基本上就是一条直线了。

出处:https://datascienceplus.com/time-dependent-roc-for-survival-prediction-models-in-r/

END
本站仅提供存储服务,所有内容均由用户发布,如发现有害或侵权内容,请点击举报
打开APP,阅读全文并永久保存 查看更多类似文章
猜你喜欢
类似文章
【热】打开小程序,算一算2024你的财运
基于BP神经网络的跨既有线高速铁路桥梁施工安全风险评估
竞争风险数据的统计分析方法 (PART Ⅲ)
钠盐摄入与老年人的10年死亡率、心脑血管疾病(CVD)和心力衰竭(HF)事件发生无关。
R语言解读自回归模型
R语言多项式线性模型:最大似然估计二次曲线
生存分析与R
更多类似文章 >>
生活服务
热点新闻
分享 收藏 导长图 关注 下载文章
绑定账号成功
后续可登录账号畅享VIP特权!
如果VIP功能使用有故障,
可点击这里联系客服!

联系客服