Regressing Logistic Regression Residuals on other Regressors

With OLS regression applied to continuous response, one can build up the multiple regression equation by sequentially running regressions of the residuals on each covariate. My question is, is there a way to do this with logistic regression via logistic regression residuals?

That is, if I want to estimate $$\Pr(Y = 1 | x, z)$$ using the standard generalized linear modeling approach, is there a way to run logistic regression against $$x$$ and get pseudo-residuals $$R_1$$, then regress $$R_1$$ on $$z$$ to get an unbiased estimator of the logistic regression coefficients. References to textbooks or literature would be appreciated.

• My guess is that this isn't going to work for the same reason that REML doesn't extend to GLMs; the magic of least squares is lost. I wonder if it would work in a fully bayesian context where you sampled the latent variable as part of the sim. The reason I wanted to do this was so I could run glmnet over different classes of variables and get different amounts of regularization for the classes - of course there are other ways to get this effect. – Ben Ogorek Sep 25 '18 at 12:39
• Is this very similar to using a back-fitting algorithm for logistic regression? – usεr11852 Sep 28 '18 at 20:40
• I mentioned this in a comment below, but in many implementations you can pass a 'base' prediction (offset parameter in glmnet), so maybe this would be possible after regressing the dependent vars. @BenOgorek do you want to add the purpose in the main text – seanv507 Oct 2 '18 at 18:38
• @seanv507 I'm worried that adding in the regularization part would increase the scope too much, especially now that there are some good answers below. After this Q&A wraps up I'll create a separate question where offset might indeed be our friend. – Ben Ogorek Oct 3 '18 at 0:23
• This is not an answer but I don't have enough reputation to comment. The question is about regressing the residual on the other regressor (i.e., predictors) rather than regressing residual on residuals. I am confused by the answers. – T Wu Mar 20 at 3:51

In standard multiple linear regression, the ability to fit ordinary-least-squares (OLS) estimates in two-steps comes from the Frisch–Waugh–Lovell theorem. This theorem shows that the estimate of a coefficient for a particular predictor in a multiple linear model is equal to the estimate obtained by regressing the response residuals (residuals from a regression of the response variable against the other explanatory variables) against the predictor residuals (residuals from a regression of the predictor variable against the other explanatory variables). Evidently, you are seeking an analogy to this theorem that can be used in a logistic regression model.

For this question, it is helpful to recall the latent-variable characterisation of logistic regression:

$$Y_i = \mathbb{I}(Y_i^* > 0) \quad \quad \quad Y_i^* = \beta_0 + \beta_X x_i + \beta_Z z_i + \varepsilon_i \quad \quad \quad \varepsilon_i \sim \text{IID Logistic}(0,1).$$

In this characterisation of the model, the latent response variable $$Y_i^*$$ is unobservable, and instead we observe the indicator $$Y_i$$ which tells us whether or not the latent response is positive. This form of the model looks similar to multiple linear regression, except that we use a slightly different error distribution (the logistic distribution instead of the normal distribution), and more importantly, we only observe an indicator showing whether or not the latent response is positive.

This creates an issue for any attempt to create a two-step fit of the model. This Frisch-Waugh-Lovell theorem hinges on the ability to obtain intermediate residuals for the response and predictor of interest, taken against the other explanatory variables. In the present case, we can only obtain residuals from a "categorised" response variable. Creating a two-step fitting process for logistic regression would require you to use response residuals from this categorised response variable, without access to the underlying latent response. This seems to me like a major hurdle, and while it does not prove impossibility, it seems unlikely to be possible to fit the model in two steps.

Below I will give you an account of what would be required to find a two-step process to fit a logistic regression. I am not sure if there is a solution to this problem, or if there is a proof of impossibility, but the material here should get you some way towards understanding what is required.

What would a two-step logistic regression fit look like? Suppose we want to construct a two-step fit for a logistic regression model where the parameters are estimated via maximum-likelihood estimation at each step. We want the process to involve an intermediate step that fits the following two models:

