25 C-index的计算

C-statistic（C-index，C指数）是评价模型区分度的指标之一，在逻辑回归模型中，C-statistic就是AUC，在生存资料中，C-statistic和AUC略有不同。

下面给大家分别介绍logistic和cox回归的C-statistic计算方法。

25.1 logistic回归

今天学习C-index的4种计算方法，在二分类变量中，C-statistic就是AUC，二者在数值上是一样的。

使用lowbirth数据集，这个数据集是关于低出生体重儿是否会死亡的数据集，其中dead这一列是结果变量，0代表死亡，1代表存活，其余列都是预测变量（获取lowbirth数据请在公众号：医学和生信笔记后台回复20220520，或者在粉丝qq群文件自取）。

我们首先对这个数据做一些预处理（参考：Chapter 7）

rm(list = ls())
lowbirth <- read.csv("./datasets/lowbirth.csv")

library(dplyr)

tmp <- lowbirth %>% 
  mutate(across(where(is.character),as.factor),
         vent = factor(vent),
         #dead = factor(dead), 下面介绍的方法2不能是因子型，所以不转换了
         race = case_when(race %in% 
                            c("native American","oriental") ~ "other",
                          .default = race),
         race = factor(race))

str(tmp)
## 'data.frame':    565 obs. of  10 variables:
##  $ birth   : num  81.5 81.6 81.6 81.6 81.6 ...
##  $ lowph   : num  7.25 7.06 7.25 6.97 7.32 ...
##  $ pltct   : int  244 114 182 54 282 153 229 182 361 378 ...
##  $ race    : Factor w/ 3 levels "black","other",..: 3 1 1 1 1 1 3 1 3 3 ...
##  $ bwt     : int  1370 620 1480 925 1255 1350 1310 1110 1180 970 ...
##  $ delivery: Factor w/ 2 levels "abdominal","vaginal": 1 2 2 1 2 1 2 2 1 2 ...
##  $ apg1    : int  7 1 8 5 9 4 6 6 6 2 ...
##  $ vent    : Factor w/ 2 levels "0","1": 1 2 1 2 1 1 2 1 1 2 ...
##  $ sex     : Factor w/ 2 levels "female","male": 1 1 2 1 1 1 2 2 2 1 ...
##  $ dead    : int  0 1 0 1 0 0 0 0 0 1 ...

25.1.1 方法1：rms

使用rms包构建模型，模型结果中Rank Discrim.下面的C就是C-Statistic，本模型中C-Statistic = 0.879。

library(rms)

dd <- datadist(tmp)
options(datadist="dd")

fit2 <- lrm(dead ~ birth + lowph + pltct + bwt + vent + race,
            data = tmp)

fit2
## Logistic Regression Model
## 
## lrm(formula = dead ~ birth + lowph + pltct + bwt + vent + race, 
##     data = tmp)
## 
##                        Model Likelihood      Discrimination    Rank Discrim.    
##                              Ratio Test             Indexes          Indexes    
## Obs           565    LR chi2     167.56      R2       0.432    C       0.879    
##  0            471    d.f.             7      R2(7,565)0.247    Dxy     0.759    
##  1             94    Pr(> chi2) <0.0001    R2(7,235.1)0.495    gamma   0.759    
## max |deriv| 1e-06                            Brier    0.095    tau-a   0.211    
## 
##            Coef    S.E.    Wald Z Pr(>|Z|)
## Intercept  39.5788 11.0070  3.60  0.0003  
## birth      -0.1201  0.0914 -1.31  0.1890  
## lowph      -4.1451  1.1881 -3.49  0.0005  
## pltct      -0.0017  0.0019 -0.91  0.3644  
## bwt        -0.0031  0.0006 -5.14  <0.0001 
## vent=1      2.7526  0.7436  3.70  0.0002  
## race=other -1.1974  0.8448 -1.42  0.1564  
## race=white -0.3377  0.2953 -1.14  0.2529

25.1.2 方法2：Hmisc

使用Hmisc包。结果中的C就是C-Statistic。

# 先计算线性预测值
tmp$predvalue<-predict(fit2)

library(Hmisc)
somers2(tmp$predvalue, tmp$dead)
##           C         Dxy           n     Missing 
##   0.8793875   0.7587749 565.0000000   0.0000000

