Which pseudo-$R^2$ measure is the one to report for logistic regression (Cox & Snell or Nagelkerke)?

Question

I have SPSS output for a logistic regression model. The output reports two measures for the model fit, Cox & Snell and Nagelkerke.

So as a rule of thumb, which of these $R^²$ measures would you report as the model fit?

Or, which of these fit indices is the one that is usually reported in journals?

---

Some Background: The regression tries to predict the presence or absence of a bird (capercaillie) from some environmental variables (e.g., steepness, vegetation cover, ...). Unfortunately, the bird did not appear very often (35 hits to 468 misses) so the regression performs rather poorly. Cox & Snell is .09, Nagelkerke, .23.

The subject is environmental sciences or ecology.

Answer 1

Normally I wouldn't report $R^2$ at all. Hosmer and Lemeshow, in their textbook Applied Logistic Regression (2nd Ed.), explain why:

In general, [$R^2$ measures] are based on various comparisons of the predicted values from the fitted model to those from [the base model], the no data or intercept only model and, as a result, do not assess goodness-of-fit. We think that a true measure of fit is one based strictly on a comparison of observed to predicted values from the fitted model.

[At p. 164.]

Concerning various ML versions of $R^2$, the "pseudo $R^2$" stat, they mention that it is not "recommended for routine use, as it is not as intuitively easy to explain," but they feel obliged to describe it because various software packages report it.

They conclude this discussion by writing,

...low $R^2$ values in logistic regression are the norm and this presents a problem when reporting their values to an audience accustomed to seeing linear regression values. ... Thus [arguing by reference to running examples in the text] we do not recommend routine publishing of $R^2$ values with results from fitted logistic models. However, they may be helpful in the model building state as a statistic to evaluate competing models.

[At p. 167.]

My experience with some large logistic models (100k to 300k records, 100 - 300 explanatory variables) has been exactly as H & L describe. I could achieve relatively high $R^2$ with my data, up to about 0.40. These corresponded to classification error rates between 3% and 15% (false negatives and false positives, balanced, as confirmed using 50% hold-out datasets). As H & L hinted, I had to spend a lot of time disabusing the client (a sophisticated consultant himself, who was familiar with $R^2$) concerning $R^2$ and getting him to focus on what mattered in the analysis (the classification error rates). I can warmly recommend describing the results of your analysis without reference to $R^2$, which is more likely to mislead than not.

Answer 2

Both indices are measures of strength of association (i.e. whether any predictor is associated with the outcome, as for an LR test), and can be used to quantify predictive ability or model performance. A single predictor may have a significant effect on the outcome but it might not necessarily be so useful for predicting individual response, hence the need to assess model performance as a whole (wrt. the null model). The Nagelkerke $R^2$ is useful because it has a maximum value of 1.0, as Srikant said. This is just a normalized version of the $R^2$ computed from the likelihood ratio, $R^2_{\text{LR}}=1-\exp(-\text{LR}/n)$, which has connection with the Wald statistic for overall association, as originally proposed by Cox and Snell. Other indices of predictive ability are Brier score, the C index (concordance probability or ROC area), or Somers' D, the latter two providing a better measure of predictive discrimination.

The only assumptions made in logistic regression are that of linearity and additivity (+ independence). Although many global goodness-of-fit tests (like the Hosmer & Lemeshow $\chi^2$ test, but see my comment to @onestop) have been proposed, they generally lack power. For assessing model fit, it is better to rely on visual criteria (stratified estimates, nonparametric smoothing) that help to spot local or global departure between predicted and observed outcomes (e.g. non-linearity or interaction), and this is largely detailed in Harrell's RMS handout. On a related subject (calibration tests), Steyerberg (Clinical Prediction Models, 2009) points to the same approach for assessing the agreement between observed outcomes and predicted probabilities:

