Inference about the outcomes of two logistic regressions

Question

I have run two separate logistic regressions and would like to assess which model fits the data better. Each model has 1 predictor.

Here's the output for both:

Model 1:

Deviance Residuals: 
Min      1Q  Median      3Q     Max  
-1.280  -1.046  -1.046   1.078   1.315  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)  
(Intercept)  -0.3185     0.1643  -1.938   0.0526 .
L1            0.5564     0.2317   2.402   0.0163 *

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 421.32  on 303  degrees of freedom
Residual deviance: 415.49  on 302  degrees of freedom
AIC: 419.49

Model 2:

glm(formula = Response ~ L2, family = binomial, data = log)

Deviance Residuals: 
Min       1Q   Median       3Q      Max  
-1.6651  -0.7235  -0.7235   0.7585   1.7138  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.2068     0.1927  -6.264 3.75e-10 ***
 L2            2.3054     0.2687   8.580  < 2e-16 ***

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 421.32  on 303  degrees of freedom
Residual deviance: 334.99  on 302  degrees of freedom
AIC: 338.99

Number of Fisher Scoring iterations: 4

Is there something similar to an R square for logistic regression that I could use to assess how well it fits, and also that I can compare to the other model?

Similarly if I run the lrm command from the rms package, which discrimination index is best to look at and how can I compare the two models?

Answer 1

As your models are run on the same data, but aren't nested, you're probably best off using the AIC (Akaike information criterion), which is a measure of model fit - arguably a better one than R square, which has all manner of problems.

For AIC, a lower number is better, but there's no formal tests for it - you'll still have to make a subjective evaluation of what you think "better" means if you're choosing between two models, but the AIC values should at least provide some guidance.

You should also look at the residuals of your model fits and the like, because one of the things AIC will not tell you is if all your models fit poorly. It can just tell you if one fits relatively better than another.

Answer 2

Comparing models

In terms of comparing the models, because they have the same number of parameters, we can simply compare the residual deviance. The models are equally "complex", so we need not penalise for complexity we comparing them to each other. Now, because the residual deviance is twice the log-likelihood ratio between the fitted model and the saturated model, the decrease in deviance is twice the log-likelihood ratio between the two models. So we have a drop in deviance of $415.49-334.99=80.5$. To convert to an odds ratio we exponentiate half this number, to give $\exp(80.5/2)=3.02\times 10^{17}$. This means that the data are $3.02\times 10^{17}$ times more likely to have occured under model 2 compared to model 1. So I would have no hesitation is declaring model 2 the better model, compared to model 1. But there may be a still better model, "model 3", which you haven't thought of yet: that is where goodness of fit tests come in handy.

Assess goodness of fit

For assessing the goodness of fit for the models, one way is to use the fact that, for the exponential family (of which the binomial is a member), we have approximately:

$$\frac{D_{res}}{\phi} \sim\chi^2_{df}$$

Where $D_{res}$ is the residual deviance, output for your model, and $df$ is the degrees of freedom, and $\chi^2_{df}$ is a chi-square distribution with $df$ degrees of freedom. Now for a binomial glm, $\phi=1$ (which is the comment "dispersion parameter for binomial family taken to be $1$"). This means that we would expect to see a residual deviance about equal to the degrees of freedom. For your example, this means $E[D_{res}]\approx 302$. We can use the standard deviation of the chi-square distributions, equal to $\sqrt{Var[D_{res}]}\approx\sqrt{2 df}$, and create an upper bound on an acceptable residual deviance as $df+3\sqrt{2 df}$ ("3 sigma" rule). For your example, this gives $302+3\sqrt{604}=351$. This shows that your first model exhibits much greater dispersion than what you would expect from a binomial model, as the residual deviance is $415>>351$. If we use the chi-square distribution to generate a p-value, we get $Pr(D_{res}>415)\approx 0.00002$ (the approx. is due to using the chi-square). Your second model slightly over-dispersed, but not overly inconsistent with a binomial distribution, because $302<D_{res}<351$. The corresponding p-value is $Pr(D_{res}>335)\approx 0.09$. However, given that $D_{res}>df$, you may find a better one variable model by transforming the covariate $L2$. The maximum amount of support theoretically possible for this "model 3" over model 2 is equal to $\exp(335/2)=5.55\times 10^{72}$, although this is very hard to achieve in practice (physical constraints aside). You are more likely to be in the range of $\exp([335-302]/2)=1.47\times 10^7$ - i.e. the best one-variable alternative is likely to give residual deviance equal to the degrees of freedom.

You can reasonably trust both deviance tests because your deviance residuals are all between $\pm 2$, so this indicates there are no observations which have been poorly fitted (outliers). Also, looking at your beta estimates suggests that the response is reasonably common in your data set, which also helps make the chi-square approximation error small (ala normal approximation to the binomial).

Answer 3

The second model appears to fit the data better based on the considerably smaller value of the AIC (and since there is the same number of parameters - this implies a considerably larger value of the converged log-likelihood of the model). Here I am implicitly using the log-likelihood as a measure of fit although there are alternatives (e.g. the Hosmer-Lemeshow index of goodness of fit for binary data - http://en.wikipedia.org/wiki/Hosmer%E2%80%93Lemeshow_test)

Note: You can see the residual deviance is much smaller for the second model but one should be cautioned that, while it is tempting, this should not be interpreted as an index of goodness of fit for binary data, because the residual deviance reduces to

$$ \sum_{i} \hat{p}_{i} \log \Big( \hat{p}_{i}/(1-\hat{p}_{i}) \Big) + \log(1-\hat{p}_{i}) $$

which only involves the model implies probabilities $\hat{p}_{i}$, and not the observed outcomes, so it can't be a proper measure of fit (although it is for binomial data in general other than binary).