# Why aren't residuals in Poisson regression on the scale of the response variable?

I am running models in R using an OLS as well as using a GLM with a Poisson distribution and log link function.

df <- data.frame(x1=runif(100), x2=runif(100))

df$$ols_y <- rnorm(100) df$$glm_y <- rpois(100, lambda=5)

ols_mod <- lm(ols_y ~ x1 + x2, data=df)
glm_mod <- glm(glm_y ~ x1 + x2, data=df, family=poisson)


Now, I can recover the original y values if I use the ols model to predict the training data and add back the residual (rounding to avoid floating-point issues). I can also use fitted.values instead of re-running the model on the training data:

all(round(ols_mod$fitted.values + ols_mod$residuals, 5) == round(df$ols_y, 5)) > TRUE all(round(predict(ols_mod, df) + ols_mod$residuals, 5) == round(df$ols_y, 5)) > TRUE  But if I do it with my Poisson model, even using type="response", the values don't match. all(round(predict(glm_mod, df, type="response") + glm_mod$residuals, 5) == round(df$glm_y, 5)) > FALSE all(round(glm_mod$fitted.values + glm_mod$residuals, 5) == round(df$glm_y, 5))
> FALSE


Assuming this was due to the log-link, I also tried to exponentiate the residuals to get them on the scale of the linear predictor, but this also fails.

all(round(predict(glm_mod, df, type="response") + exp(glm_mod$residuals), 5) == round(df$glm_y, 5))


What am I missing here? Why aren't the residuals on the scale of the response variable? Why isn't $$\hat{y}_i + r_i = y_i$$ ?

• Take a look at the help page at ?residuals.glm. What you are looking at are deviance residuals Dec 11, 2023 at 10:08
• @Knarpie, the glm\$residuals object contains working residuals. It's still correct that you need residuals(glm_mod, type="response"). Dec 11, 2023 at 10:14
• Mathematical error: exponentiation couldn't possibly retrieve negative residuals, regardless of what you feed it. Dec 11, 2023 at 11:24
• What happens if you exponentiate the absolute value, somehow keeping track of which were negative? Dec 11, 2023 at 11:59