Calibration is related to goodness-of-fit, which relates to the ability of a model to fit a given set of data. Typically, there is no single goodness-of-fit test that has good power against all kinds of lack of fit of a prediction model. Examples of lack of fit are missed non-linearities, interactions, or an inappropriate link function between the linear predictor and the outcome. Goodness-of-fit can be tested with a $\chi^2$ statistic. (p. 274)

He also suggests to rely on the absolute difference between smoothed observed outcomes and predicted probabilities either visually, or with the so-called Harrell's E statistic.

More details can be found in Harrell's book, Regression Modeling Strategies (pp. 203-205, 230-244, 247-249). For a more recent discussion, see also

Steyerberg, EW, Vickers, AJ, Cook, NR, Gerds, T, Gonen, M, Obuchowski, N, Pencina, MJ, and Kattan, MW (2010). Assessing the Performance of Prediction Models, A Framework for Traditional and Novel Measures. Epidemiology, 21(1), 128-138.

Answer 3

I would have thought the main problem with any kind of $R^2$ measure for logistic regression is that you are dealing with a model which has a known noise value. This is unlike standard linear regression, where the noise level is usually treated as unknown. For we can write a glm probability density function as:

$$f(y_i|\mu_i,\phi)=\exp\left(\frac{y_ib(\mu_i)-c(\mu_i)}{\phi}+d(y_i,\phi)\right)$$

Where $b(.),\ c(.),\ d(.;.)$ are known functions, and $\mu_i=g^{-1}(x_i^T\beta)$ for inverse link function $g^{-1}(.)$. If we define the usual GLM deviance residuals as

\begin{align} d_i^2 &= 2\phi\left(\log[f(y_i|\mu_i=y_i,\phi)]-\log[f(y_i|\mu_i=\hat{\mu}_i,\phi)]\right) \\ &= 2\phi \left[y_ib(y_i)-y_ib(\hat{\mu}_i)-c(y_i)+c(\hat{\mu}_i)\right] \end{align} The we have (via likelihood ratio chi-square, $\chi^2=\frac{1}{\phi}\sum_{i=1}^{N}d_i^2$)

$$E\left(\sum_{i=1}^{N}d_i^2\right)=E(\phi\chi^2)\approx (N-p)\phi$$

Where $p$ is the dimension of $\beta$. For logistic regression we have $\phi=1$, which is known. So we can use this to decide on a definite level of residual that is "acceptable" or "reasonable". This usually cannot be done for OLS regression (unless you have prior information about the noise). Namely, we expect each deviance residual to be about $1$. Too many $d_i^2\gg1$ and it is likely that an important effects are missing from the model (under-fitting); too many $d_i^2\ll1$ and it is likely that there are redundant or spurious effects in the model (over-fitting). (these could also mean model mispecification).

