除了Logistic回歸能夠計(jì)算預(yù)測值與真實(shí)值的綜合一致程度，也即AUC外，生存分析同樣能計(jì)算AUC。但生存分析假設(shè)不同，生存函數(shù)的構(gòu)建也就不盡相同。因而，不同方法計(jì)算得到的生存資料的AUC也不盡相同。同時(shí)，生存函數(shù)ROC曲線的繪制也更加復(fù)雜。下面提供幾種常用的基于R計(jì)算生存資料AUC的方法，并提供兩種時(shí)間依賴的生存函數(shù)的AUC計(jì)算示例。

但目前來說，還沒有方法可以直接比較兩個(gè)生存函數(shù)的AUC是否存在顯著的統(tǒng)計(jì)學(xué)差異。

install.packages(c("clinfun","CPE","risksetROC", "timeROC",

? ?? ?? ?? ?? ?? ? "survivalROC","survC1","survIDINRI"))



library(survival)



set.seed(666)

age <- rnorm(400, 50, 10)

bp??<- rnorm(400,120, 15)

d.time <- rexp(400)

cens? ?<- runif(400,.5,2)

death??<- d.time <= cens

d.time <- pmin(d.time, cens)



fit <- coxph(Surv(d.time,death) ~ age + bp)

summary(fit) 

# Concordance = 0.502??(se = 0.019 ) 



# Compute the concordance between a right-censored 

# survival time and a single continuous covariate

survConcordance(Surv(d.time,death) ~ predict(fit))

# Concordance = 0.502101??(se = 0.01883032 )

# Concordance = concordant/(concordant+discordant)



library(Hmisc)

# Computes the c index and the corresponding generalization of 

# Somers' Dxy rank correlation for a censored response variable.

rcorrcens(Surv(d.time,death) ~ predict(fit))

# C index: 1 - 0.498 = 0.502 

# C index = (1+aDxy)/2 = (1+0.004)/2 = 0.502



library(risksetROC)

tmax <- max(d.time)

# Create and plot AUC based on incident/dynamic definition of Heagerty

AUC <- risksetAUC(Stime=d.time, status=death,

? ?? ?? ?? ?? ?? ?marker=age + bp, method="Cox", tmax=tmax)

AUC$Cindex

# 0.5020833



library(clinfun)

# Calculate Gonen & Heller concordance probability estimate (CPE) for 

# for the Cox proportional hazards model.

coxphCPE(fit)

# CPE= 0.50884256, se.CPE = 0.01678487 



library(CPE)

# Calculate Gonen & Heller concordance probability estimate (CPE) 

# for the Cox proportional hazards model.

phcpe(fit, CPE.SE=TRUE,out.ties=TRUE)

# CPE= 0.5088426, se.CPE = 0.01678487 





library(survivalROC)

# Time-dependent ROC curve estimation from censored survival data

AUC2 <- survivalROC(Stime=d.time, status=death,? ???

? ?? ?? ?? ?? ?? ???marker = eta,? ???

? ?? ?? ?? ?? ?? ???predict.time =??265, method="KM")

AUC2$AUC

# 0.6124465 

plot(AUC2$FP, 

? ???AUC2$TP, 

? ???type="l", 

? ???xlim=c(0,1), ylim=c(0,1),? ?

? ???xlab=paste( "FP", "\n", "AUC = ",round(AUC2$AUC,3)), 

? ???ylab="TP",main="AUC2, Method = KM \n Year = 1")

abline(0,1)



library(timeROC)

# Time-dependent ROC curve estimation

ROC<-timeROC(T=d.time,

? ?? ?? ?? ? delta=death,

? ?? ?? ?? ? marker=eta,

? ?? ?? ?? ? other_markers=as.matrix(bp),

? ?? ?? ?? ? cause=1,

? ?? ?? ?? ? weighting="marginal",

? ?? ?? ?? ? times=quantile(d.time,probs=seq(0.2,0.8,0.1)),

? ?? ?? ?? ? ROC = TRUE,

? ?? ?? ?? ? iid = TRUE)

ROC

confint(ROC)