# 或者用
rcorr.cens(tmp$predvalue, tmp$dead)
##        C Index            Dxy           S.D.              n        missing 
##   8.793875e-01   7.587749e-01   3.417295e-02   5.650000e+02   0.000000e+00 
##     uncensored Relevant Pairs     Concordant      Uncertain 
##   5.650000e+02   8.854800e+04   7.786800e+04   0.000000e+00

25.1.3 方法3：ROCR

ROCR包计算AUC，logistic回归的AUC就是C-statistic。这种方法和SPSS得到的一样。

library(ROCR)

tmp$predvalue<-predict(fit2)

# 取出C-Statistics，和上面结果一样
pred <- prediction(tmp$predvalue, tmp$dead)

auc <- round(performance(pred, "auc")@y.values[[1]],digits = 4)
auc
## [1] 0.8794

这个包也是用来画ROC曲线常用的包，可以根据上面的结果直接画出ROC曲线，这里就不重复演示了。

25.1.4 方法4：pROC

pROC包计算AUC，这个包也是画ROC曲线常用的R包，但是这个包在使用时需要注意，这部分内容会在后面详细介绍。

library(pROC)

# 计算AUC，也就是C-statistic
roc(tmp$dead, tmp$predvalue, legacy.axes = T, print.auc = T, print.auc.y = 45)
## 
## Call:
## roc.default(response = tmp$dead, predictor = tmp$predvalue, legacy.axes = T,     print.auc = T, print.auc.y = 45)
## 
## Data: tmp$predvalue in 471 controls (tmp$dead 0) < 94 cases (tmp$dead 1).
## Area under the curve: 0.8794

也是可以直接画法ROC曲线的，请参考前面的文章。

25.2 cox回归

cox回归的C-statistic可以用survival包计算，需要注意，生存分析的C-statistic和AUC是不一样的，生存分析的C指数一般是指Harrell’s C指数（Harrell就是rms包的作者，统计学大佬）。

使用survival包自带的lung数据集进行演示。

rm(list = ls())
library(survival)
library(dplyr)

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

str(df1)
## 'data.frame':    228 obs. of  10 variables:
##  $ inst     : num  3 3 3 5 1 12 7 11 1 7 ...
##  $ time     : num  306 455 1010 210 883 ...
##  $ status   : num  0 0 1 0 0 1 0 0 0 0 ...
##  $ age      : num  74 68 56 57 60 74 68 71 53 61 ...
##  $ sex      : num  1 1 1 1 1 1 2 2 1 1 ...
##  $ ph.ecog  : num  1 0 0 1 0 1 2 2 1 2 ...
##  $ ph.karno : num  90 90 90 90 100 50 70 60 70 70 ...
##  $ pat.karno: num  100 90 90 60 90 80 60 80 80 70 ...
##  $ meal.cal : num  1175 1225 NA 1150 NA ...
##  $ wt.loss  : num  NA 15 15 11 0 0 10 1 16 34 ...

R语言自带的coxph函数即可给出C-index，非常简单：

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

summary(cox_fit1)
## Call:
## coxph(formula = Surv(time, status) ~ age + sex + ph.ecog, data = df1, 
##     x = T, y = T)
## 
##   n= 227, number of events= 63 
##    (1 observation deleted due to missingness)
## 
##             coef exp(coef) se(coef)      z Pr(>|z|)  
## age     -0.02193   0.97831  0.01496 -1.465   0.1428  
## sex      0.63989   1.89626  0.26588  2.407   0.0161 *
## ph.ecog -0.19941   0.81922  0.20708 -0.963   0.3356  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##         exp(coef) exp(-coef) lower .95 upper .95
## age        0.9783     1.0222    0.9500     1.007
## sex        1.8963     0.5274    1.1261     3.193
## ph.ecog    0.8192     1.2207    0.5459     1.229
## 
## Concordance= 0.584  (se = 0.042 )
## Likelihood ratio test= 10.28  on 3 df,   p=0.02
## Wald test            = 10.14  on 3 df,   p=0.02
## Score (logrank) test = 10.43  on 3 df,   p=0.02

Concordance就是C-statistic，本次示例中为0.584。

也可以用以下代码提取结果：

summary(cox_fit1)$concordance
##          C      se(C) 
## 0.58390352 0.04179314

以上就是C-statistic的计算。