# Compare predicted versus actual outcomes in a GLM

I read somewhere that you could compute a "residual value" for a GLM by taking the actual values of your response variable divided by the predicted value of that response variable.

For example, suppose the response variable y represents number of cars, and $$x_1$$ represents the age of. a car. We would fit a glm model and calculate the "residual value", denoted residual below, for every observation in our data set with something like the following in R:

library(dplyr)
dataset <- dataset %>% mutate(pred_value = predict(m,type='response))
dataset %>% mutate(residual = y/pred_value)


I'm wondering if that actually makes sense, since unlike linear regression the GLM equation generally doesn't contain a residual term in it.

If not, what would be the best way to compare predicted versus expected values? The goal is the see if one can derive a 2nd predictor from the noise not modeled by the GLM model.

• This isn't a residual. A typical residual we use for GLM's is the deviance residual. You can get this using resid(model, type = 'deviance') – Demetri Pananos Apr 14 at 18:19
• Ah what I mean is that I want to – platypus17 Apr 14 at 19:19
• Ah what I mean by residual is a measure of the differences between the actual minus expected value of a response variable in a GLM, which is the same as a response residual in a linear regression context. So does the definition of a deviance residual measure the differences between the actual minus expected value of a response variable in a GLM? – platypus17 Apr 14 at 19:24
• In a sense, yes the deviance residual is kind of like that. – Demetri Pananos Apr 14 at 19:27
• I appreciate the response. Can you clarify why this is true a bit further? What is the explicit relationship between a deviance residual and how it measures the differences between the actual and expected values for a response variable? – platypus17 Apr 14 at 19:38

For a poisson glm, you can get the residuals, its similar to what @Demetri commented (y_observed - y_predicted) :

library(MASS)
data(Insurance)
dataset = Insurance
fit = glm(Claims ~ Age,offset=log(Holders),data=dataset)
dataset$$residuals = fit$$residuals


You can compare it against the outcome:

plot(fitted(fit),dataset\$residuals)


Now if we want to explore whether the residuals can be explained, I don't think you can fit a poisson glm again (might be wrong), so maybe we explore that with a regression tree:

library(randomForest)
f2 = randomForest(residuals ~ Group + Age + District,data=dataset,importance=TRUE)

importance(f2,type=1)
%IncMSE
Group    11.898473
Age      -1.058334
District 10.104154


You can see age has no more effect.. And in this dataset, unfortunately the other remaining factors have an effect. Maybe you can also check this

I am guessing your intend is to control / regress out the efficient of certain variables from your response, and fit them to another model. You can consider fitting a full model with all the covariates, regress out the so called "nuisance" parameters and fitting everything again.

I'm not going to go too much into the deviance goodness of fit test. Most books on GLM cover it quite thoroughly.

In short, the deviance goodness of fit test is a way to test your model against a so called saturated model; one which can perfectly predict the data. If the deviance between the saturated model and your model is not too large, then we can choose our model over the saturated model on the grounds that it is simpler and hence more parsimonious. See here for more.

The deviance goodness of fit statistic is comprised of individual contributions to the deviance. We can plot each of these contributions with the appropriate sign of $$y_i - \mathbb{E}(y\vert x)$$. Mathematically, the deviance goodness of fit statistic is $$D(\boldsymbol{y}, \hat{\boldsymbol{\mu}})=2 \sum_{i=1}^{N} \omega_{i}\left(y_{i}\left(\tilde{\theta}_{i}-\hat{\theta}_{i}\right)-\left(b\left(\tilde{\theta}_{i}\right)-b\left(\hat{\theta}_{i}\right)\right)\right)=\sum_{i=1}^{N} d_{i}$$

Where

$$d_{i}=2 \omega_{i}\left(y_{i}\left(\tilde{\theta}_{i}-\hat{\theta}_{i}\right)-\left(b\left(\tilde{\theta}_{i}\right)-b\left(\hat{\theta}_{i}\right)\right)\right)$$

and $$\tilde{\theta}_{i}$$ and $$\hat{\theta}_{i}$$ are estimates of the natural parameter under the saturated model and your model respectively.

Luckily, we don't actually have to compute those ourselves. R handles them using the command I mentioned in the comments. What R returns is soemthing a little different. R calculates

$$\operatorname{sign}\left(y_{i}-\hat{\mu}_{i}\right) \sqrt{d_{i}}$$

allowing residuals to be negative when the estimated outcome is larger than the observed data. We can plot these residuals against the predicted values and hope to see a cloud of them centered around the y=0 line, much like a linear model.