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<- 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)

Logistic回归中C-Statistics计算方法

计算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的可信区间。

发表评论

匿名网友