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?

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.
