26  C-index的比较

前面介绍了C-index的计算，下面学习C-index的比较。

26.1 二分类资料C-index的比较

二分类资料的AUC和C-index是一样的，所以可以参考Chapter 23关于ROC曲线的显著性检验。

26.2 生存资料C-index的比较

可以使用compareC包，专门用来比较生存资料的C-index。

rm(list = ls())
library(compareC)

还是用之前推文的数据，获取数据可以查看历史推文。

load(file = "./datasets/timeROC.RData")
str(df2)
## 'data.frame':    297 obs. of  8 variables:
##  $ event    : num  0 0 1 0 0 1 0 0 0 0 ...
##  $ age      : int  59 63 65 73 59 66 56 42 61 48 ...
##  $ riskScore: num  -0.249 -0.511 -0.211 -0.427 0.279 ...
##  $ futime   : num  3.03 1.16 1.82 1.52 1.34 ...
##  $ gender   : num  2 2 2 1 2 2 1 2 2 2 ...
##  $ t        : num  4 4 4 3 3 3 5 3 NA 4 ...
##  $ n        : num  1 5 1 1 1 1 3 1 NA 1 ...
##  $ m        : num  1 1 1 1 1 3 1 1 3 3 ...

只要提供4个参数：time，status，第一个指标，第二个指标，即可。

compareC(df2$futime,
         df2$event,
         df2$riskScore,
         df2$age
         )
## $est.c
##       Cxy       Cxz 
## 0.3383690 0.3894508 
## 
## $est.diff_c
## [1] -0.05108181
## 
## $est.vardiff_c
## [1] 0.002124315
## 
## $est.varCxy
## [1] 0.001046384
## 
## $est.varCxz
## [1] 0.00112578
## 
## $est.cov
## [1] 2.392402e-05
## 
## $zscore
## [1] -1.108299
## 
## $pval
## [1] 0.2677329

第1行给出了2个C指数；第2行是2个C指数相减的差值；第3行是2个C指数方差相减的差值；第4、5行是方差；第6行是协方差；第7、8行是z值和p值。

26.3 两个cox模型的比较

下面再多说一点两个cox模型的比较，也是很简单的，方差分析即可，使用anova()函数。

我们用lung数据集进行演示。

library(survival)
library(dplyr)

df1 <- lung %>% 
  mutate(status=ifelse(status == 1,1,0))

建立两个cox模型：

cox_fit1 <- coxph(Surv(time, status) ~ age + sex + ph.ecog + ph.karno + pat.karno,
                  data = df1,x = T, y = T)

cox_fit2 <- coxph(Surv(time, status) ~ ph.ecog + ph.karno + pat.karno,
                  data = df1,x = T, y = T)

直接使用anova()即可：

anova(cox_fit1,cox_fit2)
## Analysis of Deviance Table
##  Cox model: response is  Surv(time, status)
##  Model 1: ~ age + sex + ph.ecog + ph.karno + pat.karno
##  Model 2: ~ ph.ecog + ph.karno + pat.karno
##    loglik  Chisq Df Pr(>|Chi|)  
## 1 -257.97                       
## 2 -261.94 7.9486  2    0.01879 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

rms包也带了一个lrtest()函数，也是直接构建两个cph模型即可进行比较，简单演示一下。

suppressMessages(library(rms))
dd <- datadist(lung)
options(datadist="dd")

cox_fit1 <- cph(Surv(time, status) ~ age + sex + ph.ecog + ph.karno + pat.karno,
                data = df1,x = T, y = T)

cox_fit2 <- cph(Surv(time, status) ~ ph.ecog + ph.karno + pat.karno,
                data = df1,x = T, y = T)

lrtest(cox_fit1,cox_fit2)
## 
## Model 1: Surv(time, status) ~ age + sex + ph.ecog + ph.karno + pat.karno
## Model 2: Surv(time, status) ~ ph.ecog + ph.karno + pat.karno
## 
## L.R. Chisq       d.f.          P 
## 7.94861288 2.00000000 0.01879233

可以看到P值是一样的~