Deriving the Distributions of Residuals in Linear Regression

Question

This is a follow up to my question here: Deriving Distributions of Linear Regression

I am trying to manually derive the distribution of the observed residuals (I call the $\epsilon$ as the theoretical errors and $e$ as the observed residuals) in a linear regression model. I keep getting confused on how to do this and would appreciate some guidance on this.

Here is the theoretical model (I found this very difficult to understand, the conditional distributions are normally distributed - not the marginal distributions). I like to write these as conditional on $x_i$, $\beta$ - because $y_i$ and $\epsilon_i$ could not have existed without the existence of a regression model (i.e. involving $x_i$, $\beta$). This might be a very redundant way to write this, but I hope its correct:

$$y_i = x_i\beta + \epsilon_i$$ $$P(\epsilon_i \mid x_i, \beta) \sim N(0, \sigma^2)$$ $$P(y_i \mid x_i, \beta) \sim N(X_i\beta, \sigma^2)$$

Here is the applied model:

$$\hat{y_i} = x_i\hat{\beta}$$ $$e_i = y_i - \hat{y_i} = x_i\beta + \epsilon_i - x_i \hat{\beta} = x_i(\beta-\hat{\beta}) + \epsilon_i$$

• Using the Law of Total Expectations (https://en.wikipedia.org/wiki/Law_of_total_expectation), we can write:

$$E(\epsilon_i) = E[E(\epsilon_i \mid x_i)] = E[0] = 0$$

• We also know that $\hat{\beta}$ is Unbiased (using the fact that $E(\epsilon)=0$):

$$\hat{\beta} = (X^TX)^{-1}X^Ty$$ $$ \text{substitute } y = X\hat{\beta} \text{, then: } \hat{\beta} = (X^TX)^{-1}X^T(X\beta + \epsilon)$$

$$ E(\hat{\beta}) = E((X^TX)^{-1}X^T(X\beta + \epsilon)) = E[(X^TX)^{-1}X^TX\beta + (X^TX)^{-1}X^T\epsilon] = E[I\beta + (X^TX)^{-1}X^T\epsilon] = \beta + E[(X^TX)^{-1}X^T\epsilon] = \beta + (X^TX)^{-1}X^TE(\epsilon) = \beta$$

• Again, using this unbiased result [i.e. $\hat{\beta} = \beta$ and $E(\hat{\beta} - \beta) = 0]$ as well as the result from the Total Law of Expectations [i.e. $E(\epsilon)=0$], we can write:

$$E[e \mid x] = E[x(\beta - \hat{\beta}) + \epsilon] = E[x(\beta - \hat{\beta})] + E(\epsilon) = x \cdot E[(\beta - \hat{\beta})] + E(\epsilon)=0$$

$$E(e) = E[E(e \mid x)] = E(0) = 0 $$

Ideally from here, I would try to find out the variance:

$$Var(e) = E(e^2) - [E(e)]^2 = E(e^2) - 0 = E(e^2)$$

However, here is where I run into a problem: I don't know how to calculate $E(e^2)$ and I also don't know how to determine if $e$ in general will have a Normal Distribution.

I tried to approach this a different way:

$$ e = y - \hat{y} = y - x \hat{\beta} = y - x (X^TX)^{-1}X^Ty $$

In the above expression $y - x (X^TX)^{-1}X^Ty$, if we factor $Y$ out, we can write:

$$y - x (X^TX)^{-1}X^Ty = [I - X (X^TX)^{-1}X^T]Y = MY$$ $$M = I - X (X^TX)^{-1}X^T$$

From here, we can also see that:

$$MX = IX - X (X^TX)^{-1}X^TX = IX - XI = 0$$

$M$ is also an idempotent matrix - this can be manually verified as well, i.e. $M M^T = M$.

Now, we substitute $Y = X\beta + \epsilon$ into the $MY$:

$$e=M(X\beta + \epsilon) = M(X\beta) + M\epsilon = M\epsilon $$

Doing this, we have basically transformed the residuals $e$ into a function of the theoretical errors $\epsilon$.

From here, we can see that Expected Value of $e$ is (using the Law of Total Expectations from above):

$$E(e) = E(M\epsilon) = M\cdot E(\epsilon) = M \cdot 0 = 0$$

We can also calculate the Variance of $e$ as (by invoking the properties of idempotent matrices):

$$Var(e) = Var(M\epsilon) = M \cdot Var (\epsilon) \cdot M^T = \sigma^2 \cdot M = \sigma^2 \cdot M$$

And finally, we can see that $e$ also has to be Normally Distributed, since $\epsilon$ is Normally Distributed and $e$ is simply $M$ (a constant) multiplied by a $\epsilon$.

Therefore, the distribution of the residuals $e$ is:

$$ P(e_i | x_i, \beta) \sim N(0, M \cdot \sigma^2) $$

Here, I think that $M \cdot \sigma^2$ is a placeholder for the estimate $\hat{ \sigma^2}$ (i.e. $M$ contains the influence of measurements and plays the role of a scaling/correction factor to adjust for estimation biases).

Is the correct derivation?

Note: An additional derivation directly from $e = M\cdot \epsilon$

We know that $E(e) = E(M\epsilon) = M\cdot E(\epsilon) = M \cdot 0 = 0$.

In general, $Var(e) = E(e^2) - [E(e)]^2$. Since $E(e)=0$, this means that $Var(e) = [E(e)]^2$.

Since we are dealing with matrices, $Var(e) = E(e \cdot e^T)$

Further manipulating $e = M\cdot \epsilon$, we can write $ee^T = (M\cdot \epsilon) ( M\cdot \epsilon)^T = M \epsilon \epsilon^T M$

Taking the expectation, we can write $Var(e) = E(ee^T) = M E(\epsilon \epsilon^T) M^T$

We have to remember that $Var(e)$ is a matrix of the following form:

$$e e^T = \begin{bmatrix} e_{1} \\ e_{2} \\ \vdots \\ e_{n} \end{bmatrix} \begin{bmatrix} e_{1} & e_{2} & \cdots & e_{n} \end{bmatrix} = \begin{bmatrix} e_{1}e_{1} & e_{1}e_{2} & \cdots & e_{1}e_{n} \\ e_{2}e_{1} & e_{2}e_{2} & \cdots & e_{2}e_{n} \\ \vdots & \vdots & \ddots & \vdots \\ e_{n}e_{1} & e_{n}e_{2} & \cdots & e_{n}e_{n} \end{bmatrix}$$

$$E(ee^T) = Var(e) = \begin{bmatrix} Var(e_1) & Cov(e_1, e_2) & \cdots & Cov(e_1, e_n) \\ Cov(e_2,e_1) & Var(e_2) & \cdots & Cov(e_2, e_n) \\ \vdots & \vdots & \ddots & \vdots \\ Cov(e_n, e_1) & Cov(e_n, e_2) & \cdots & Var(e_n) \end{bmatrix} = \begin{bmatrix} Var(e_{1}) & 0 & \cdots & 0 \\ 0 & Var(e_{2}) & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & Var(e_{n}) \end{bmatrix}$$

Using the same logic as we did for residuals, the variance of errors can be found out the same way:

$$Var(\epsilon) = E(\epsilon^2) - [E(\epsilon)]^2 = E(\epsilon^2) - 0 = E(\epsilon \epsilon^T)$$

$$\epsilon \epsilon^T = \begin{bmatrix} \epsilon_{1} \\ \epsilon_{2} \\ \vdots \\ \epsilon_{n} \end{bmatrix} \begin{bmatrix} \epsilon_{1} & \epsilon_{2} & \cdots & \epsilon_{n} \end{bmatrix} = \begin{bmatrix} \epsilon_{1}\epsilon_{1} & \epsilon_{1}\epsilon_{2} & \cdots & \epsilon_{1}\epsilon_{n} \\ \epsilon_{2}\epsilon_{1} & \epsilon_{2}\epsilon_{2} & \cdots & \epsilon_{2}\epsilon_{n} \\ \vdots & \vdots & \ddots & \vdots \\ \epsilon_{n}\epsilon_{1} & \epsilon_{n}\epsilon_{2} & \cdots & \epsilon_{n}\epsilon_{n} \end{bmatrix}$$

$$E( \epsilon \epsilon^T) = Var(\epsilon) = \begin{bmatrix} Var(\epsilon_1) & Cov(\epsilon_1, \epsilon_2) & \cdots & Cov(\epsilon_1, \epsilon_n) \\ Cov(\epsilon_2,\epsilon_1) & Var(\epsilon_2) & \cdots & Cov(\epsilon_2, \epsilon_n) \\ \vdots & \vdots & \ddots & \vdots \\ Cov(\epsilon_n, \epsilon_1) & Cov(\epsilon_n, \epsilon_2) & \cdots & Var(\epsilon_n) \end{bmatrix} = \begin{bmatrix} Var(\epsilon_{1}) & 0 & \cdots & 0 \\ 0 & Var(\epsilon_{2}) & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & Var(\epsilon_{n}) \end{bmatrix} = \begin{bmatrix} \sigma^2 & 0 & \cdots & 0 \\ 0 & \sigma^2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \sigma^2) \end{bmatrix} = \sigma^2 \cdot I $$

