Understanding the assumptions of linear regression: residuals or data must come from a normal distribution? I'm trying to better understand linear regression and have always heard that in order to meet the assumption of normality the data has to be normally distributed. However, I've also heard it's more correct to say the residuals should be normally distributed.
Is this because if a data point falls on the linear regression line it is then safe to assume it came from a normal distribution, and in order to make conclusions about the data as a whole, the residuals (the rest of the points that don't fall on the line) must also have come from a normal distribution? But, if the data were all collected together then they must all come from the same distribution, no? Can different parts of the data follow different distributions? Why not just keep saying the data, rather then the residuals, must come from a normal distribution? What is the importance of making this distinction?
 A: tl;dr: We need normal assumption for the standard $t$-test, $F$-test, and $CI$s. We make the assumption about how the errors are distributed randomly and then do all analysis on the distribution of our estimates with this starting point. We condition on $X$ so all of the randomness in $Y|X$ comes from the errors, but since the errors are unobserved we want to be able to test the assumption with something we do observe since the residuals are observed and they are just an orthogonal projection of the errors, we can use them to test the assumption of the normality of the errors.
Full answer:
First note that the residuals are not the points that are not on the regression line, they are the distances of each point $Y_i$ from the regression line $X\hat{\beta}$. This is fundamentally different since each observation will just be a single value based on distance in the $Y$ direction, not a point in the space of $X$.
Second, note that under standard assumptions, we say that
$Y = X\beta + \epsilon$, $\epsilon$ is independent of $X$, and $\epsilon_i$ are independent of each other and identically distributed with mean $0$ and variance $\sigma^2$ (note we made no assumption on the actual distribution of $\epsilon$).
These are the assumptions made to find the OLS estimator i.e. the set of parameters $\beta$ that will minimize $\|Y - X\beta\|^2$. We also get that the OLS esitmator $\hat{\beta} = (X^TX)^{-1}X^TY$ is conditionally unbiased, given $X$ and that an unbiased estimate of the covariance matrix of $\hat{\beta}$ is $\hat{{\rm cov}}(\hat{\beta}|X) = \hat{\sigma}^2(X^TX)^{-1}$, where $\hat{\sigma}^2 = \frac{1}{n-p}\sum_i e_i^2$ and $e$ are the residuals.
The normal assumption comes in when you want to do various tests to see if some of the coefficients are significant or evaluate their confidence intervals. The $t$ and $F$ test values that you get in a standard regression output come from the assumption of the normality of the residuals, while the coefficient estimates themselves do not generally need it (and there are ways to do estimates of the CIs without it!).
Note that the assumption of normality for these purposes comes wanting to know the distribution, not just the value estimates, of our covariance and OLS estimates.
Where does the normal assumption come in? In our linear model we have that $Y = X\beta + \epsilon$. If we want to make distributional claims about our estimate, we need to know the distributions of $Y|X$ and $\epsilon|X$, since we assumed that $\epsilon$ is independent of $X$ we know that $\epsilon|X$ is just distributed as $\epsilon$ so once we have an assumption on the distribution of $\epsilon$ we get all of the information we need.
We want to know the distribution of $\hat{\beta}|X$, so we know that this is the distribution of $A_XY$—where $A_X = (X^TX)^{-1}X^T$ but since we are conditioning on $X$ this can be treated as fixed. Again noting that $X$ is being conditioned on, if we assume that the errors are distributed as a normal then we can say that $Y|X \sim X\beta + N(0, \sigma^2I)$ but since $X\beta$ are fixed we know that all the randomness is coming from $\epsilon$(!) and the $X\beta$ term is just changing the mean. So $Y|X \sim N(X\beta, \sigma^2I)$. Note that the normality of $Y$ is conditional normality, that is why we don't care about the normality of the data as a whole since marginally the distribution of $Y$ might not be normal even if the assumptions are correct.
Now we go back to
$\hat{\beta}|X = A_XY|X$ and using the well known formulas that multiplying a normal RV by a fixed matrix will multiply the mean and multiple the covariance by the matrix and the matrix transposed we get
$$\hat{\beta}|X \sim N((X^TX)^{-1}X^TX\beta, \sigma^2(X^TX)^{-1}X^TX(X^TX)^{-1}) = N(\beta, \sigma^2(X^TX)^{-1})$$
So using the normal assumption on the error and conditioning on $X$ we got the distribution of $\hat{\beta}$, but in the distribution for $\hat{\beta}$ we have a $\sigma^2$ term. This we estimate using the residuals ($\hat{\sigma}^2$ above), but in order to make distributional claims about $\hat{\sigma}^2$, we further need the normal assumption of the errors since that will give us the distribution of the residuals. Under the assumption of normality of the errors, we get that the residuals are also normally distributed (why?
\begin{align}
e &= (I - X(X^TX)^{-1}X^T)Y \\
&= (I - X(X^TX)^{-1}X^T)(X\beta + \epsilon) \\
&= X\beta - X(X^TX)^{-1}X^TX\beta + (I - X(X^TX)^{-1}X^T)\epsilon \\
&= X\beta - X\beta + (I - X(X^TX)^{-1}X^T)\epsilon \\
&= (I - X(X^TX)^{-1}X^T)\epsilon \\
&= \tilde{H}\epsilon
\end{align} which is again just an orthogonal projection times the errors) and the sum of squares residuals is $\chi^2$ distributed. We can actually observe the residuals, unlike the errors, so we can evaluate how well the model fits by using tests on the normality of the residuals since they will also be normal if the model is correct.
A: First we should distinguish between errors and residuals. If $X_1,\ldots,X_n$ are sampled from a population whose average is $\mu,$ then $X_1-\mu,\ldots,X_n-\mu$ are the errors and $X_1-\overline X,\ldots,X_n-\overline X,$ where $\overline X = (X_1+\cdots+X_n)/n$ is the mean of the sample rather than of the population, are the residuals.
A simple linear regression model says $Y_i=\beta_0+\beta_1 x_i + \varepsilon_i$ for $i=1,\ldots ,n$ and that the errors $\varepsilon_i$ are normally distributed with expected value $0$ are variance $1$ and are mutually indpendent, the parameters $\beta_0,\beta_1$ are not random (and in this context "constant" usually means not random), the parameters $\beta_0,\beta_2$ are not observable but instead must be estimated based on the observed data $(x_i,Y_i),\, i=1,\ldots,n.$
Notice that I have written capital $Y_i$ and lower-case $x_i;$ that is because the $x_i$ are treated as non-random. The word "random" in this context in effect means this: If we toss out the data $(x_i, Y_i),\,i=1,\ldots,n$ and collect new data, that which changes when that is done is "random"; that which does not change is not "random". However, realistically, in many (most) cases, the $x_i$ do change. But that are treated by the model as described above as if they were not random since what is considered is the conditional distribution of $Y_i$ given $x_i.$ (In designed experiments the $x_i$ are in some instances chosen by the experimenter and would remain the same if the experiment is replicated.)
Are the "data" normally distributed? The $Y_i$ are normally distributed because the $\varepsilon_i$ are normally distributed. The errors $\varepsilon_i$ all come from the same distribution; the $Y_i$ come from different distributions because the $x_i$ differ from each other.
But that's just a commonplace and often useful model. Is it "safe to assume" the errors are normally distributed? Certainly not. Statisticians look at the residuals to assess that.
The Gauss–Markov theorem makes weaker assumptions: the errors $\varepsilon_i$ are assumed uncorrelated rather than independent; they are assumed to have expected value $0$ and equal variances but they are not assumed to be identically distributed; hence not assumed to be normally distributed. The theorem states that among all unbiased estimators of $\beta_0,\beta_1$ that are linear as a function of $(Y_1,\ldots,Y_n),$ the one with the smallest variance is the least-squares estimator.
