94
$\begingroup$

In answering this question John Christie suggested that the fit of logistic regression models should be assessed by evaluating the residuals. I'm familiar with how to interpret residuals in OLS, they are in the same scale as the DV and very clearly the difference between y and the y predicted by the model. However for logistic regression, in the past I've typically just examined estimates of model fit, e.g. AIC, because I wasn't sure what a residual would mean for a logistic regression. After looking into R's help files a little bit I see that in R there are five types of glm residuals available, c("deviance", "pearson", "working","response", "partial"). The help file refers to:

I do not have a copy of that. Is there a short way to describe how to interpret each of these types? In a logistic context will sum of squared residuals provide a meaningful measure of model fit or is one better off with an Information Criterion?

$\endgroup$
10
  • 2
    $\begingroup$ There are elements to this question that remain unanswered, e.g. the nature of the "pearson", "working","response", and "partial" residuals, but for now I will accept Thylacoleo's answer. $\endgroup$ Commented Aug 12, 2010 at 7:12
  • 2
    $\begingroup$ I find the binnedplot function in the R package arm gives a very helpful plot of residuals. It's described nicely on p.97-101 of Gelman and Hill 2007. $\endgroup$ Commented Nov 1, 2010 at 14:03
  • 2
    $\begingroup$ One really easy way to check model fit is a plot of the observed vs the predicted proportions. But this won't work if you have bernoulli regression (i.e. all of your observations have unique combinations of the independent variables, so that $n_i=1$), because you will just see a line of zeros and ones. $\endgroup$ Commented Feb 17, 2011 at 13:17
  • $\begingroup$ Yeah - sadly I usually am using a Bernoulli DV. $\endgroup$ Commented Feb 17, 2011 at 20:05
  • 2
    $\begingroup$ See also Understanding glm$residuals and resid(glm) on Stack Overflow. $\endgroup$ Commented Mar 17, 2018 at 19:18

5 Answers 5

49
$\begingroup$

The easiest residuals to understand are the deviance residuals as when squared these sum to -2 times the log-likelihood. In its simplest terms logistic regression can be understood in terms of fitting the function $p = \text{logit}^{-1}(X\beta)$ for known $X$ in such a way as to minimise the total deviance, which is the sum of squared deviance residuals of all the data points.

The (squared) deviance of each data point is equal to (-2 times) the logarithm of the difference between its predicted probability $\text{logit}^{-1}(X\beta)$ and the complement of its actual value (1 for a control; a 0 for a case) in absolute terms. A perfect fit of a point (which never occurs) gives a deviance of zero as log(1) is zero. A poorly fitting point has a large residual deviance as -2 times the log of a very small value is a large number.

Doing logistic regression is akin to finding a beta value such that the sum of squared deviance residuals is minimised.

This can be illustrated with a plot, but I don't know how to upload one.

$\endgroup$
5
  • $\begingroup$ I've corrected an error in my original answer. I first wrote p=logit(Xbeta). In fact the predicted probability is the inverse logit of the linear combination, p=inv-logit(Xbeta). In R this is calculated as p<-plogit(Xbeta), which is p=exp(Xbeta)/(1+exp(X*beta)). $\endgroup$
    – Thylacoleo
    Commented Aug 10, 2010 at 1:07
  • 1
    $\begingroup$ Which R package is plogit from? It wasn't clear if you were defining it here or getting it from somewhere else. $\endgroup$
    – Amyunimus
    Commented Jan 31, 2013 at 18:54
  • 1
    $\begingroup$ @Amyunimus plogit is in R (stats), no package required (at least not anymore) $\endgroup$ Commented Jul 7, 2014 at 7:26
  • 1
    $\begingroup$ Maximum likelihood estimation with binary Y does not minimize any sort of sum of squares. $\endgroup$ Commented Sep 5, 2020 at 12:06
  • $\begingroup$ @FrankHarrell +1 to your comment, not sure how this answer got accepted. Maximum likelihood minimizes cross entropy for binary response. $\endgroup$
    – iggy
    Commented Nov 22, 2022 at 20:06
26
$\begingroup$

Response: $$y_i - \hat\mu_i$$ response residuals are inadequate for assessing a fitted glm, because GLMs are based on distributions where (in general) the variance depends on the mean.

Pearson:

The most direct way to handle the non-constant variance is to divide it out: $$ \frac{y_i - \hat\mu_i}{\sqrt{V(\mu_i)|_{\hat\mu_i}}}$$ where $V()$ is the (GLM) variance function ($Var(y_i) = a(\phi)*V(\mu_i)$)

Under "Small dispersion asymptotics" conditions, the Pearson residuals have an approximate normal distribution.