Now, going back to our original equation $Var(e) = E(ee^T) = M E(\epsilon \epsilon^T) M^T$, we can see that (note that $M M^T = M$):

$$Var(e) = M \sigma^2 I M^T = \sigma^2 M M^T = \sigma^2 M$$

A final point to be noted is that the diagonal elements in the matrix $M = I - X (X^TX)^{-1}X^T$ are not necessarily identical . Given that $e = M \epsilon$, this means that the variance of the residuals at different $x$ values could be different. This means that in real life, the residuals likely have a mixture of Normal Distributions.

Answer 1

For the vector of residuals $e$ we can write $e= M\epsilon=(I-H)\epsilon$ where $H$ is the so called "hat" matrix, being $X(X'X)^{-1}X'$. For the $N$ residuals in the sample, the $N*N$ (co)variance matrix $E(ee')$ is equal to: $E(ee') = \sigma^2M=\sigma^2(I-H)$. The diagonal terms of this matrix contain the variances of all individuals' residuals. This means that there is not just one single variance which holds for all cases in the sample, but that for each (combination of) $X$ value(s), there is a particular variance of the corresponding estimated residual. This variance given $X$ for case $i$ can be expressed as $\sigma^2(1-h_{ii})$ and thus depends on the case's diagonal element in the hat matrix, which in turn depends on the $X$ value(s), and on the value of $\sigma^2$ which is unknown, but can be estimated by $s^2=SSE/(N-K)$ in the usual way.

So, talking of just one variance of the estimated residuals $e$ is problematic somehow. The separate residuals $e$ given $X$ are all normally distributed, given that the true errors are. But if we would consider all estimated residuals together and wonder which one-and-only distribution they come from, that is another question, the answer of which is not so simple, unfortunately.

Further it seems relevant to me that the residuals often function as an instrument for evaluating the assumptions made in regression. Non-normal residuals can point to non-normal error terms and could be caused by e.g. omitted variables. With that in mind, we should be careful with stating that the residuals have a normal distribution. Only if the true errors have, the residuals will also have a "quite" normal distribution (given what was said before), but not as a rule.

It's illustrative to show things in an example. In the following R script I generate data in such a way that 2% of all cases are extreme in terms of their X values.

set.seed(12345)
library(summarytools)
library(ggplot2)

# Generate X data, 49% being -1, 49% being +1, 1% being 100, 1% being -100.
xmedium <- c(rep(-1,4900), rep(1,4900))
xlow    <- rep(-100, 100)
xhigh   <- rep(100, 100)
x <- c(xlow, xmedium, xhigh)

# Add a factor for the relative position of the X value.
xposition <- factor(c(rep(1,100), rep(2,4900), rep(3,4900), rep(4,100)))

# Generate the Y data with normal error.
y <- 0 + x + rnorm(length(x), 0, 10)

# Estimate the model and extract the residuals.
model <- lm(y ~ x)
resid <- residuals(model)

# Show the residuals against X, and also a histogram.
ggplot()+aes(x,resid,colour=xposition)+geom_point() 
hist(resid, breaks=30)

# Ask for the hat values and show these against X.
# The hat values are rounded to 10 decimals.
hatv <- round(hatvalues(model),10)
table(x, hatv)

Look at the plot of the residuals against X:

It is clear that for the two extreme X values the variance of the residuals is smaller then for the middle values of X. This is because the hat values at the extremes are higher, i.e. closer to 1. This is shown below in the crosstable of the hat values against X.

But first look at the histogram of the residuals:

It looks very "normal"! But actually it isn't: for each of the four X values we have a normal distribution of the residuals, each with mean 0 but with a different variance. So, the combined residuals have as distribution a so called mixture of four normals! Here we have the following mixture:

$0.01*Norm(0,s^2_1) + 0.49*Norm(0,s^2_2) + 0.49*Norm(0,s^2_3) + 0.01*Norm(0,s^2_4)$

with the $s^2$ terms being the variances for the four different X values. However, this mixture is not a normal distribution, as is explained here.

The hat values corresponding to the four X values are:

table(x,hatv)
      hatv
x      0.0001004976 0.0050756195
  -100            0          100
  -1           4900            0
  1            4900            0
  100             0          100

The hat values are rounded to 10 decimals. E.g. those of the two extreme X values show very small differences, due to imprecision of calculations needed to obtain the hat matrix, but theoretically the two hat values for X=-100 and X=100 should be equal. Given that there are only two different hat values, the mixture above can be simplified to:

$0.02*N(0,s^2_1) + 0.98*N(0,s^2_2)$

So here, the distribution of the estimated residuals $e$ is a mixture of two normals with equal mean 0, but different variances. In general, there would be a mixture of as many normals as there are different hat values, each weighted by the proportion of cases with the given X value. But only with the caveat that the true errors $\epsilon$ are normally distributed.

It is relevant to note that the way the (mixture) distribution of residuals was derived above started from the idea that for each (combination of) X value(s) there exists a distribution of residuals from which for our sample we only have drawn a "small" number of cases. The final mixture is based on the underlying distributions of residuals for the given X value(s) in the sample. So, we do not just consider the limited set of $N$ observed residuals in the sample, but tried to generalize to the underlying distributions of possible residuals that we assume to exist.