Logistic回归中C-Statistics计算方法

来源：临床研究与医学统计评论5,555

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包

  1. library(foreign) 
  2. library(rms)
  3. ## Loading required package: Hmisc
  4. ## Loading required package: lattice
  5. ## Loading required package: survival
  6. ## Loading required package: Formula
  7. ## Loading required package: ggplot2
  8. ## Attaching package: 'Hmisc'
  9. ## The following objects are masked from 'package:base':
  10. ## 
  11. ##     format.pval, units
  12. ## Loading required package: SparseM
  13. ## Attaching package: 'SparseM'
  14. ## The following object is masked from 'package:base':
  15. ## 
  16. ##     backsolve

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

  1. mydata<- read.spss("lweight.sav") 
  2. mydata<- as.data.frame(mydata)  
  3. head(mydata) 
  4. ##   id      low age lwt     race  smoke ptl       ht ui ftv  bwt
  5. ## 1 85 正常体重  19 182   黑种人 不吸烟   0 无妊高症 有   0 2523
  6. ## 2 86 正常体重  33 155 其他种族 不吸烟   0 无妊高症 无   3 2551
  7. ## 3 87 正常体重  20 105   白种人   吸烟   0 无妊高症 无   1 2557
  8. ## 4 88 正常体重  21 108   白种人   吸烟   0 无妊高症 有   2 2594
  9. ## 5 89 正常体重  18 107   白种人   吸烟   0 无妊高症 有   0 2600
  10. ## 6 91 正常体重  21 124 其他种族 不吸烟   0 无妊高症 无   0 2622

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

  1. mydata$low <- ifelse(mydata$low =="低出生体重",1,0) 
  2. mydata$race1 <- ifelse(mydata$race =="白种人",1,0)
  3. mydata$race2 <- ifelse(mydata$race =="黑种人",1,0)
  4. mydata$race3 <- ifelse(mydata$race =="其他种族",1,0)

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

  1. dd<-datadist(mydata)
  2. options(datadist='dd')

拟合Logistic回归模型

  1. 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。

  1. fit1 
  2. ## Logistic Regression Model
  3. ##  
  4. ##  lrm(formula = low ~ age + ftv + ht + lwt + ptl + smoke + ui + 
  5. ##      race1 + race2, data = mydata, x = T, y = T)
  6. ##  
  7. ##                       Model Likelihood     Discrimination    Rank Discrim.    
  8. ##                          Ratio Test           Indexes           Indexes       
  9. ##  Obs           189    LR chi2     31.12    R2       0.213    C       0.738    
  10. ##   0            130    d.f.            9    g        1.122    Dxy     0.476    
  11. ##   1             59    Pr(> chi2) 0.0003    gr       3.070    gamma   0.477    
  12. ##  max |deriv| 7e-05                         gp       0.207    tau-a   0.206    
  13. ##                                            Brier    0.181                     
  14. ##  
  15. ##             Coef    S.E.   Wald Z Pr(>|Z|)
  16. ##  Intercept   1.1427 1.0873  1.05  0.2933  
  17. ##  age        -0.0255 0.0366 -0.69  0.4871  
  18. ##  ftv         0.0321 0.1708  0.19  0.8509  
  19. ##  ht=妊高症   1.7631 0.6894  2.56  0.0105  
  20. ##  lwt        -0.0137 0.0068 -2.02  0.0431  
  21. ##  ptl         0.5517 0.3446  1.60  0.1094  
  22. ##  smoke=吸烟  0.9275 0.3986  2.33  0.0200  
  23. ##  ui=有       0.6488 0.4676  1.39  0.1653  
  24. ##  race1      -0.9082 0.4367 -2.08  0.0375  
  25. ##  race2       0.3293 0.5339  0.62  0.5374  

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

计算模型预测概率

  1. mydata$predvalue<-predict(fit1) #

载入ROCR包

  1. library(ROCR) 
  2. ## Loading required package: gplots
  3. ## Attaching package: 'gplots'
  4. ## The following object is masked from 'package:stats':
  5. ## 
  6. ##     lowess
  7. pred <- prediction(mydata$predvalue, mydata$low) 
  8. perf<- performance(pred,"tpr","fpr")

绘制ROC曲线

  1. plot(perf)
  2. abline(0,1, col = 3, lty = 2)

Logistic回归中C-Statistics计算方法

计算auc即是C-statistics = 0.7382008

  1. auc <- performance(pred,"auc")
  2. auc
  3. ## An object of class "performance"
  4. ## Slot "x.name":
  5. ## [1] "None"
  6. ## 
  7. ## Slot "y.name":
  8. ## [1] "Area under the ROC curve"
  9. ## 
  10. ## Slot "alpha.name":
  11. ## [1] "none"
  12. ## 
  13. ## Slot "x.values":
  14. ## list()
  15. ## 
  16. ## Slot "y.values":
  17. ## [[1]]
  18. ## [1] 0.7382008
  19. ## 
  20. ## Slot "alpha.values":
  21. ## list()

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

  1. library(Hmisc)
  2. somers2(mydata$predvalue, mydata$low)
  3. ##           C         Dxy           n     Missing 
  4. ##   0.7382008   0.4764016 189.0000000   0.0000000

4. 总结与讨论

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

发表评论

匿名网友

拖动滑块以完成验证
加载失败