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<- read.spss("lweight.sav")
mydata<- as.data.frame(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 <- ifelse(mydata$low =="低出生体重",1,0) mydata$race1 <- ifelse(mydata$race =="白种人",1,0) mydata$race2 <- ifelse(mydata$race =="黑种人",1,0) mydata$race3 <- ifelse(mydata$race =="其他种族",1,0)
把数据加载到当前工作环境并打包数据
dd<-datadist(mydata) options(datadist='dd')
拟合Logistic回归模型
fit1<-lrm(low~age+ftv+ht+lwt+ptl+smoke+ui+race1+race2,data=mydata,x=T,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<-predict(fit1) #
载入ROCR包
library(ROCR) ## Loading required package: gplots ## Attaching package: 'gplots' ## The following object is masked from 'package:stats': ## ## lowess pred <- prediction(mydata$predvalue, mydata$low) perf<- performance(pred,"tpr","fpr")
绘制ROC曲线
plot(perf) abline(0,1, col = 3, lty = 2)
计算auc即是C-statistics = 0.7382008
auc <- performance(pred,"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(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的可信区间。