12
$\begingroup$

I performed multivariate logistic regression with the dependent variable Y being death at a nursing home within a certain period of entry and got the following results (note if the variables starts in A it is a continuous value while those starting in B are categorical):

Call:
glm(Y ~ A1 + B2 + B3 + B4 + B5 + A6 + A7 + A8 + A9, data=mydata, family=binomial)
Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.0728  -0.2167  -0.1588  -0.1193   3.7788  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  20.048631  6.036637   3.321 0.000896 ***
A1           0.051167   0.016942   3.020 0.002527 ** 
B2          -0.664940   0.304299  -2.185 0.028878 *  
B3          -2.825281   0.633072  -4.463 8.09e-06 ***
B4          -2.547931   0.957784  -2.660 0.007809 ** 
B5          -2.862460   1.385118  -2.067 0.038774 *  
A6          -0.129808   0.041286  -3.144 0.001666 ** 
A7           0.020016   0.009456   2.117 0.034276 *  
A8          -0.707924   0.253396  -2.794 0.005210 ** 
A9           0.003453   0.001549   2.229 0.025837 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 485.10  on 2206  degrees of freedom
Residual deviance: 417.28  on 2197  degrees of freedom
AIC: 437.28

Number of Fisher Scoring iterations: 7

 (Intercept)           A1           B2           B3           B4           B5           A6           A7           A8           A9 
5.093426e+08 1.052499e+00 5.143045e-01 5.929197e-02 7.824340e-02 5.712806e-02 8.782641e-01 1.020218e+00 4.926657e-01 1.003459e+00 

                   2.5 %       97.5 %
(Intercept) 3.703525e+03 7.004944e+13
A1          1.018123e+00 1.088035e+00
B2          2.832698e-01 9.337710e-01
B3          1.714448e-02 2.050537e-01
B4          1.197238e-02 5.113460e-01
B5          3.782990e-03 8.627079e-01
A6          8.099945e-01 9.522876e-01
A7          1.001484e+00 1.039302e+00
A8          2.998207e-01 8.095488e-01
A9          1.000416e+00 1.006510e+00

As you can see, all of the variables are "significant" in that their p values are below the usual threshold of 0.05. However looking at the coefficients, I'm not quite sure what to make of these results. It seems that although these variables contribute to the model, looking at the odds ratios, they don't seem to really seem to have much predictive power. Of note, when I calculated the AUC, I got approximately 0.8.

Can I say that this model is better at predicting against mortality (e.g. predicting that seniors will live past the prescribed period) compared to predicting for mortality?

$\endgroup$
2
  • 4
    $\begingroup$ This is not a “multivariate” model as you only have one outcome/dependent variable. Fitting a model with several predictors is called “multiple” regression. $\endgroup$ – Gala Jul 19 '13 at 6:41
  • $\begingroup$ I really appreciate the comments and answers - I will definitely take a peek at the rms package. However my question still lingers- most of the odds ratios are <1, would that imply that those variables in this model are better at predicting against the outcome? $\endgroup$ – oort Jul 20 '13 at 11:37
22
$\begingroup$

I would suggest that you use Frank Harrell's excellent rms package. It contains many useful functions to validate and calibrate your model. As far as I know, you cannot assess predictive performance solely based on the coefficients. Further, I would suggest that you use the bootstrap to validate the model. The AUC or concordance-index (c-index) is a useful measure of predictive performance. A c-index of $0.8$ is quite high but as in many predictive models, the fit of your model is likely overoptimistic (overfitting). This overoptimism can be assessed using bootstrap. But let me give an example:

#-----------------------------------------------------------------------------
# Load packages
#-----------------------------------------------------------------------------

library(rms)

#-----------------------------------------------------------------------------
# Load data
#-----------------------------------------------------------------------------

mydata <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")

mydata$rank <- factor(mydata$rank)

#-----------------------------------------------------------------------------
# Fit logistic regression model
#-----------------------------------------------------------------------------

mylogit <- lrm(admit ~ gre + gpa + rank, x=TRUE, y=TRUE, data = mydata)
mylogit

                      Model Likelihood     Discrimination    Rank Discrim.    
                         Ratio Test            Indexes          Indexes       
Obs           400    LR chi2      41.46    R2       0.138    C       0.693    
 0            273    d.f.             5    g        0.838    Dxy     0.386    
 1            127    Pr(> chi2) <0.0001    gr       2.311    gamma   0.387    
max |deriv| 2e-06                          gp       0.167    tau-a   0.168    
                                           Brier    0.195                     

          Coef    S.E.   Wald Z Pr(>|Z|)
Intercept -3.9900 1.1400 -3.50  0.0005  
gre        0.0023 0.0011  2.07  0.0385  
gpa        0.8040 0.3318  2.42  0.0154  
rank=2    -0.6754 0.3165 -2.13  0.0328  
rank=3    -1.3402 0.3453 -3.88  0.0001  
rank=4    -1.5515 0.4178 -3.71  0.0002 

