R|timeROC-分析 - 码农教程

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.0pbc<-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.