imd 统计教程 · 生存分析与竞争风险
← 教程首页 返回本模块

time-dependent ROC-生存数据ROC

聂志强实战医学统计2022-07-06 12:43:39广东


    在生存分析中，疾病状态和因素取值均会随时间发生变化。而标准的ROC曲线分析将个体的疾病状态和因素取值视作固定值，未将时间因素考虑在分析之中。在这种情况下，使用传统的ROC会忽略疾病状态或markers的时间依赖性，此时需用随时间变化的Time-dependent ROC。一旦应用了时间相关设置，就会在每个时间点观察疾病状态，从而在整个研究中产生不同的敏感性和特异性值。

一、原理

    简而言之，按不同时间点(通常固定1~5年)计算不同ROC。时间依赖性ROC定义：Ti表示疾病发作的时间，Xi是个体i的marker值（通常是基线值），X称为“marker”，X也可以表示从回归或其他模型计算的风险评分pre值。定义观察到的事件时间，Zi?=? min ( Ti ,? Ci )，其中Ci是censor时间，令δi是censor指标，如果事件（疾病）发生则取值为1，否则为0。令Di( t ) 为时间t的疾病状态，取值 1 或 0。对于给定的阈值cutoff (公式简写c)，时间相关的敏感性、特异性、AUC可以分别定义为：

上述公式计算不同cutoff下的敏感性和特异性我们即可得到时间t下的ROC曲线。个体被划分到病例组还是对照组会随着时间t取值的变化而发生变化。在此情形下，相应的ROC曲线也会发生变化。时间依赖性ROC曲线有(1) cumulative/dynamic(C/D)（即Cumulative sensitivity/dynamic specificity）、(2) incident/dynamic(I/D)和(3) incident/static(I/S)三种定义，其中cumulative/dynamic(C/D)定义适合生存分析，当前大多数研究使用此方法。对于任意时间t，每一个个体会按照其在时间t的状态被划分到病例组或对照组。在cumulative/dynamic中，如果一个个体在时间0和时间t之间发病，那么其会被划分到病例组（图1，A、B和E）;如果一个个体在时间0和时间t之间没有发病，那么其会被划分到对照组（图1，C、D和F）。

二、SCI案例

  1. Eur Heart J 2020，PMID：31620788

2.J Clin Oncol 2015，PMID：25605857

3.Diabetes Care 2019，PMID：30967432

三、R代码

#时间依赖性ROC，推荐 R中timeROC包 (可计算AUC95%CI，竞争风险)， 其他包还有：survivalROC(结果不一致), tdROC, timereg, risksetROC和survAUC。
#R4.1.3
#---------------  ovarian数据(survival包自带)  ---------------  
pacman::p_load(tidyverse,survival,survminer)
#data(ovarian)
ovarian <- ovarian
## 生存曲线720天后不变；单组1，多组group
ggsurvplot(survfit(Surv(futime, fustat) ~ 1, data = ovarian),
           risk.table = TRUE,break.time.by = 182.5)
#协变量(age, resid.ds, rx, ecog.ps)拟合Cox回归模型，并基于线性预测变量构建风险评分
# 拟合 Cox model
coxph1 <- coxph(formula = Surv(futime, fustat) ~ pspline(age, df = 4) +
                  factor(resid.ds) + factor(rx) + factor(ecog.ps),data = ovarian)
#Cox预测值 type=lp
ovarian$lp <- predict(coxph1, type = "lp")
#------推荐方法1---------  timeROC包  ---------------  
#可计算AUC95%CI和竞争风险，cox或km算法
pacman::p_load(timeROC,survival)
#K-M分析法，分析了在不同时间截点时，预测概率lp对患者的生存预测能力
#简化版data表达
time_roc_km <- with(ovarian, timeROC(
  T = futime,
  delta = fustat,
  marker = lp,  #lp 为cox预测值
  #默认的是更大的marker值意味着结果的更高风险；如果ROC反向，则反过来负号"-ovarian$lp“
  cause = 1,
  weighting="marginal", #marginal= km法，cox=cox法
  times = c(1 * 365, 2 * 365, 3 * 365), # 时间点，选取1年，2年和3年的生存率
  #也可设定%，times=quantile(dat$time,probs=seq(0.2,0.8,0.1)
  ROC = TRUE,other_markers = NULL,#保存TP和FP值
  iid = TRUE
))

#123年 AUC 点估计值
time_roc_km$AUC
#AUC 95%CI;confint()计函数计算95%的置信区间
confint(time_roc_km, level = 0.95)$CI_AUC
#区间2.5%  97.5% t=365  69.58 105.86