Now this means that the problem for the pseudo-$R^2$ is that it fails to take into account that the level of binomial variation is predictable (provided the binomial error structure isn't being questioned). Thus even though Nagelkerke ranges from $0$ to $1$, it is still not scaled properly. Additionally, I can't see why these are called pseudo $R^2$ if they aren't equal to the usual $R^2$ when you fit a "GLM" with an identity link and normal error. For example, the equivalent cox-snell R-squared for normal error (using REML estimate of variance) is given by:

$$R^2_{CS}=1-\exp\left(-\frac{N-p}{N}\cdot \frac{R^2_{OLS}}{1-R^2_{OLS}}\right)$$

Which certainly looks strange.

I think the better "Goodness of Fit" measure is the sum of the deviance residuals, $\chi^2$. This is mainly because we have a target to aim for.

Answer 4

I found Tue Tjur's short paper "Coefficients of Determination in Logistic Regression Models - A New Proposal: The Coefficient of Discrimination" (2009, The American Statistician) on various proposals for a coefficient of determination in logistic models quite enlightening. He does a good job highlighting pros and cons - and of course offers a new definition. Very much recommended (though I have no favorite myself).

Answer 5

I was also going to say 'neither of them', so i've upvoted whuber's answer.

As well as criticising R^2, Hosmer & Lemeshow did propose an alternative measure of goodness-of-fit for logistic regression that is sometimes useful. This is based on dividing the data into (say) 10 groups of equal size (or as near as possible) by ordering on the predicted probability (or equivalently, the linear predictor) then comparing the observed to expected number of positive responses in each group and performing a chi-squared test. This 'Hosmer-Lemeshow goodness-of-fit test' is implemented in most statistical software packages.

Answer 6

Despite the arguments against using pseudo-r-squareds, some people will for various reasons want to continue using them at least at certain times. What I have internalized from my readings (and I'm sorry I cannot provide citations at the moment) is that

• if both C&S and Nag. are below .5, C&S will be a better gauge;
  if they're both above .5, Nag. will; and
  if they straddle .5, punt.

Also, a formula whose results often fall between these two, mentioned by Scott Menard in Applied Logistic Regression Analysis (Sage), is

[-2LL0 - (-2LL1)]/-2LL0.

This is denoted as "L" in the chart below.

Answer 7

I would agree in general that just using R2 is not good. But also see that point of @rolando2 comments that focusing on classification metrics could be not enough while comparing the models.

I guess my contribution to the discussion is that I think that several measures are to be reported to assess different model qualities.

For example, one may not only want to know what happens at the optimum threshold which separates cases from non-cases and respective error rates (false positives etc and integrated measure of this which is c-statistic), but also how good model output in terms of probabilities of being a case or non-case is close to reality (or corresponds to the actual rate of cases in the pool of observations with given parameters) across ALL ranges of output. I.e. that if it says an observation is 20% likely to have output of 1, that around 20% of observations with similar risk factors are cases. This is what referred to as calibration quality as opposed to discrimination (e.g. in Steyerberg et al Assessing the performance of prediction models: a framework for traditional and novel measures.) And this is exactly what calibration plot visually assesses, while Hosmer & Lemeshow quantified it in their goodness-of-fit statistics. The critics of H-L is that it depends on how you group observations, and playing with it I see that the value can change quite a bit on different data. Calibration slopes and intercept may be good alternative if one also does cross-validation or evaluates the model on a hold-out data as in predictive modelling. Finally, some kind of R2, preferably Brier score or scaled Brier score, can be used to assess the overall fit. This measure is the same as R2 for linear regression, where error is defined as diff in probability and binary output, but also takes into account known variance of the binary output and normalises by q*(1-q). Also, for a narrower discussion on checking of adding a new predictor makes a better model, IDI - integrated discrimination improvement, which is similar to the difference in scaled Brier scores, could be a very good addition to, say, the change in c-statistics or how good reclassification was - as it checks how re-classification got better across all thresholds. (Pencina MJ, D’Agostino RB, D’Agostino RB, Vasan RS. Evaluating the added predictive ability of a new marker: from area under the ROC curve to reclassification and beyond)

Answer 8

I would prefer the Nagelkerke as this model fit attains 1 when the model fits perfectly giving the reader a sense of how far your model is from perfect fit. The Cox & Shell does not attain 1 for perfect model fit and hence interpreting a value of 0.09 is a bit harder. See this url for further info on Pseudo RSquared for an explanation of various types of fits.

Answer 9

Instead of the Nagelkerke way of scaling $R^2$ to allow a 1.0 to be attained, I prefer to substitute the effective sample size for $N$ in the $R^2$ formula. This will not reach 1.0 for perfect binary predictions but this approach translates to other settings such as survival analysis where often the effective $N$ is the number of events, and to ordinal regression. See https://hbiostat.org/bib/r2.html.

Of those my favorite is the modified Maddala-Cox-Snell $R^{2}_{m,p}$ which uses effective sample size $m$ and penalizes for $p$ covariates. In the normal linear model this is almost exactly the traditional $R^{2}_{\mathrm{adj}}$.