Deviance: $$sign(y_i-\hat\mu_i)*\sqrt{d_i}$$ where $d_i$ is the unit deviance, i.e. $d_i = 2(t(y_i,y_i)-t(y_i,\hat\mu_i))$

The deviance statistic (sum of squared unit-deviances) has an approximate chi-square distribution (when the saddlepoint approximation applies and under "Small dispersion asymptotics" conditions). Under these same conditions, the deviance residuals have an approximate normal distribution.

Working: $$z_i - \eta_i $$ where $z_i$ are the working responses $\eta_i + \frac{d\eta_i}{d\mu_i}(y_i-\hat\mu_i)$ and $\eta_i$ is the linear predictor. Meaning you get that the residual is $\frac{d\eta_i}{d\mu_i}(y_i-\hat\mu_i)$.

The model coefficients are fitted using Fisher scoring algorithm / Iterative Reweighted Least Square (IRLS). And it can be shown that each iteration of this algorithm is equivalent to doing ordinary least-squares on the working responses as defined here.

To test the link function - plotting the linear predictor against the working responses should come out linear if the right link function was used.

Partial: $$z_i - \eta_i + X^*\beta$$

where $X^*$ is the centered $X$. Partial residuals can be used to determine if a covariate/predictor is on an inappropriate scale.

Quantile: $$\Phi^{-1}(F(y_i))$$

Where $F(y_i)$ is the CDF of $y_i$, and $\Phi^{-1}$ is the quantile function of standard normal (inverse CDF). For discrete $y_i$'s you take $u \sim Unif(F(y_i-1), F(y_i))$ and $\Phi^{-1}(u)$.

Here is an example code to calculate these residuals:

Y = c(0,0,0,0,1,1,1,1,1)
x1 = c(1,2,3,1,2,2,3,3,3)
x2 = c(1,0,0,1,0,0,0,0,0)

fit = glm(y ~ x1 + x2, family = 'binomial')

lp = predict(fit)
mu = exp(lp)/(1+exp(lp))

# manually calculating the 1st response residual
resid(fit, type="response")[1]
Y[1] - mu[1]

# manually calculating the 1st pearson residual
resid(fit, type="pearson")[1]
(Y[1]-mu[1]) / sqrt(mu[1]*(1-mu[1]))

# manually calculating the 1st deviance residual
resid(fit, type="deviance")[1]
sqrt(-2*log(1-mu[1]))*sign(Y[1]-mu[1])  # shortcut, since y_1=0

# manually calculating the 1st working residual
resid(fit, type="working")[1]
(Y[1]-mu[1]) / (mu[1]*(1-mu[1]))

# manually calculating the 1st partial residual
resid(fit, type="partial")[1,1]
(Y[1]-mu[1]) / (mu[1]*(1-mu[1])) + fit$coefficients[2]*(x1[1] - mean(x1))
resid(fit, type="partial")[1,2]
(Y[1]-mu[1]) / (mu[1]*(1-mu[1])) + fit$coefficients[3]*(x2[1] - mean(x2))

# manually calculating the 1st quantile residual
library(statmod)
qresid(fit)[1] # results are random (uniformly), so won't come the same
a = pbinom(Y[1]-1, 1, mu[1]) 
b = pbinom(Y[1], 1, mu[1])
qnorm(runif(1, a, b)) # results are random (uniformly), so won't come the same
n = 10000
mean(replicate(n, qresid(fit)[1]))
mean(qnorm(runif(1000, a, b))) # should be close

For more information I suggest you check this book: Generalized Linear Models With Examples in R: working response - section 6.3, working residuals - section 6.7, response residuals - section 8.3.1, pearson residuals - section 8.3.2, deviance residuals - section 8.3.3, partial residuals - section 8.7.3

So,

will sum of squared residuals provide a meaningful measure of model fit ?

For Deviance/Pearson - I think so.

But more generally inspecting the residuals can be a bit tricky. In many cases neither the Pearson nor deviance residuals can be guaranteed to have distributions close to normal, especially for discrete distributions. "Small dispersion asymptotics" need to hold (see section 7.5 in the book), so some rule of thumbs are used. For Binomial distributions, and the deviance residual $\min(n_i y_i) > 3$ as well as $\min(n_i(1-y_i)) > 3$. There are also the Quantile Residuals that can be used when these conditions are not met. Check section 8.3.4 of the book.