#cox分析法，校正cov(本例age)，分析了在不同时间截点时，预测概率lp对患者的生存预测能力
#weighting="marginal"时,iid为TRUE; weighting="cox"时,iid是默认的FALSE
time_roc_cox<-timeROC(T=ovarian$futime,
                      delta=ovarian$fustat,marker=ovarian$lp,
                      ROC = TRUE, #保存TP和FP值
                      other_markers=as.matrix(ovarian[,c("age")]),
                      cause=1,weighting="cox",
                      times = c(1 * 365, 2 * 365, 3 * 365),iid=FALSE)
#默认的是更大的marker值意味着结果的更高风险,反过来负号；
time_roc_cox
#原始图
plot(time_roc_km, time=1 * 365, col = "red", title = FALSE)  
plot(time_roc_km, time=2 * 365, add=TRUE, col="blue") 
plot(time_roc_km, time=3 * 365, add=TRUE, col="green") 
legend("bottomright",c("1 Years" ,"3 Years", "5 Years"),col=c("red", "blue", "green"), lty=1, lwd=2)

#美化方法1
#TP=sen FP=1-spec
time_ROC_df <- data.frame(
  TP_1year = time_roc_km$TP[, 1],
  FP_1year = time_roc_km$FP[, 1],
  TP_2year = time_roc_km$TP[, 2],
  FP_2year = time_roc_km$FP[, 2],
  TP_3year = time_roc_km$TP[, 3],
  FP_3year = time_roc_km$FP[, 3]
)
library(ggplot2)
#theme_bw()带格子；theme_classic 不带右上角；theme_test 带右上角
ggplot(data = time_ROC_df) +
  geom_line(aes(x = FP_1year, y = TP_1year), size = 1, color = "#BC3C29FF") +
  geom_line(aes(x = FP_2year, y = TP_2year), size = 1, color = "#0072B5FF") +
  geom_line(aes(x = FP_3year, y = TP_3year), size = 1, color = "#E18727FF") +
  geom_abline(slope = 1, intercept = 0, color = "grey", size = 1, linetype = 2) +
  theme_test( ) +
  annotate("text", x = 0.75, y = 0.25, size = 4.5,
           label = paste0("AUC at 1 years = ", sprintf("%.3f", time_roc_km$AUC[[1]])), color = "#BC3C29FF") +
  annotate("text",x = 0.75, y = 0.15, size = 4.5,
           label = paste0("AUC at 2 years = ", sprintf("%.3f", time_roc_km$AUC[[2]])), color = "#0072B5FF") +
  annotate("text",x = 0.75, y = 0.05, size = 4.5,
           label = paste0("AUC at 3 years = ", sprintf("%.3f", time_roc_km$AUC[[3]])), color = "#E18727FF") +
  labs(x = "False positive rate", y = "True positive rate") +
  theme(
    axis.text = element_text(face = "bold", size = 11, color = "black"),
    axis.title.x = element_text(face = "bold", size = 14, color = "black", margin = margin(c(15, 0, 0, 0))),
    axis.title.y = element_text(face = "bold", size = 14, color = "black", margin = margin(c(0, 15, 0, 0)))
  )
#美化方法2，推荐！0,0坐标紧贴XY轴
result <- with(ovarian, timeROC(
  T = futime,
  delta = fustat,
  marker = lp,  #lp 为cox预测值
  #默认的是更大的marker值意味着结果的更高风险；如果ROC反向，则反过来负号"-ovarian$lp“
  cause = 1,
  weighting="marginal", #marginal= km法，cox=cox法
  times = c(1 * 365, 2 * 365, 3 * 365), # 时间点，选取1年，2年和3年的生存率
  #也可设定%，times=quantile(dat$time,probs=seq(0.2,0.8,0.1)
  ROC = TRUE,other_markers = NULL,#保存TP和FP值
  iid = TRUE
))
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,3),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()


#---------第1年的ROC曲线阈值cutoff；指标youden --------
#方法1 dt数据
#提取各生存时间点的TPR、FPR值：
dt <- data.frame(FPR = as.numeric(time_roc_km$FP),
                 TPR = as.numeric(time_roc_km$TP),
                 time = rep(as.factor(c(365,2*365,3*365)), each = nrow(time_roc_km$TP)))
years1 <- dt[which(dt$time == 365),]
years1$youden <- as.numeric(years1$TPR) - as.numeric(years1$FPR)#TP-FP
cutoff1 <- ovarian$lp[which.max(years1$youden)] #YOUDEN最大值
cutoff1 
#方法2 推荐直接法；-0.04782439
ovarian$lp[which.max(time_ROC_df$TP_1year - time_ROC_df$FP_1year)] 
ovarian$lp[which.max(time_ROC_df$TP_2year - time_ROC_df$FP_2year)]
ovarian$lp[which.max(time_ROC_df$TP_3year - time_ROC_df$FP_3year)]