$$\begin{matrix} Y_i = \mathbb{I}(Y_i^{**} > 0) & & & Y_i^{**} = \alpha_0 + \alpha_X x_i + \tau_i & & & \tau_i \sim \text{IID Logistic}(0,1), \\[6pt] & & & \text{ } \text{ } Z_i = \gamma_0 + \gamma_X x_i + \delta_i & & & \delta_i \sim \text{IID } g. \quad \quad \quad \quad \quad \\ \end{matrix}$$

We estimate the coefficients of these models (via MLEs) and this yields intermediate fitted values $$\hat{\alpha}_0, \hat{\alpha}_X, \hat{\gamma}_0, \hat{\gamma}_X$$. Then in the second step we fit the model:

$$Y_i = \text{logistic}(\hat{\alpha}_0 + \hat{\alpha}_1 x_i) + \beta_Z (z_i - \hat{\gamma}_0 - \hat{\gamma}_X x_i) + \epsilon_i \quad \quad \quad \epsilon_i \sim \text{IID } f.$$

As specified, the procedure has a lot of fixed elements, but the density functions $$g$$ and $$f$$ in these steps are left unspecified (though they should be zero-mean distributions that do not depend on the data). To obtain a two-step fitting method under these constraints we need to choose $$g$$ and $$f$$ to ensure that the MLE for $$\beta_Z$$ in this two-step model-fit algorithm is the same as the MLE obtained from the one-step logistic regression model above.

To see if this is possible, we first write all the estimated parameters from the first step:

\begin{aligned} \ell_{\mathbf{y}| \mathbf{x}} (\hat{\alpha}_0, \hat{\alpha}_X) &= \underset{\alpha_0, \alpha_X}{\max} \sum_{i=1}^n \ln \text{Bern}(y_i | \text{logistic}(\alpha_0 + \alpha_X x_i)), \\[10pt] \ell_{\mathbf{z}| \mathbf{x}} (\hat{\gamma}_0, \hat{\gamma}_X) &= \underset{\gamma_0, \gamma_X}{\max} \sum_{i=1}^n \ln g( z_i - \gamma_0 - \gamma_X x_i ). \end{aligned}

Let $$\epsilon_i = y_i - \text{logistic}(\hat{\alpha}_0 - \hat{\alpha}_1 x_i) + \beta_Z (z_i - \hat{\gamma}_0 - \hat{\gamma}_X x_i)$$ so that the log-likelihood function for the second step is:

$$\ell_{\mathbf{y}|\mathbf{z}|\mathbf{x}}(\beta_Z) = \sum_{i=1}^n \ln f(y_i - \text{logistic}(\hat{\alpha}_0 - \hat{\alpha}_1 x_i) + \beta_Z (z_i - \hat{\gamma}_0 - \hat{\gamma}_X x_i)).$$

We require that the maximising value of this function is the MLE of the multiple logistic regression model. In other words, we require:

$$\underset{\beta_X}{\text{arg max }} \ell_{\mathbf{y}|\mathbf{z}|\mathbf{x}}(\beta_Z) = \underset{\beta_X}{\text{arg max }} \underset{\beta_0, \beta_Z}{\max} \sum_{i=1}^n \ln \text{Bern}(y_i | \text{logistic}(\beta_0 + \beta_X x_i + \beta_Z z_i)).$$

I leave it to others to determine if there is a solution to this problem, or a proof of no solution. I suspect that the "categorisation" of the latent response variable in a logistic regression will make it impossible to find a two-step process.

• Hi @Ben, thanks for teaching me about the Frisch–Waugh–Lovell theorem. I blew it on the bounty - thought "expired" meant it just stopped being advertised. Sorry about that. I like your likelihood based idea. Might try it out or something similar and post below. – Ben Ogorek Oct 6 '18 at 16:05
• @Ben Ogorek: No worries on the bounty. Glad the answer helped. – Ben Oct 6 '18 at 23:07
• @Ben Ogorek: (To make up for the lost 25 points of bounty, which vanish into the ether, just go around the site and up-vote any 3 answers. Then your karma is restored!) – Ben Oct 6 '18 at 23:11
• Done! (And I did read them first). – Ben Ogorek Oct 7 '18 at 15:07

I may be misinterpreting the question. I doubt you can build up the linear regression equation by regression on residuals in the way OP specified. OP's method would only work if the predictors are independent of each other.

