Inverse Probability of Censoring Weighting (IPCW) estimation of Cumulative/Dynamic time-dependent ROC curve。
一 、用法:
timeROC(T, delta, marker, other_markers = NULL, cause, weighting = "marginal", times, ROC = TRUE, iid = FALSE)
T:事件时间
delta:事件状态. 删失数据编码为0.
marker :计算ROC的biomaker,默认是marker值越大,事件越可能发生;反之的话,前面加-号。
other_markers:矩阵形式输入,可多个marker,类似协变量. 默认值other_markers=NULL.
cause:所关心的事件结局。没有竞争风险(Without competing risks)中,必须是非删失数据的编码方式,一般为1。
存在竞争风险(With competing risks)中,和所关心的事件结局一致,通常为1 or 2.
weighting:计算方法,默认是weighting="marginal",KM模型;weighting="cox" 和weighting="aalen"分别为COX模型和additive Aalen 模型。
times:想计算的ROC曲线的时间节点。
ROC:默认值ROC = TRUE,保存sensitivities 和 specificties值。
iid: 默认值iid = FALSE。iid = TRUE 才会保存置信区间,但是样本量大了后,耗时耗资源。
二、Time-Dependent ROC 实现
2.1 -------------Without competing risks-------------------
载入R包和数据
library(timeROC) library(survival) data(pbc) head(pbc) id time status trt age sex ascites hepato spiders edema bili chol albumin copper alk.phos 1 1 400 2 1 58.76523 f 1 1 1 1.0 14.5 261 2.60 156 1718.0 2 2 4500 0 1 56.44627 f 0 1 1 0.0 1.1 302 4.14 54 7394.8 3 3 1012 2 1 70.07255 m 0 0 0 0.5 1.4 176 3.48 210 516.0 4 4 1925 2 1 54.74059 f 0 1 1 0.5 1.8 244 2.54 64 6121.8 5 5 1504 1 2 38.10541 f 0 1 1 0.0 3.4 279 3.53 143 671.0 6 6 2503 2 2 66.25873 f 0 1 0 0.0 0.8 248 3.98 50 944.0 pbc<-pbc[!is.na(pbc$trt),] # select only randomised subjects pbc$status<-as.numeric(pbc$status==2) # create event indicator: 1 for death, 0 for censored
1) Kaplan-Meier estimator(default)
ROC.bili.marginal<-timeROC(T=pbc$time, delta=pbc$status,marker=pbc$bili, cause=1,weighting="marginal", times=quantile(pbc$time,probs=seq(0.2,0.8,0.1)), iid=TRUE) ROC.bili.marginal
结果如下:
Time-dependent-Roc curve estimated using IPCW (n=312, without competing risks). Cases Survivors Censored AUC (%) se t=999.2 53 249 10 83.96 2.91 t=1307.4 68 218 26 85.66 2.56 t=1839.5 86 156 70 88.03 2.25 t=2555.7 102 94 116 83.41 3.17 t=3039 108 63 141 80.79 3.48
AUC : vector of time-dependent AUC estimates at each time points.
TP : matrix of time-dependent True Positive fraction (sensitivity) estimates.
FP : matrix of time-dependent False Positive fraction (1-specificity) estimates.
2) Cox model
###可添加协变量 (with covariates bili, chol and albumin) ROC.bili.cox<-timeROC(T=pbc$time, delta=pbc$status,marker=pbc$bili, other_markers=as.matrix(pbc[,c("chol","albumin")]), cause=1,weighting="cox", times=quantile(pbc$time,probs=seq(0.2,0.8,0.1))) ROC.bili.cox
2.2 -------------With competing risks-------------------
data(Paquid) ##Example with Paquid data # evaluate DDST cognitive score as a prognostic tool for dementia onset, accounting for death without dementia competing risk. ROC.DSST<-timeROC(T=Paquid$time,delta=Paquid$status, marker=-Paquid$DSST,cause=1, weighting="cox", other_markers=as.matrix(Paquid$MMSE), times=c(3,5,10),ROC=TRUE) ROC.DSST Time-dependent-Roc curve estimated using IPCW (n=2561, with competing risks). Cases Survivors Other events Censored AUC_1 (%) AUC_2 (%) t=3 70 2117 194 180 80.83 79.85 t=5 122 1834 313 292 79.55 77.65 t=10 318 1107 545 591 76.40 71.93
注1:竞争风险 (见文末)
三、绘制ROC曲线
plot(ROC.DSST,time=5) plot(ROC.DSST,time=3,add=TRUE,col="blue") plot(ROC.DSST,time=10,add=TRUE,col="grey50") legend("bottomright",c("Y-5","Y-3","Y-10"),col=c("red","blue","grey50"),lty=1,lwd=2)
四、输出置信区间
confint(object, parm=NULL, level = 0.95,n.sim=2000, ...)
注:object 必须是 timeROC function 得到的,且参数 weighting="marginal" , iid = TRUE.
#confint(ROC.bili.marginal)$CB_AUC confint(ROC.bili.marginal)$CI_AUC 2.5% 97.5% t=999.2 78.27 89.66 t=1307.4 80.64 90.68 t=1548.4 83.31 92.15 ... ... t=3039 73.96 87.62
confint(ROC.bili.cox) ##iid=FALSE (default) 所以报错
注2:见文章结尾。
五、比较两个Time-Dependent AUC
compare(x, y, adjusted = FALSE, abseps = 1e-06)
其中 x y 都必须是timeROC function ,使用 weighting="marginal" 和 iid = TRUE参数得到的。
ROC.albumin<-timeROC(T=pbc$time, delta=pbc$status,marker=-pbc$albumin, cause=1,weighting="marginal", times=quantile(pbc$time,probs=seq(0.2,0.8,0.1)), iid=TRUE)
we compare albumin and bilirubin as prognostic biomarkers.
compare(ROC.albumin,ROC.bili.marginal) #compute p-values of comparison tests
当adjusted=TRUE 会输出矫正后的P值,以及相关系数矩阵
compare(ROC.albumin,ROC.bili.marginal,adjusted=TRUE) $p_values_AUC t=999.2 t=1307.4 t=1548.4 t=1839.5 t=2234.2 t=2555.7 t=3039 Non-adjusted 0.02294328 0.01568158 0.0006036388 0.001167184 0.02265649 0.06409548 0.2022485 Adjusted 0.09739715 0.06915425 0.0032942671 0.005991873 0.09638428 0.23790873 0.5874903
相关 系数矩阵
$Cor 略
六 绘制Time-Dependent AUC 曲线
plotAUCcurve(object, FP = 2, add = FALSE, conf.int = FALSE, conf.band = FALSE, col = "black")
当weighting="marginal" and iid = TRUE的前提下,可以添加置信区间(两种)
conf.int:the bands of pointwise confidence intervals
conf.band:the simultaneous confidence bands.
# # plot AUC curve for albumin only with pointwise confidence intervals and simultaneous confidence bands plotAUCcurve(ROC.albumin,conf.int=TRUE,conf.band=TRUE)
# # plot AUC curve for albumin and bilirunbin with pointwise confidence intervals plotAUCcurve(ROC.albumin,conf.int=TRUE,col="red") plotAUCcurve(ROC.bili.marginal,conf.int=TRUE,col="blue",add=TRUE) legend("bottomright",c("albumin","bilirunbin"),col=c("red","blue"),lty=1,lwd=2)
七、参考资料
https://www.rdocumentation.org/packages/timeROC/
八、注释
注1:竞争风险指研究对象除了会出现研究者感 兴趣的结局(如CSS) ,还会出现如车祸等其他意外结局,它的出现会导致感兴趣的事件永远不会 发生(这被认为是与右删失数据(right2censord data)的最大差别) ,即出现了竞争。
注2:若是With competing risks,结果中会分1 2 :
CI_AUC_1 : a matrix.the pointwise confidence intervals of AUC with definition (i) of controls.
CB_AUC_1 : a matrix. the simultaneous confidence band of the AUC curve with definition (i) of controls.
C.alpha.1 : a numeric value corresponding to the quantile required for simultaneous confidence bands computation CB_AUC_1 (estimated by simulations).
注3:definition (i)和definition (ii)的区分如下:
With competing risks, two definitions of controls were suggested: (i) a control is defined as a subject i that is free of any event, i.e with Ti > t, and (ii) a control is defined as a subject i that is not a case, i.e with Ti > t or with Ti ≤ t and δi 6= 1.