How to calculate pseudo-$R^2$ from R's logistic regression?

Question

Christopher Manning's writeup on logistic regression in R shows a logistic regression in R as follows:

ced.logr <- glm(ced.del ~ cat + follows + factor(class), 
  family=binomial)

Some output:

> summary(ced.logr)
Call:
glm(formula = ced.del ~ cat + follows + factor(class),
    family = binomial("logit"))
Deviance Residuals:
Min            1Q    Median       3Q      Max
-3.24384 -1.34325   0.04954  1.01488  6.40094

Coefficients:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)   -1.31827    0.12221 -10.787 < 2e-16
catd          -0.16931    0.10032  -1.688 0.091459
catm           0.17858    0.08952   1.995 0.046053
catn           0.66672    0.09651   6.908 4.91e-12
catv          -0.76754    0.21844  -3.514 0.000442
followsP       0.95255    0.07400  12.872 < 2e-16
followsV       0.53408    0.05660   9.436 < 2e-16
factor(class)2 1.27045    0.10320  12.310 < 2e-16
factor(class)3 1.04805    0.10355  10.122 < 2e-16
factor(class)4 1.37425    0.10155  13.532 < 2e-16
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 958.66 on 51 degrees of freedom
Residual deviance: 198.63 on 42 degrees of freedom
AIC: 446.10
Number of Fisher Scoring iterations: 4

He then goes into some detail about how to interpret coefficients, compare different models, and so on. Quite useful.

However, how much variance does the model account for? A Stata page on logistic regression says:

Technically, $R^2$ cannot be computed the same way in logistic regression as it is in OLS regression. The pseudo-$R^2$, in logistic regression, is defined as $1 - \frac{L1}{L0}$, where $L0$ represents the log likelihood for the "constant-only" model and $L1$ is the log likelihood for the full model with constant and predictors.

I understand this at the high level. The constant-only model would be without any of the parameters (only the intercept term). Log likelihood is a measure of how closely the parameters fit the data. In fact, Manning sort of hints that the deviance might be $-2 \log L$. Perhaps null deviance is constant-only and residual deviance is $-2 \log L$ of the model? However, I'm not crystal clear on it.

Can someone verify how one actually computes the pseudo-$R^2$ in R using this example?

Answer 1

Don't forget the rms package, by Frank Harrell. You'll find everything you need for fitting and validating GLMs.

Here is a toy example (with only one predictor):

set.seed(101)
n <- 200
x <- rnorm(n)
a <- 1
b <- -2
p <- exp(a+b*x)/(1+exp(a+b*x))
y <- factor(ifelse(runif(n)<p, 1, 0), levels=0:1)
mod1 <- glm(y ~ x, family=binomial)
summary(mod1)

This yields:

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.8959     0.1969    4.55 5.36e-06 ***
x            -1.8720     0.2807   -6.67 2.56e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 258.98  on 199  degrees of freedom
Residual deviance: 181.02  on 198  degrees of freedom
AIC: 185.02

Now, using the lrm function,

require(rms)
mod1b <- lrm(y ~ x)

You soon get a lot of model fit indices, including Nagelkerke $R^2$, with print(mod1b):

Logistic Regression Model

lrm(formula = y ~ x)

                      Model Likelihood     Discrimination    Rank Discrim.    
                         Ratio Test            Indexes          Indexes       

Obs           200    LR chi2      77.96    R2       0.445    C       0.852    
 0             70    d.f.             1    g        2.054    Dxy     0.705    
 1            130    Pr(> chi2) <0.0001    gr       7.801    gamma   0.705    
max |deriv| 2e-08                          gp       0.319    tau-a   0.322    
                                           Brier    0.150                     


          Coef    S.E.   Wald Z Pr(>|Z|)
Intercept  0.8959 0.1969  4.55  <0.0001 
x         -1.8720 0.2807 -6.67  <0.0001 

Here, $R^2=0.445$ and it is computed as $\left(1-\exp(-\text{LR}/n)\right)/\left(1-\exp(-(-2L_0)/n)\right)$, where LR is the $\chi^2$ stat (comparing the two nested models you described), whereas the denominator is just the max value for $R^2$. For a perfect model, we would expect $\text{LR}=2L_0$, that is $R^2=1$.

By hand,

> mod0 <- update(mod1, .~.-x)
> lr.stat <- lrtest(mod0, mod1)
> (1-exp(-as.numeric(lr.stat$stats[1])/n))/(1-exp(2*as.numeric(logLik(mod0)/n)))
[1] 0.4445742
> mod1b$stats["R2"]
       R2 
0.4445742 

Ewout W. Steyerberg discussed the use of $R^2$ with GLM, in his book Clinical Prediction Models (Springer, 2009, § 4.2.2 pp. 58-60). Basically, the relationship between the LR statistic and Nagelkerke's $R^2$ is approximately linear (it will be more linear with low incidence). Now, as discussed on the earlier thread I linked to in my comment, you can use other measures like the $c$ statistic which is equivalent to the AUC statistic (there's also a nice illustration in the above reference, see Figure 4.6).

Answer 2

To easily get a McFadden's pseudo $R^2$ for a fitted model in R, use the "pscl" package by Simon Jackman and use the pR2 command. http://cran.r-project.org/web/packages/pscl/index.html

