Logistic Regression - Getting Pearson Standardized Residuals in R vs Stata

Question

I am working on an assignment involving a logistic regression model, where I need to plot the pearson standardized residuals against one of the predictors. Here's the basic setup:

model <- glm(outcome ~ predictor1 + predictor2, family=binomial(logit))
res <- residuals(model, "pearson")

When looking at the residuals' distribution, I see something totally different than my colleagues who use Stata (using predict and rstandard). Their residuals are more or less normal, whereas in mine there is a gap in the values (not a singe residual is between -0.05 and 1.15). That does make sense in the context of logistic regression, especially that the maximum predicted probability is not so high (38%).

I'd like to understand what's happening here... What is Stata doing that R isn't, with those residuals?

Answer 1

For logistic regression, Stata defines residuals and related quantities to be those you'd get if you grouped all the observations with the same values for all the predictor variables, counted up the successes and failures for those observations, and fitted a logistic regression model to the resulting binomial data instead of the original Bernoulli data. This is a useful thing to do as (if there are multiple observations with the same covariate pattern) the resulting residuals behave more like those you're used to from least squares.

To get the same residuals from R, I suspect you will need to group the data and fit the model to the grouped data. But I'm not clear whether R is using the same definition of 'standardized residuals' as Stata as I don't presently have access to the numerous textbooks that the R documentation references.

Here's an excerpt from 'Methods and formulas' section of the Stata manual entry for 'logistic postestimation' (one thing I like about Stata is that the manuals provide the full formulas for everything):

Define $M_j$ for each observation as the total number of observations sharing $j$’s covariate pattern. Define $Y_j$ as the total number of positive responses among observations sharing $j$’s covariate pattern.

The Pearson residual for the $j$th observation is defined as $$r_j = \frac{Y_j - M_j p_j}{\sqrt{M_j p_j(1 - p_j)}}$$ ...
The unadjusted diagonal elements of the hat matrix $h_{Uj}$ are given by $h_{Uj} = (\mathbf{XVX}')_{jj}$, where $\mathbf{V}$ is the estimated covariance matrix of parameters. The adjusted diagonal elements $h_j$ created by hat are then $h_j = M_j p_j(1 - p_j)h_{Uj}$.
The standardized Pearson residual $r_{Sj}$ is $r_j / \sqrt{1 - h_j}.$

Answer 2

Pearson residuals are obtained by dividing the each observation's raw residual by the square root of the corresponding variance. The idea is to get something that has variance 1, approximately. In your example, try this;

set.seed(3141)
x1 <- rnorm(100)
x2 <- rnorm(100)
y <- rbinom(100, 1, 0.25)
glm1 <- glm(y~x1+x2, family=binomial)
f1 <- fitted(glm1) # the fitted probability of y=1, for each observation
plot( residuals(glm1, "pearson"), (y-f1)/sqrt(f1*(1-f1)))
abline(0,1)        # they match

The 'gap' occurs because the residuals where $Y=1$ are on one side, and those with $Y=0$ are on the other. Standardized residuals are a different animal; they divide by the estimated standard deviation of the residual; you can obtain them in R using rstandard(), though for non-linear GLMs it uses a linear approximation in the calculation.

NB residuals of any form tend not to be terribly helpful in logistic regression. With independent binary data, the only real concern is whether we've specified the mean correctly - and with modest sample sizes, plots of residuals typically provide little power to assess that.