$\endgroup$
6
  • 1
    $\begingroup$ In my book Regression Modeling Strategies I downplay the use of residuals in logistic regression because (1) logistic regression makes no distributional assumptions and (2) there are more direct ways to not only assess model fit but to make the fit more flexible in the needed directions through the use of splines and interactions. $\endgroup$ Commented Sep 5, 2020 at 12:05
  • 1
    $\begingroup$ @FrankHarrell I saw you also wrote somewhere as a comment to me on this site that small dispersion asymptotics is a myth. Unfortunately I'm not (yet) expert enough to decide on these matters :-). But I will try to check your book in the future. $\endgroup$ Commented Sep 5, 2020 at 12:22
  • $\begingroup$ @MaverickMeerkat above definition of error will also be okay with Fractional logistic fit? en.wikipedia.org/wiki/Fractional_model $\endgroup$ Commented Mar 25, 2022 at 16:41
  • $\begingroup$ There is a typo in the R code, it should be the upper case Y, fit = glm(Y ~ x1 + x2, family = 'binomial') $\endgroup$
    – Dan Li
    Commented Nov 30, 2023 at 18:17
  • $\begingroup$ Could you please explain a little bit what is the $t(y_i,y_i)$ and $t(y_i,\hat\mu_i)$ in the Deviance: $$sign(y_i-\hat\mu_i)*\sqrt{d_i}$$ where $d_i$ is the unit deviance, i.e. $d_i = 2(t(y_i,y_i)-t(y_i,\hat\mu_i))$ $\endgroup$
    – Dan Li
    Commented Nov 30, 2023 at 21:16
9
$\begingroup$

On Pearsons residuals,

The Pearson residual is the difference between the observed and estimated probabilities divided by the binomial standard deviation of the estimated probability. Therefore standardizing the residuals. For large samples the standardized residuals should have a normal distribution.

From Menard, Scott (2002). Applied logistic regression analysis, 2nd Edition. Thousand Oaks, CA: Sage Publications. Series: Quantitative Applications in the Social Sciences, No. 106. First ed., 1995. See Chapter 4.4

$\endgroup$
1
  • 9
    $\begingroup$ this is not entirely correct about large samples. It is rather that you require large binomial cell counts $n_i$, or what is the same thing, a large amount of replication of covariates. The pearson residuals are far from normally distributed for any observation where $n_i<5$. $\endgroup$ Commented Nov 24, 2011 at 9:26
6
$\begingroup$

The working residuals are the residuals in the final iteration of any iteratively weighted least squares method. I reckon that means the residuals when we think its the last iteration of our running of model. That can give rise to discussion that model running is an iterative exercise.

$\endgroup$
6
$\begingroup$
  1. Pearson Residuals

As tosonb1 points out, "The Pearson residual is the difference between the observed and estimated probabilities divided by the binomial standard deviation of the estimated probability".

I just wanted to mention that Pearson residual is mostly useful with grouped data i.e, say, there are $n_i$ trials at setting i of the explanatory variables (many observations for the same value of predictors) and let $y_i$ denote the number of “successes” for $n_i$ trials. Also, let $\hat{π_i} $ denote the estimated probability of success for the logistic regression model we have fit.

Pearson residual = $e_i = \frac{y_i - \hat{π_i}}{\sqrt {n_i \hat{π_i}(1 - \hat{π_i})}}$

For ungrouped binary data and often when explanatory variables are continuous, each $n_i$ = 1. Then, $yi$ can equal only 0 or 1, and a residual can assume only two values and is usually uninformative. Plots of residuals also then have limited use, consisting merely of two parallel lines of dots.

From, An Introduction to Categorical Data Analysis, 2nd Edition by Alan Agresti - vide chapter 5, section 5.2.4

  1. Deviance Residuals

(I am not entirely sure about this one, please point out errors, if any)

The i-th deviance residual can be computed as square root of twice the difference between loglikelihood of the ith observation in the saturated model and loglikelihood of the ith observation in the fitted model. Saturated Model is the model that predicts each observation perfectly and for all purposes in logistic regression, loglikelihood of saturated model $= ln(1) = 0$. Finally, we add a sign '+' in front of the residual if the observed response is 1 and put '-' if the observed response is 0.

Hence, deviance residual for the ith observation, $$d_i = (-1)^{y_i + 1}\sqrt{-2 (y_i ln(\hat{π_i}) + (1-y_i) ln(1 - \hat{π_i}))}$$ $y_i \in $ {0,1}

The sum of squares of deviance residuals add up to the residual deviance which is an indicator of model fit.

If a deviance residual is unusually large (which can be identified after plotting them) you might want to check if there was a mistake in labelling that data point.

$\endgroup$
1
  • $\begingroup$ Thank you! I was looking for an explicit formula of the deviance residuals and it was hard to find. This is exactly what I was looking for and is, indeed, the correct answer! $\endgroup$
    – Toby
    Commented Feb 29 at 18:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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