题记：本章为预测模型系列的第6篇文章，本章主要介绍Logistic回归模型中C-Statistics的三种
计
算方法，请各位斧正~
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"lweight.sav")
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 "低出生体重",1,0)
mydata$race1 "白种人",1,0)
mydata$race2 "黑种人",1,0)
mydata$race3 "其他种族",1,0)
把数据加载到当前工作环境并打包数据
dd'dd')
拟合Logistic回归模型
fit1T,y=T)
方法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 "tpr","fpr")
绘制ROC曲线
plot(perf)
abline(0,1, col = 3, lty = 2)
计算auc即是C-statistics = 0.7382008
auc "auc")
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.7382008
library
