ROC after logistic regression

References

Load subarachnoid hemorrhage dataset (see ?aSAH for details)

library(pROC)
data(aSAH)
head(aSAH)

##    gos6 outcome gender age wfns s100b  ndka
## 29    5    Good Female  42    1  0.13  3.01
## 30    5    Good Female  37    1  0.14  8.54
## 31    5    Good Female  42    1  0.10  8.09
## 32    5    Good Female  27    1  0.04 10.42
## 33    1    Poor Female  42    3  0.13 17.40
## 34    1    Poor   Male  48    2  0.10 12.75

Fit multiple logistic regression

## Logistic regression
fitLogistic <- glm(formula = outcome ~ gender + age + s100b + ndka,
                   family  = binomial(link = "logit"),
                   data    = aSAH)
summary(fitLogistic)

## 
## Call:
## glm(formula = outcome ~ gender + age + s100b + ndka, family = binomial(link = "logit"), 
##     data = aSAH)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.974  -0.753  -0.504   0.707   2.239  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -3.3322     1.0110   -3.30  0.00098 ***
## genderFemale  -1.1025     0.5049   -2.18  0.02900 *  
## age            0.0335     0.0183    1.83  0.06663 .  
## s100b          4.9955     1.2758    3.92  0.00009 ***
## ndka           0.0290     0.0168    1.73  0.08374 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 148.04  on 112  degrees of freedom
## Residual deviance: 113.08  on 108  degrees of freedom
## AIC: 123.1
## 
## Number of Fisher Scoring iterations: 5


## Easier view
library(epicalc)
logistic.display(fitLogistic)

## 
## Logistic regression predicting outcome : Poor vs Good 
##  
##                        crude OR(95%CI)        adj. OR(95%CI)          P(Wald's test) P(LR-test)
## gender: Female vs Male 0.46 (0.21,1.02)       0.33 (0.12,0.89)        0.029          0.026     
##                                                                                                
## age (cont. var.)       1.03 (1,1.06)          1.03 (1,1.07)           0.067          0.061     
##                                                                                                
## s100b (cont. var.)     134.87 (12.9,1409.56)  147.75 (12.12,1800.75)  < 0.001        < 0.001   
##                                                                                                
## ndka (cont. var.)      1.02 (0.99,1.05)       1.03 (1,1.06)           0.084          0.076     
##                                                                                                
## Log-likelihood = -56.5391
## No. of observations = 113
## AIC value = 123.0783

AUC from epicalc::lroc()

resLroc <- lroc(fitLogistic)

plot of chunk unnamed-chunk-4 

resLroc$auc

## [1] 0.8096

AUC from pROC::roc()

resRoc <- roc(aSAH$outcome ~ fitLogistic$fitted)
plot(resRoc, legacy.axes = TRUE)

plot of chunk unnamed-chunk-5 

## 
## Call:
## roc.formula(formula = aSAH$outcome ~ fitLogistic$fitted)
## 
## Data: fitLogistic$fitted in 72 controls (aSAH$outcome Good) < 41 cases (aSAH$outcome Poor).
## Area under the curve: 0.81

AUC (C in Rank Discrim. Indexes) from rms::lrm()

library(rms)
resLrm <- lrm(formula = outcome ~ gender + age + s100b + ndka,
              data    = aSAH,
              x = TRUE, y = TRUE)
resLrm

## 
## Logistic Regression Model
## 
## lrm(formula = outcome ~ gender + age + s100b + ndka, data = aSAH, 
##     x = TRUE, y = TRUE)
## 
##                          Model Likelihood     Discrimination    Rank Discrim.    
##                             Ratio Test            Indexes          Indexes       
## Obs              113    LR chi2      34.96    R2       0.364    C       0.809    
##  Good             72    d.f.             4    g        1.844    Dxy     0.619    
##  Poor             41    Pr(> chi2) <0.0001    gr       6.320    gamma   0.619    
## max |deriv| 0.000003                          gp       0.285    tau-a   0.289    
##                                               Brier    0.165                     
## 
##               Coef    S.E.   Wald Z Pr(>|Z|)
## Intercept     -3.3322 1.0111 -3.30  0.0010  
## gender=Female -1.1025 0.5049 -2.18  0.0290  
## age            0.0335 0.0183  1.83  0.0666  
## s100b          4.9955 1.2758  3.92  <0.0001 
## ndka           0.0290 0.0168  1.73  0.0838