1. 背景知识

上一节中我们讲解了R语言在Cox回归模型中计算C-index的方法，参见：预测模型系列 05 -- Cox回归中C-index的两种常用计算方法。本节我们将介绍用R语言计算Logistic回归中C-Statistics。上一节我们提到在Logistic回归模型中ROC曲线下面积=C-Statistics，所以SPSS软件也是可以计算C-Statistics，本文不再赘述。下面我们将以一个经典的Logistic回归的案例讲解R语言计算C-statistics的方法。

2. 案例分析

Hosmer和 Lemeshow于1989年研究了低出生体重婴儿的影响因素。结果变量为是否娩出低出生体重儿(变量名为“low”，二分类变量，1=低出生体重，即婴儿出生体重<2500g;0=非低出生体重)，考虑的影响因素(自变量)有：产妇妊娠前体重(lwt，磅)；产妇年龄(age，岁)；产妇在妊娠期间是否吸烟(smoke，0＝未吸、1＝吸烟)；本次妊娠前早产次数(ptl，次)；是否患有高血压(ht，0＝未患、1＝患病)；子宫对按摩、催产素等刺激引起收缩的应激性(ui，0＝无、1＝有)；妊娠前三个月社区医生随访次数(ftv，次)；种族(race，1＝白人、2＝黑人、3＝其他民族)。

本案例因变量是二分类变量(是否低出生体重儿)，研究目的是探讨低出生体重儿的独立影响因素，符合二元Logistic回归的应用条件。我们构建一个以“age+ftv+ht+lwt+ptl+smoke+ui+race”为自变量，以“low”为因变量的Logistic回归方程。基于此Logistic回归模型我们有三种方法可以计算其C-Statistics。

方法1. 利用 {rms} 包中的 lrm 函数构建Logistic回归模型，直接读取模型Rank Discrim.参数 C，即为C-Statistics。

方法2. 构建Logistic回归模型，predict函数计算模型预测概率，然后利用ROCR包根据此预测概率画ROC曲线，并计算曲线下面积AUC，此即为C-Statistics。注：此方法与SPSS中的计算方法一致。

方法3. 构建Logistic回归模型，predict函数计算模型预测概率，利用Hmisc包中somers2函数直接计算ROC曲线下面积AUC。注：此方法与SPSS中的计算方法一致。

3. R计算过程

载入foreign包与rms包library(foreign)

library(rms)

## Loading required package: Hmisc

## Loading required package: lattice

## Loading required package: survival

## Loading required package: Formula

## Loading required package: ggplot2

## Attaching package: 'Hmisc'

## The following objects are masked from 'package:base':

##

##     format.pval, units

## Loading required package: SparseM

## Attaching package: 'SparseM'

## The following object is masked from 'package:base':

##

##     backsolve

导入sav格式的数据，并把导入数据转化为数据框结构，展示数据的前6行mydata

mydata

head(mydata)

##   id      low age lwt     race  smoke ptl       ht ui ftv  bwt

## 1 85 正常体重  19 182   黑种人 不吸烟   0 无妊高症 有   0 2523

## 2 86 正常体重  33 155 其他种族 不吸烟   0 无妊高症 无   3 2551

## 3 87 正常体重  20 105   白种人   吸烟   0 无妊高症 无   1 2557

## 4 88 正常体重  21 108   白种人   吸烟   0 无妊高症 有   2 2594

## 5 89 正常体重  18 107   白种人   吸烟   0 无妊高症 有   0 2600

## 6 91 正常体重  21 124 其他种族 不吸烟   0 无妊高症 无   0 2622

设置结局变量为二分类，把race变量设置为哑变量mydata$low

mydata$race1

mydata$race2

mydata$race3

把数据加载到当前工作环境并打包数据dd

options(datadist='dd')

拟合Logistic回归模型fit1

方法1.直接读取模型参数中Rank Discrim.参数C，即是C-Statistics = 0.738。fit1

## Logistic Regression Model

##

##  lrm(formula = low ~ age + ftv + ht + lwt + ptl + smoke + ui +

##      race1 + race2, data = mydata, x = T, y = T)

##

##                       Model Likelihood     Discrimination    Rank Discrim.

##                          Ratio Test           Indexes           Indexes

##  Obs           189    LR chi2     31.12    R2       0.213    C       0.738

##   0            130    d.f.            9    g        1.122    Dxy     0.476

##   1             59    Pr(> chi2) 0.0003    gr       3.070    gamma   0.477

##  max |deriv| 7e-05                         gp       0.207    tau-a   0.206

##                                            Brier    0.181

##

##             Coef    S.E.   Wald Z Pr(>|Z|)

##  Intercept   1.1427 1.0873  1.05  0.2933

##  age        -0.0255 0.0366 -0.69  0.4871

##  ftv         0.0321 0.1708  0.19  0.8509

##  ht=妊高症   1.7631 0.6894  2.56  0.0105

##  lwt        -0.0137 0.0068 -2.02  0.0431

##  ptl         0.5517 0.3446  1.60  0.1094

##  smoke=吸烟  0.9275 0.3986  2.33  0.0200

##  ui=有       0.6488 0.4676  1.39  0.1653

##  race1      -0.9082 0.4367 -2.08  0.0375

##  race2       0.3293 0.5339  0.62  0.5374

方法2. ROCR包计算AUC，代码如下：

计算模型预测概率mydata$predvalue

载入ROCR包library(ROCR)

## Loading required package: gplots

## Attaching package: 'gplots'

## The following object is masked from 'package:stats':

##

##     lowess

pred

perf

绘制ROC曲线plot(perf)

abline(0,1, col = 3, lty = 2)

计算auc即是C-statistics = 0.7382008auc

auc

## An object of class "performance"

## Slot "x.name":

## [1] "None"

##

## Slot "y.name":

## [1] "Area under the ROC curve"

##

## Slot "alpha.name":

## [1] "none"

##

## Slot "x.values":

## list()

##

## Slot "y.values":

## [[1]]

## [1] 0.7382008

##

## Slot "alpha.values":

## list()

Hmisc包somers2() 函数计算AUC = 0.7382008library(Hmisc)

somers2(mydata$predvalue, mydata$low)

##           C         Dxy           n     Missing

##   0.7382008   0.4764016 189.0000000   0.0000000

4. 总结与讨论

至此本章关于Logistic回归中C-Statistics的计算三种方法演示完毕。不管那种方法都未给出标准误，所以可信区间的计算就很麻烦，如果一定要报告C-Statistics可信区间，可以考虑使用SPSS软件进行ROC分析，软件可以直接给出AUC的可信区间。