# Inference about the outcomes of two logistic regressions

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?

-
 Are the models nested? Meaning can one model be considered a subset of the other model, or are they entirely different? Similarly, are they being run on the same data? – EpiGrad Aug 24 '11 at 23:09 maybe @Frank Harrel can answer my question. very nice package btw – Dbr Aug 24 '11 at 23:10 @EpiGrad the models are different, but run on the same data – Dbr Aug 24 '11 at 23:10

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.

-
 You're right that AICs are not formal tests but in this case, since the dimension of the parameter space is the same, comparing the AICs is actually equivalent to comparing the converged log-likelihoods, which is a defensible measure of fit (particularly when one model's converged log-likelihood is 70 points higher from a moderate sample size of around 300). – Macro Aug 25 '11 at 0:15 @EpiGrad that's exactly what I would like to see, whether the model fits poorly or not. How can I calculate the log likelihood? – Dbr Aug 25 '11 at 8:12 @ Daniel. Bluntly, in R, I have no idea. I'm a SAS man myself, outside of the occasional visualization problem. But nothing you've talked about (including R squared) is a singular solution to "How well does the model fit". Even the log likelihood is mainly useful to compare two models. Seeing how well a model "fits" really does involve multiple steps, including looking at your residuals for patterns of poor fits, particularly poor observations, etc. – EpiGrad Aug 25 '11 at 19:18 @EpiGrad - couldn't you just take the AIC for the saturated model for an "absolute" comparison. This basically means doing the deviance residual test I describe. For the AIC difference between the current model with, say $k$, degrees of freedom and a model with one parameter for each of, say $N$, observations we have that $AIC_{sat}-AIC_{model}=2df-deviance$. This means that AIC will favour the current model over the saturated model if $deviance<2df$ – probabilityislogic Mar 23 '12 at 11:08

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

-

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

-
 This is plain wrong about the deviance residuals not incorporating the model fit, as the deviance is given by $2\sum_i y_i\log\left(\frac{y_i}{n_i\hat{p}_i}\right) + (n_i-y_i)\log\left(\frac{n_i-y_i}{n_i-n_i\hat{p}_i}\right)$. Now if all the $n_i=1$, then the deviance reduces to $-2\sum_{i:y_i=0}\log(\hat{p}_i)-2\sum_{i:y_i=1}\log(1-\hat{p}_i)$. Basically, you measure how close to zero the probabilities for observed $0$s are, and how close to one the probabilities for observed $1$s are. How is this not a proper measure of fit? – probabilityislogic Nov 24 '11 at 9:40 Sorry, that should be the other way round, first sum should be over $y_i=1$ and second should be over $y_i=0$. So you measure how close to zero the probabilities for observed $1$s (which makes more sense, as these are residuals). – probabilityislogic Nov 24 '11 at 9:54