# Why GLM don't have an error term and why shouldn't residuals be i.i.d?

I've read dozens on post on the subject but I cannot figure this out. From what I've gathered, GLMS don't include an error term in their formulation unlike linear models (LM). I was wondering why (or why do LM include that error term if they're the odd ones) ? I feel like this is due to the additive properties of the normal distribution which makes it possible in LM to model the distribution of Y by focusing on its mean and adding a normally distributed variable centered on 0 (the error term) with a variance that fits the gap between the variance of $$\sum_i X_i$$ and the variance of Y. But isn't the link function supposed to turn non-normal distributed variables into normally distributed ones and therefore make it possible to use an error term ?

Note : if you can provide examples with continuous distributions in your answers and not just discrete ones it would be awesome because it's harder for me to understand this problem with discrete distributions.

Thanks

GLMs specify a conditional distribution at each combination of predictors.

When that conditional distribution is normal with constant variance, then the errors about the mean will be i.i.d. normal and the marginal distribution of residuals will (usually) be approximately so.

Outsider that special case, those conditional distributions don't have the same spread and typically don't even have the same shape. So when you subtract the collection of conditional "population" (process) means, whenever the means differed, the resulting error terms have a different distribution.

For a first continuous example consider a gamma GLM. The shape of each conditional distribution is the same but the spreads differ. If you subtract the mean, the support of each is different so the "errors" differ.

See the displays here: https://stats.stackexchange.com/a/224253/805

which contrast an additive model with a particular gamma-distributed error (one with shape=1) against a GLM with a response is conditionally gamma distributed.

With the inverse Gaussian it's even worse - the shape changes as the mean changes.

However, perhaps best simple example is a Poisson GLM (I know you prefer not to have discrete examples but it's still just about ideal as a base case and you're not the only person who this explanation might help). The spread and shape of each Poisson with a different mean is different. When you subtract the means, you have a collection of different distributions, with different support. If you scale for the difference in spreads, they still have different support and different shapes.

Your post seems to be assuming that a link function is used as a data transformation. If that's the way you were thinking about it, that premise is false. Even if you did use it as a transformation of the response, it doesn't produce normality in general, not even approximately. For natural link GLMs the linearizing transform is different from the variance-stabilizing transform which is in turn different from the symmetrizing transform; indeed within the exponential dispersion family these transforms get progressively weaker.

• Thanks for your answer. Digging into this makes me realise there is soo much I don't understand or have a crooked understading of, even for simple linear regression. Your answer does clarify a lot of stuff though. I'm not familiar with all the transfromations you mention in your last paragraph and I realise "transformation" might have a specific technical sense that I'm not aware of. I do have one very basic question however. What is the purpose of a link function ? In my mind, the linear model is useful because of the additive proprieties of the normal distribution (rest in the next comment). Commented May 20 at 20:54
• One can add the effect of different variables (X) in different proportions (the coefficients) and reproduce the distribution of the observed variable Y. However that does not work with other distributions. I thought that using the link function was a way to solve that. For example, you take the logit of odds between 0 and 1 and you get normally distributed values, same if you take the log of count data. And then you could model that new, somewaht normally distributed variable, with a linear model, aka a sum of normally distributed variables. What am I missing ? Commented May 20 at 20:54
• "For example, you take the logit of odds between 0 and 1 and you get normally distributed values, same if you take the log of count data." -- You should explain what led you to think these two things would be the case. Neither of those claims is true, and I already directly said that the premise they arise from is false in my answer, in the first two sentences of my last paragraph., while the final sentence is an explanation of why. I'll think about how to explain it in more detail. Commented May 20 at 23:10
• @Boussens-DumonGrégoire The relationship is whatever it is. The link function allows a particular type of nonlinearity, however, to enter the model. For instance, if a logistic regression model is correct, the relationship between the features and conditional expectation is nonlinear. However, the relationship between the features and the transformed conditional expectation is linear.
– Dave
Commented May 21 at 15:21
• @Glen_b I guess I understood it wrong when it was taught to me. To be honest I don't understand anything about your last paragraph. I'll look on the Internet what all those transformations are but right now I don't know what they mean nor how one is weaker or stronger than the other. I was reading the wikipedia page for GLMs and here is another try. Is the purpose of the link function only to transform a relation that is not linear (the relation between E(Y) and the predictor) into a linear relation and satisfy the assumptions of the distribution assumed for E(Y|X) ? Commented May 21 at 18:48

GLMs are not supposed to generalize $$y = X\beta + \varepsilon$$. They generalize the equivalent$$^{\dagger}$$ notion of $$\mathbb E\left[Y\vert X=x\right] = X\beta$$. There is no error term in this about which we can make assumptions.

The generalization comes from noticing that the previous equation is equivalent to $$g\left(\mathbb E\left[Y\vert X=x\right]\right) = X\beta$$ when $$g$$ is the identify function.

So GLMs do not include an error term, or at least are not typically written in terms of an error term, because that deviates from how they generalize linear models. They generalize linear models by starting at a form of linear models that does not include an error term.

