Why doesn't correlation of residuals matter when testing for normality?

Question

When $Y = AX + \varepsilon$ (i.e., $Y$ comes from linear regression model), $$\varepsilon \sim \mathcal{N}(0, \sigma^2 I) \hspace{1em} \Rightarrow \hspace{1em} \hat{e} = (I - H) Y \sim \mathcal{N}(0, (I - H) \sigma^2_{})$$ and in that case residuals $\hat{e}_1, \ldots, \hat{e}_n$ are correlated and not independent. But when we do regression diagnostics and want to test the assumption $\varepsilon \sim \mathcal{N}(0, \sigma^2 I)$, every textbook suggests to use Q–Q plots and statistical tests on residuals $\hat{e}$ that were designed to test whether $\hat{e} \sim \mathcal{N}(0, \sigma^2 I)$ for some $\sigma^2 \in \mathbb{R}$.

How come it doesn't matter for these tests that residuals are correlated, and not independent? It is often suggested to use standardised residuals: $$\hat{e}_i' = \frac{\hat{e}_i}{\sqrt{1 - h_{ii}}},$$ but that only makes them homoscedastic, not independent.

To rephrase the question: Residuals from OLS regression are correlated. I understand that in practice, these correlations are so small (most of the time? always?), they can be ignored when testing whether residuals came from normal distribution. My question is, why?

Answer 1

In your notation, $H$ is the projection an the column space of $X$, i.e. the subspace spanned of all regressors. Therefore $M:=I_{n}-H$ is the projection on everything orthogonal to the subspace spanned by all regressors.

If $X\in\mathbb{R}^{n\times k}$, then $\hat{e}\in\mathbb{R}^{n}$ is singular normal distributed and the elements are correlated, as you state.

The errors $\varepsilon$ are unobservable and are in general not orthogonal to the subspace spanned by $X$. For the sake of argument, assume that the error $\varepsilon\perp\operatorname{span}\left(X\right)$. If this was true, we would have $y=X\beta+\varepsilon=\tilde{y}+\varepsilon$ with $\tilde{y}\perp\varepsilon$. Since $\tilde{y}=X\beta\in\operatorname{span}\left(X\right)$, we could decompose $y$ and get the true $\varepsilon$.

Assume we have a basis $b_{1},\ldots,b_{n}$ of $\mathbb{R}^{n}$, where the first $b_{1},\ldots,b_{k}$ basis vector span the subspace $\operatorname{span}\left(X\right)$ and the remaining $b_{k+1},\ldots,b_{n}$ span $\operatorname{span}\left(X\right)^{\perp}$. In general, the error $\varepsilon=\alpha_{1}b_{1}+\ldots+\alpha_{n}b_{n}$ will have non-zero components $\alpha_{i}$ for $i\in\left\{1,\ldots,k\right\}$. This non-zero components will get mixed up with $X\beta$ and therefore can not be recovered by projection on $\operatorname{span}\left(X\right)$.

Since we can never hope to recover the true errors $\varepsilon$ and $\hat{e}$ are correlated singular $n$-dimensional normal, we could transform $\hat{e}\in\mathbb{R}^{n}\mapsto e^{*}\in\mathbb{R}^{n-k}$. There we can have that \begin{equation} e^{*}\sim\mathcal{N}_{n-k}\left(0,\sigma^{2}I_{n-k}\right) \textrm{,} \end{equation} i.e. $e^{*}$ is non-singular uncorrelated and homoscedastic normal distributed. The residuals $e^{*}$ are called Theil's BLUS residuals.

In the short paper On the Testing of Regression Disturbances for Normality you find a comparison of OLS and BLUS residuals. In the tested Monte Carlo setting the OLS residuals are superior to BLUS residuals. But this should give you some starting point.