To make it work, assume $$y$$ is the outcome vector, $$X$$ is the model matrix for the predictors already in the model and you want to include $$x_1$$. You need to regress the residual of the regression of $$y$$ on $$X$$ against the residual of the regression of $$x_1$$ on $$X$$ to obtain the OLS coefficient for $$x_1$$.

Here's a simple example:

set.seed(12345)
n <- 5000
x1 <- rnorm(n)
x2 <- .5 * x1 + rnorm(n) # Correlated predictors
y <- x1 + x2 + rnorm(n)


Fit model with OLS:

coef(lm(y ~ x1 + x2))
(Intercept)          x1          x2
0.001653707 1.037426007 0.996259446


Regression on residuals:

coef(lm(residuals(lm(y ~ x1)) ~ x2))
(Intercept)          x2
0.001219232 0.818774874


This is wrong, you need to fit:

coef(lm(residuals(lm(y ~ x1)) ~ residuals(lm(x2 ~ x1))))
(Intercept) residuals(lm(x2 ~ x1))
-6.707350e-17           9.962594e-01


Which returns the right coefficient for x2, this aligns with expected differences in y given differences in x2, holding x1 constant (taking it out of both y and x1).

That aside, in logistic regression, it would even be more problematic because logistic regression coefficients suffer from omitted variable bias even in the absence of confounded relations, see here and here, so unless all predictors of the outcome are in the model, one cannot obtain unbiased estimates of the true population parameters. Moreover, I do not know of any residuals from the model that would be amenable to a second logistic regression with all values lying between 0 and 1.

Some references on regression on residuals:

• Maxwell, S. E., Delaney, H. D., & Manheimer, J. M. (1985). Anova of Residuals and Ancova: Correcting an Illusion by Using Model Comparisons and Graphs. Journal of Educational Statistics, 10(3), 197–209. Retrieved from http://journals.sagepub.com/doi/pdf/10.3102/10769986010003197
• Freckleton, R. P. (2002), On the misuse of residuals in ecology: regression of residuals vs. multiple regression. Journal of Animal Ecology, 71, 542-545. doi:10.1046/j.1365-2656.2002.00618.x
• I think your first couple of paragraphs are slightly misleading/unclear...it would be better if you started with how you actually do 'linear regression with residuals'..(+1) and you can find it in elements of statistical learning (multiple regression from single regressions subsection?) – seanv507 Oct 2 '18 at 7:48
• In many implementations you can pass a 'base' prediction (offset parameter in glmnet), so maybe this would be possible after regressing the dependent vars – seanv507 Oct 2 '18 at 8:21
• @seanv507 I already include it in my answer. It's the last code demonstration I have. It's just not possible in the way OP described, regressing residuals on a predictor. But I could rewrite it to show the proper way from the start if that's what you mean. – Heteroskedastic Jim Oct 2 '18 at 10:22
• Yes I meant rewrite it to show the proper way from the start, – seanv507 Oct 2 '18 at 18:36
• @seanv507 don't know what you mean by you can pass a base prediction? And regressing the dependent variables? – Heteroskedastic Jim Oct 3 '18 at 2:03

I hope I am not misinterpreting your question, as my answer is going to change somewhat the wording of how you phrased your subject.

I think what you are trying to do is build your regression model by adding one independent variable at a time. And, you do that by observing which prospective variable has the highest correlation with the residual of your first regression between Y and X1. So, the variable with the highest correlation with this first residual will be X2. So, now you have a model with two independent variables X1 & X2. And, you continue this exact process to select X3, X4, etc. This is a stepwise forward process.

You can do the exact same thing with Logistic Regression for the simple reason that Logistic Regression is pretty much an OLS Regression where the dependent variable is the log of the odd (or logit). But, whether Y is a logit or not does not affect the stepwise forward process mentioned above.

OLS minimizes the sum of the square errors to fit the actual data. Logit regression uses a maximum likelihood process that generates a fit that is not all that different than OLS. And, that too (the fitting mechanism) should not affect the stepwise forward process that allows you to build your multiple regression model, whether the latter is an OLS Regression or a Logit Regression.