Sure, we can write the model in terms of an expected value plus by how much the expected value misses the observed value, but then the connection to a decent distribution is lost. When you use a Gamma distribution, the expected values can be higher or lower than the observed values. Therefore, the residuals are some weird mixture distribution of both positive and negative values, losing sight of the conditional Gamma distribution.

Then how do you estimate the model parameters? In linear models, we might choose to use OLS, which was developed before the discovery of Gaussian maximum likelihood estimation and does not require conditional Gaussian distributions (though the equivalence with maximum likelihood estimation is nice). Do we go with minimization of square loss in a generalized linear model? That is an option, but then we would want to know what properties that estimator has. Under what, if any, conditions is that estimator unbiased? Under what, if any, conditions is that estimator consistent? Are we operating under those conditions? How efficient is that estimator? How robust is that estimator? How do we write confidence intervals for the regression coefficients or test regression coefficients?

But there is a Gamma distribution in there. Thus, if we really believe in that Gamma likelihood,$$^{\ddagger}$$, we have an obvious estimator through maximum likelihood estimation. We know many theorems about maximum likelihood estimators and how great their properties are. We have theorems about how to calculate parameter confidence intervals. We have theorems about how to test nested models.

So the options are:

1. Shoehorn the GLM into a format that explicitly uses residuals with a nasty distribution that does not lead to an obvious estimation strategy with known, desirable properties, nor are there clear ways to calculate confidence intervals or test nested models.

2. Use the format that does not use residuals, leading to the usual types of point estimation, interval estimation, and hypothesis testing that we know and love from classical linear regression (OLS with a Gaussian assumption, equivalent to maximum likelihood estimation).

The field of statistics, I think understandably, has chosen to go with the second option instead of working out GLM theory in terms of residuals.

$$^{\dagger}$$Equivalent under the usual assumptions, such as Gauss-Markov (though even Gauss-Markov is more restrictive than is required)

$$^{\ddagger}$$This could be a dubious assumption for the Gamma. However, consider a logistic regression where the outcome is binary. The conditional distribution is Bernoulli.

• Hi @Dave. Thanks for your answer. I've been reading posts about that all day and my head is about to explode but I think I't slowly starting to understand. One things I'm still uncertain about. Let's say the distribution of $Y|X$ is a Gamma. Then what is the distribution of $g(Y|X)$ ? Commented May 20 at 22:55
• @Boussens-DumonGrégoire The point of the link function is to transform the (unobserved) expected value, not the entire distribution.
– Dave
Commented May 20 at 23:36
• And the goal of this transformation is to get a linear relation between the predictor and Y right ? Commented May 21 at 1:03
• @Boussens-DumonGrégoire The expected value
– Dave
Commented May 21 at 1:10
• @Boussens-DumonGrégoire I think Hardy's comment was directed to me, and I edited in response.
– Dave
Commented May 21 at 19:06

Arguably, regression concerns modeling the expected response. If we focus on this aspect of modeling, the "error" is less something intrinsic to a model than it is a way of calculating the uncertainty of our estimates. You seem to already sense that the classic assumption of OLS - that the error is independent of the fitted values - is not attainable in many settings.

For instance, in a binary response, no matter what transformation, the distribution of the "error" defined as the observed value minus the predicted value will always be dependent on that predicted value. For $$p$$ (true probability) close to 0.5, you would find the error to be symmetric, whereas $$p$$ close to 1 would have negatively skewed errors, and $$p$$ close to 0 the opposite regardless of the scale (normal, logarithmic, logistic, etc.)

In fact, if you followed the approaches of handling heteroscedasticity by inverse variance weighting, you would derive the Newton-Raphson method to solve the likelihood equations as a kind of method of moments estimator - thus reinforcing the idea this is actually about modeling expectations.

In the OLS framework you can write two models equivalently:

$$E[Y|X] = \alpha + \beta X$$

or

$$Y = \alpha + \beta X + \epsilon$$

Neither really reveals a prudent estimation routine, nor does either suggest any different approach unless we stipulate additional information on $$\epsilon$$.

• I think I understand your example for the binomial distribution. You want fitted values to be either 0 or 1. Therefore if the observed value is 1 you error will be 0 or negative and if it is 0, it will be 0 or positive. Hence the distribution of the error depends on the fitted value. Could you give me an example why we have the same problem with continuous distributions such as Gamma distributions for example ? Commented May 20 at 17:38
• Also, when you use the logit kink function, you go from a discrete distribution contained within [0-1] to a continuous unbounded distribution. It you try to fit data on the logit scale (estimate $\text{logit}(y_i)$), your errors are not bounded either and therfore couldn't you have an error term in your model which would be iid ? Commented May 20 at 17:45
• @Boussens-DumonGrégoire to your second comment, you keep speaking of "distributions". What distribution are you referring to? Note, if I take any $\text{Binomial}(n,p)$ distribution, then subtract $p$ from it, that distribution is still discrete. In fact, it is quite possible to guess both the $n$ and $p$ depending on the sample size. Commented May 20 at 20:24