#---------survminer包里的surv_cutpoint函数，选出让time=1年，高低两组间差异最显著的截断值
library(survminer)
#因为timeROC的计算中，生存时间超过1年，最终结局为死亡的病人，在1年时生存状态为活着
ovarian2 = ovarian
ovarian2$fustat[ovarian2$futime>365 & ovarian2$fustat==1] = 0
res.cut <- surv_cutpoint(ovarian2, time = "futime", event = "fustat",variables = "lp")
res.cut
#1.78955
ovarian2$risk = ifelse(ovarian2$lp<res.cut$cutpoint$cutpoint,"low","high")
fit <- survfit(Surv(futime, fustat)~risk, data=ovarian2)
ggsurvplot(fit, data=ovarian2, pval=TRUE,palette = "jco",
           risk.table = TRUE, conf.int = TRUE)

#------compare函数进行多条ROC比较，并输出矫正后的P值和相关系数矩阵
compare(time_roc_km, time_roc_km2, adjusted = TRUE)

#------plotAUCcurve函数绘制不同时间节点的AUC曲线及其置信区间
plotAUCcurve(time_roc_km, conf.int=TRUE, col="red")
#plotAUCcurve(time_roc_km2, conf.int=TRUE, col="blue", add=TRUE)
legend("bottomright",c("lp", "lp2"), col = c("red","blue"), lty=1, lwd=2)
#------方法2---------  survivalROC包 ，与timeROC结果不同！有问题---------------  
pacman::p_load(survivalROC,dplyr)
#定义函数survivalROC_helper
survivalROC_helper <- function(t) {
  survivalROC(Stime        = ovarian$futime,
              status       = ovarian$fustat,
              marker       = ovarian$lp,
              predict.time = t,
              method       = "KM",
              #KM或 NNE ，NNE加span
              #span = 0.25 * nrow(ovarian)^(-0.20)
              )
}
#估计每半年生存率
survivalROC_data <- data_frame(t = 182.5 * c(1,2,3,4,5,6)) %>%
  mutate(survivalROC = map(t, survivalROC_helper),
         ## Extract scalar AUC
         auc = map_dbl(survivalROC, magrittr::extract2, "AUC"),
         ## Put cut off dependent values in a data_frame
         df_survivalROC = map(survivalROC, function(obj) {
           as_tibble(obj[c("cut.values","TP","FP")])
         })) %>%
  dplyr::select(-survivalROC) %>%
  unnest(cols = c(df_survivalROC)) %>%
  arrange(t, FP, TP)

# 绘制0.5 1 1.5 2年AUC曲线
survivalROC_data %>%
  ggplot(mapping = aes(x = FP, y = TP)) +
  geom_point() +
  geom_line() +
  geom_label(data = survivalROC_data %>% dplyr::select(t,auc) %>% unique,
             mapping = aes(label = sprintf("%.3f", auc)), x = 0.5, y = 0.5) +
  facet_wrap( ~ t) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
        legend.key = element_blank(),
        plot.title = element_text(hjust = 0.5),
        strip.background = element_blank())


四、发散思维

    尽管 C/D 是最常用的应用，但如果研究人员有一个特定的兴趣时间点，以便区分在该时间点发生事件的个体和没有事件的个体，则 I/D 或 I/S 更合适。C/D 和 I/D 通常用于单个marker，而 I/S 可以包括纵向marker。与empire-ROC一样，time-dependent-ROC的X值，可以是marker也可以是预测模型生成的预测概率pre。time-dependent-ROC应用场景：①找到连续变量最佳的cut-off 值；②用于评估变量对生存资料outcome的预测效能；③可以评价某个marker随时间变化对outcome的影响。

参考文献：

  1. Kamarudin AN, Cox T, Kolamunnage-Dona R. Time-dependent ROC curve analysis in medical research: current methods and applications. BMC Med Res Methodol. 2017 Apr 7;17(1):53. doi: 10.1186/s12874-017-0332-6. PMID: 28388943; PMCID: PMC5384160.
  2. Bader Peter,Kreyenberg Hermann,von Stackelberg Arend et al. Monitoring of minimal residual disease after allogeneic stem-cell transplantation in relapsed childhood acute lymphoblastic leukemia allows for the identification of impending relapse: results of the ALL-BFM-SCT 2003 trial.[J] .J Clin Oncol, 2015, 33: 1275-84.
  3. Kahles Florian,Rückbeil Marcia V,Mertens Robert W et al. Glucagon-like peptide 1 levels predict cardiovascular risk in patients with acute myocardial infarction.[J] .Eur Heart J, 2020, 41: 882-889.
  4. Krischer Jeffrey P,Liu Xiang,Vehik Kendra et al. Predicting Islet Cell Autoimmunity and Type 1 Diabetes: An 8-Year TEDDY Study Progress Report.[J] .Diabetes Care, 2019, 42: 1051-1060.

原创不易，转载请说明来自本公众号。