On the bottom you see the usual regression coefficients with corresponding $p$-values. On the top right, you see several discrimination indices. The C denotes the c-index (AUC), and a c-index of $0.5$ denotes random splitting whereas a c-index of $1$ denotes perfect prediction. Dxy is Somers' $D_{xy}$ rank correlation between the predicted probabilities and the observed responses. $D_{xy}$ has simple relationship with the c-index: $D_{xy}=2(c-0.5)$. A $D_{xy}$ of $0$ occurs when the model's predictions are random and when $D_{xy}=1$, the model is perfectly discriminating. In this case, the c-index is $0.693$ which is slightly better than chance but a c-index of $>0.8$ is good enough for predicting the outcomes of individuals.

As said above, the model is likely overoptimistic. We now use bootstrap to quantify the optimism:

#-----------------------------------------------------------------------------
# Validate model using bootstrap
#-----------------------------------------------------------------------------

my.valid <- validate(mylogit, method="boot", B=1000)
my.valid

          index.orig training    test optimism index.corrected    n
Dxy           0.3857   0.4033  0.3674   0.0358          0.3498 1000
R2            0.1380   0.1554  0.1264   0.0290          0.1090 1000
Intercept     0.0000   0.0000 -0.0629   0.0629         -0.0629 1000
Slope         1.0000   1.0000  0.9034   0.0966          0.9034 1000
Emax          0.0000   0.0000  0.0334   0.0334          0.0334 1000
D             0.1011   0.1154  0.0920   0.0234          0.0778 1000
U            -0.0050  -0.0050  0.0015  -0.0065          0.0015 1000
Q             0.1061   0.1204  0.0905   0.0299          0.0762 1000
B             0.1947   0.1915  0.1977  -0.0062          0.2009 1000
g             0.8378   0.9011  0.7963   0.1048          0.7331 1000
gp            0.1673   0.1757  0.1596   0.0161          0.1511 1000

Let's concentrate on the $D_{xy}$ which is at the top. The first column denotes the original index, which was $0.3857$. The column called optimism denotes the amount of estimated overestimation by the model. The column index.corrected is the original estimate minus the optimism. In this case, the bias-corrected $D_{xy}$ is a bit smaller than the original. The bias-corrected c-index (AUC) is $c=\frac{1+ D_{xy}}{2}=0.6749$.

We can also calculate a calibration curve using resampling:

#-----------------------------------------------------------------------------
# Calibration curve using bootstrap
#-----------------------------------------------------------------------------

my.calib <- calibrate(mylogit, method="boot", B=1000)

par(bg="white", las=1)
plot(my.calib, las=1)

n=400   Mean absolute error=0.016   Mean squared error=0.00034
0.9 Quantile of absolute error=0.025

LogReg Calibration

The plot provides some evidence that our models is overfitting: the model underestimates low probabilities and overestimates high probabilities. There is also a systematic overestimation around $0.3$.

Predictive model building is a big topic and I suggest reading Frank Harrell's course notes.

$\endgroup$
5
  • 5
    $\begingroup$ Thanks so much for the wonderful examples and how-to for the rms package @COOLSerdash and Nick. $\endgroup$ – Frank Harrell Jul 19 '13 at 11:18
  • $\begingroup$ Thank you very much for your answer, it also helped me a lot ! I have a quick question, what do you mean by systematic over estimation around 0.3 ? And can you please have a look at my graph ? s23.postimg.org/9cucdg6e3/calibration_curve.png How should I interpret this output $\endgroup$ – CanCeylan Jul 23 '13 at 12:59
  • 2
    $\begingroup$ @CanCeylan In my graph, there is a small curve in the line around 0.3 (the line goes below the ideal line). Regarding your graph: Please add your graph to your original question so that others can see it too. It seems as your model is underestimating the probabilities systematically over the whole space. $\endgroup$ – COOLSerdash Jul 23 '13 at 13:06
  • $\begingroup$ @COOLSerdash thanks. I've also updated my question. But from there, how should I decide whether to cancel the model or use it ? Because c value is not bad(0.7755), Mean absolute error(0.002) and Mean squared error(5e-05) are so small, but the calibration graph looks like bad. \ $\endgroup$ – CanCeylan Jul 23 '13 at 14:15
  • $\begingroup$ @CanCeylan Thanks for updating your question. Hm, I'm unsure. The calibration plot looks not too good. I would wait for more opinions (from Frank Harrell or others). I don't feel confident enough to give a definitive recommendation, sorry. $\endgroup$ – COOLSerdash Jul 23 '13 at 14:20
6
$\begingroup$

A note on interpretation of coefficients: do recall they depend on how the predictors are written as numbers. So for continuous variables they depend on the units in which they are measured; for categorical predictors, the coding scheme. Don't be tempted to think that, say, A9 is 'unimportant' just because its coefficient of 0.003453 is small—A9 might range over several orders of magnitude in some population of interest while the other predictors vary only slightly, or it may be easy to set to very high or low values while the others are hard to change much.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.