Answer 3

Before calculate the pseudo-$R^2$ for your logistic regression, I want to ask you, do you think it is appropriate to use pseudo-$R^2$ measures (McFadden’s, McKelvey-Zavoina’s)?

Firstly, the McFadden’s Pseudo-$R^2$ of logistic model do not imply the proportion of the variance of the response explained by explanatory variables at all. But one of the purposes of developing a goodness-of-fit is this property. So I want to go deeper and broader for you to think what could be a good goodness-of-fit measure for logit/probit model.

Before answering this question, I want to emphasis the “comparability” first. This is the reason why the OLS $R^2$ is such an extensively used measure of goodness of fit. It is useful because we can compare the $R^2$ of different models to get an idea about how each model performs to fit the data and what are the adequacy of them.

Another aspect of the “comparability” is when we compare models across different samples! consider the study of job satisfaction where data may be collected in three different fashions using (a) a quantitative score on a scale of 0-100; (b) a dichotomous indicator Yes/No; or (c) a five-category rating ranging from extremely unsatisfied to extremely satisfied. Although neither the samples nor the empirical models used to draw inferences are the same, “most empirical researchers are explicitly or implicitly making rough comparisons of ‘goodness of fit’ across” these models and samples, because they address similar domain questions (Veall and Zimmermann, 1996). In social studies, “the research experience in the area is far more important than any statistical criteria” such as what specific method is used to select variables (Veall and Zimmermann, 1996).

We believe it is vital to have a goodness-of-fit measure for logit/probit models that is analogous to the OLS R2, as it will ensure comparability of different models across different samples for similar research questions. There is no existing pseudo-$R^2$ known so far that can meet all the needs.

To solve this issue, we have developed a new goodness-of-fit measure to resemble the OLS R-squared that can also implies the proportion of the variance of the surrogate response explained by explanatory variables. We will have a new paper "A new goodness-of-fit measure for probit models: surrogate $R^2$" coming soon with all the details. Please check on my webpage: https://homepages.uc.edu/~zhuxr/. Then package will be also available soon.

---

Old answer:

Be careful with the calculation of Pseudo-$R^2$:

McFadden’s Pseudo-$R^2$ is calculated as $R^2_M=1- \frac{ln\hat{L}_{full}}{ln\hat{L}_{null}}$, where $ln\hat{L}_{full}$ is the log-likelihood of full model, and $ln\hat{L}_{full}$ is log-likelihood of model with only intercept.

Two approach to calculate Pseudo-$R^2$:

1. Use deviance: since $deviance = -2*ln(L_{full})$, $null.deviance = -2*ln(L_{null})$

   pR2 = 1 - mod$deviance / mod$null.deviance # works for glm

But the above approach doesn't work for out-of-sample Pseudo $R^2$

2. Use "logLik" function in R and definition(also works for in-sample)

   mod_null <- glm(y~1, family = binomial, data = insample) 1- logLik(mod)/logLik(mod_null)

This can be slightly modified to compute out-of-sample Pseudo $R^2$

Example:

out-of-sample pseudo-R

Usually, the out-of-sample pseudo-$R^2$ is calculated as $$R_p^2=1−\frac{L_{est.out}}{L_{null.out}},$$ where $L_{est.out}$ is the log likelihood for the out-of-sample period based on the estimated coefficients of in-sample period, while and $L_{null.out}$ is the log likelihood for intercept-only model for the out-of-sample period.

Codes:

pred.out.link <- predict(mod, outSample, type = "link")
mod.out.null <- glm(Default~1, family = binomial, data = outSample)
pR2.out <- 1 - sum(outSample$y * pred.out.link - log(1 + exp(pred.out.link))) / logLik(mod.out.null)

Answer 4

if deviance were proportional to log likelihood, and one uses the definition (see for example McFadden's here)

pseudo R^2 = 1 - L(model) / L(intercept)

then the pseudo-$R^2$ above would be $1 - \frac{198.63}{958.66}$ = 0.7928

The question is: is reported deviance proportional to log likelihood?

Answer 5

If its out of sample, then I believe the $R^2$ must be computed with the according log-likelihoods as $R^2=1-\frac{ll_{full}}{ll_{constant}}$, where $ll_{full}$ is the log-likelihood of the test data with the predictive model calibrated on the training set, and $ll_{constant}$ is the log-likelihood of the test data with a model with just a constant fitted on the training set, and then use the fitted constant to predict on the testing set computing the probabilities and therefore get the log-likelihood.

Note that in a linear regression, is analogous, the out of sample $R^2$ is computed as $R^2=1-\frac{\sum_{i}(y_{i}-\hat{y}_i)^2}{\sum_{i}(y_{i}-\overline{y}_{train})^2}$, where in particular if we look at the denominator term $\sum_{i}(y_{i}-\overline{y}_{train})^2$, the prediction uses the average over the training set, $\overline{y}_{train}$. This is like if we fit a model in the training data with just a constant, so we have to minimize $\sum_{i}(y_i-\beta_0)^2$, which results in $\hat{\beta}_0=\overline{y}_{train}$, then, this plain constant predictive model is the one used as benchamrk (i.e. in the denominator of the oos $R^2$ term) for the computation of the out of sample $R^2$.