GLMs (generalized linear models) include a lot of different types of models. It's more of a family than a single model. Notably, the typical linear model is a special case of the GLM. A GLM consists of a specified response distribution from the exponential family (which certainly could be normal), and a link function. Let's take a different special case of the GLM, namely logistic regression. In that case, the response is assumed to be conditionally binomial and the link is the logit (i.e., the log of the odds of 'success'). Now recognize that the variance of the residuals of a binomial are dependent on the mean. Thus, you have heteroscedasticity and non-normality. So, surely, there is at least one form of the GLM that does not have those features (actually, there are many more).

Response to edit:
The link function certainly does not turn the non-normal distribution into a normally distributed response, nor is it 'supposed' to. Again, it is constructive to consider the binomial case. Imagine a binary outcome (e.g., cancer / not) regressed onto a binary exposure (e.g., radon / not). The logit transformed $$0$$ (no cancer) is $$-\infty$$, and the logit transformed $$1$$ (cancer) is $$\infty$$. They are not normally distributed. The variance of a Bernouli is: $${\rm Var}(X) = p(1-p)$$ where $$p$$ is the probability of 'success' (i.e., the mean).

• Thanks for the edit. There's still a lot of things I'm confused about haha. The logit of 0 and 1 are clearly not normal indeed but I thought in binary regression the logit function was applied on the odds and not directly on the 0s and 1s. And when you generate a vector of numbers between 0 and 1 and apply the logit function to these numbers, the distribution of the values you get is actually normal. So by applying the inverse-logit function to a normally distributed predictor you should be able to get estimated odds from which you can derive estimated values for Y right ? Commented May 20 at 18:02
• @Boussens-DumonGrégoire, consider the binary Y, binary X example I made up. The 2 observed proportions are $.2$ and $.05$. The logits are $-1.386294$ and $-2.944439$. They are not normally distributed. Commented May 20 at 18:05

#### You can give a GLM an "error term" if you want to --- it is just not useful

To understand this issue, it is worth going back to linear regression and understanding the definition of the error term in that context. In the linear regression model we define the error vector to be the deviation of the observed response vector from the true regression function (i.e., from the true conditional expectation of the response vector given the explanatory variables):

$$\boldsymbol{\varepsilon} \equiv \mathbf{Y} - \mathbb{E}(\mathbf{Y}|\mathbf{x}).$$

By defining the error vector in this way, we can then decompose the response vector in the usual form:

$$\mathbf{Y} = \mathbb{E}(\mathbf{Y}|\mathbf{x}) + \boldsymbol{\varepsilon}.$$

Now, in a standard Gaussian linear regression we have a linear regression function $$\mathbb{E}(\mathbf{Y}|\mathbf{x}) = \mathbf{x} \boldsymbol{\beta}$$ and we have the error distribution $$\boldsymbol{\varepsilon} = \mathbf{Y} - \mathbb{E}(\mathbf{Y}|\mathbf{x}) \sim \text{N}(\mathbf{0}, \sigma^2 \mathbf{I})$$, which gives the simplified model form:

$$\mathbf{Y} = \mathbf{x} \boldsymbol{\beta} + \boldsymbol{\varepsilon} \quad \quad \quad \quad \quad \boldsymbol{\varepsilon} \sim \text{N}(\mathbf{0}, \sigma^2 \mathbf{I}).$$

This is a useful decomposition because it splits the response vector into the sum of a part attributable to the true regression function and a deviation from this function (which we call the "error vector") that has zero expected value. That decomposition is simple, intuitive and appealing.

Now, once you go into the territory of GLMs, you can do the exact same thing, but it is just not particularly useful. As with a standard regression model, a GLM has a true regression function that measures the conditional expected value of the response given the explanatory vector. In a GLM, this has the form $$\mathbb{E}(\mathbf{Y}|\mathbf{x}) = g^{-1}(\mathbf{x} \boldsymbol{\beta})$$ where $$g$$ is some link function. If you define the "error vector" as before, you get:

$$\mathbf{Y} = g^{-1}(\mathbf{x} \boldsymbol{\beta}) + \boldsymbol{\varepsilon} \quad \quad \quad \quad \quad \boldsymbol{\varepsilon} \sim F,$$

where the distribution $$F$$ is now usually some "shifted" distribution with non-zero mean. For example, if we are using logistic regression, the distribution of the error terms in this decomposition would be a shifted Bernoulli distribution where each error term $$\varepsilon_i$$ lies on the two-point support consisting of the points $$-g^{-1}(\mathbf{x} \boldsymbol{\beta})$$ and $$1-g^{-1}(\mathbf{x} \boldsymbol{\beta})$$. This is a rather strange distribution with a rather strange support, and most importantly, the support is dependent on the explanatory variables which means that each error term has a different support.

What this essentially gives you is a mathematically valid decomposition that is really just a big mess. It is not simple, intuitive or appealing, and so it doesn't really aid in understanding the model. The reason for this is that (unlike in Gaussian linear regression models) the different distributional structure in a GLM does not decompose nicely into a simple fixed distributional form for the deviation of the response variables from their true conditional expected value.