R中进行生存资料AUC的计算-经管之家官网！

除了Logistic回归能够计算预测值与真实值的综合一致程度，也即AUC外，生存分析同样能计算AUC。但生存分析假设不同，生存函数的构建也就不尽相同。因而，不同方法计算得到的生存资料的AUC也不尽相同。同时，生存函数ROC曲线的绘制也更加复杂。下面提供几种常用的基于R计算生存资料AUC的方法，并提供两种时间依赖的生存函数的AUC计算示例。
但目前来说，还没有方法可以直接比较两个生存函数的AUC是否存在显著的统计学差异。
如果对R中进行生存资料分析的基础程序包不是很熟悉的话，可以先看看之前的一个帖子：
https://bbs.pinggu.org/thread-4127756-1-1.html

1. install.packages(c("clinfun","CPE","risksetROC", "timeROC",
   
2. "survivalROC","survC1","survIDINRI"))
   
3. 
   
4. library(survival)
   
5. 
   
6. set.seed(666)
   
7. age <- rnorm(400, 50, 10)
   
8. bp<- rnorm(400,120, 15)
   
9. d.time <- rexp(400)
   
10. cens <- runif(400,.5,2)
    
11. death<- d.time <= cens
    
12. d.time <- pmin(d.time, cens)
    
13. 
    
14. fit <- coxph(Surv(d.time,death) ~ age + bp)
    
15. summary(fit)
    
16. # Concordance = 0.502(se = 0.019 )
    
17. 
    
18. # Compute the concordance between a right-censored
    
19. # survival time and a single continuous covariate
    
20. survConcordance(Surv(d.time,death) ~ predict(fit))
    
21. # Concordance = 0.502101(se = 0.01883032 )
    
22. # Concordance = concordant/(concordant+discordant)
    
23. 
    
24. library(Hmisc)
    
25. # Computes the c index and the corresponding generalization of
    
26. # Somers' Dxy rank correlation for a censored response variable.
    
27. rcorrcens(Surv(d.time,death) ~ predict(fit))
    
28. # C index: 1 - 0.498 = 0.502
    
29. # C index = (1+aDxy)/2 = (1+0.004)/2 = 0.502
    
30. 
    
31. library(risksetROC)
    
32. tmax <- max(d.time)
    
33. # Create and plot AUC based on incident/dynamic definition of Heagerty
    
34. AUC <- risksetAUC(Stime=d.time, status=death,
    
35. marker=age + bp, method="Cox", tmax=tmax)
    
36. AUC$Cindex
    
37. # 0.5020833
    
38. 
    
39. library(clinfun)
    
40. # Calculate Gonen & Heller concordance probability estimate (CPE) for
    
41. # for the Cox proportional hazards model.
    
42. coxphCPE(fit)
    
43. # CPE= 0.50884256, se.CPE = 0.01678487
    
44. 
    
45. library(CPE)
    
46. # Calculate Gonen & Heller concordance probability estimate (CPE)
    
47. # for the Cox proportional hazards model.
    
48. phcpe(fit, CPE.SE=TRUE,out.ties=TRUE)
    
49. # CPE= 0.5088426, se.CPE = 0.01678487
    
50. 
    
51. 
    
52. library(survivalROC)
    
53. # Time-dependent ROC curve estimation from censored survival data
    
54. AUC2 <- survivalROC(Stime=d.time, status=death,
    
55. marker = eta,
    
56. predict.time =265, method="KM")
    
57. AUC2$AUC
    
58. # 0.6124465
    
59. plot(AUC2$FP,
    
60. AUC2$TP,
    
61. type="l",
    
62. xlim=c(0,1), ylim=c(0,1),
    
63. xlab=paste( "FP", "\n", "AUC = ",round(AUC2$AUC,3)),
    
64. ylab="TP",main="AUC2, Method = KM \n Year = 1")
    
65. abline(0,1)
    
66. 
    
67. library(timeROC)
    
68. # Time-dependent ROC curve estimation
    
69. ROC<-timeROC(T=d.time,
    
70. delta=death,
    
71. marker=eta,
    
72. other_markers=as.matrix(bp),
    
73. cause=1,
    
74. weighting="marginal",
    
75. times=quantile(d.time,probs=seq(0.2,0.8,0.1)),
    
76. ROC = TRUE,
    
77. iid = TRUE)
    
78. ROC
    
79. confint(ROC)

复制代码