Help clarify the implication of normality in an Ordinary Least Square (OLS) Regression

Question

If the dependent variable is normally distributed for a fixed set of predictor values, then the residual values should be normally distributed with a mean of 0.

I have two questions based on the above paragraph:

First, I don't understand why the above paragraph must be true. Could someone please provide more details why the above must be true? Would be great if the answer clearly states what properties of random variables are being used to arrive at the above conclusion and any other critical assumptions not explicitly stated in the paragraph. For simplicity, please assume a simple linear regression.

The second question is I am not really sure what the phrase "...for a fixed set of predictor values.." mean. Does this mean that:

Suppose we have the following pairs of observations ($x$ is the predictor value, $y$ is the dependent value):

$$( x_1 , y_1) , (x_2,y_2), (x_3,y_3)$$

Does this mean that for each $x_i$, there will be corresponding predicted value $y_i$ and all of these $y_i$ will follow a normal distribution with possibly different means for $y_i$ but same variance, where $i \in [1,3]$?

Answer 1

Assuming you are working with some model like $$ y_i = x_i'\beta + e_i$$

If you assume: $$ y_i | x_i \sim N(\mu, \sigma^2) \qquad (1)$$

This implies: $$ e_i | x_i = (y_i - x_i'\beta) | x_i$$

Since you are conditioning on $x_i$, the $x_i'\beta$ term is not random. As long as you include an intercept, $y_i - x_i'\beta$ will have mean 0. Therefore,

$$e_i | x_i \sim N(0, \sigma^2)$$

You should be aware that (1) is a very strong assumption, and the model can be estimated under different (weaker) assumptions.

The bit about fixed regressors is an assumption that I think is more common in some fields (perhaps, lab sciences) than in others (social sciences where $x$ includes covariates that are part of a sample that are random along with $y$). This assumption gets you straight-away that the distribution of $y_i|x_i$ does not depend on $x_i$. I.e. $f(y_i|x_i) = f(y_i)$ for any distribution that $y$ follows.

Answer 2

The original claim probably implied or mentioned elsewhere that the conditional distributions have constant variance. However, without that assumption, the claim is false, as the below simulation shows.

set.seed(2024)
N <- 10000
x <- rf(N, 2, 3) 
Ey <- x + 1
e <- rnorm(N, 0, Ey)
y <- Ey + e
L <- lm(y ~ x)
r <- resid(L)
qqnorm(r)
qqline(r)

We see the residuals are non-normal, despite the fact that the error term e is composed of independent Gaussians with zero mean (but different variances).

Visualizing the empirical distribution function of the residuals vs the theoretical Gaussian distribution for the empirical mean and variance may be helpful.

library(ggplot2) # You don't have to use ggplot2 to calculate the distributions
d0 <- data.frame(
  Residuals = r,
  CDF = ecdf(r)(r),
  Distribution = "Empirical"
)
d1 <- data.frame(
  Residuals = r,
  CDF = pnorm(r, mean(r), sd(r)),
  Distribution = "Theoretical Gaussian"
)
d <- rbind(d0, d1)
ggplot(d, aes(x = Residuals, y = CDF, col = Distribution)) +
  geom_line() +
  xlim(quantile(r, c(0.01, 0.99))) +
  theme(legend.position = "bottom")

Those do not look like the same distribution, and with ten-thousand residuals, they should be pretty close if the residuals are Gaussian.

This shows the importance of requiring the conditional distributions to have equal variances, which is not the case in this simulation.

Compare this to what happens if the variance is required to be the same for each conditional distribution.

set.seed(2024)
N <- 10000
x <- rf(N, 2, 3)
Ey <- x + 1
e <- rnorm(N, 0, 1)#Ey)
y <- Ey + e
L <- lm(y ~ x)
r <- resid(L)
qqnorm(r)
qqline(r)

library(ggplot2)
d0 <- data.frame(
  Residuals = r,
  CDF = ecdf(r)(r),
  Distribution = "Empirical"
)
d1 <- data.frame(
  Residuals = r,
  CDF = pnorm(r, mean(r), sd(r)),
  Distribution = "Theoretical Gaussian"
)
d <- rbind(d0, d1)
ggplot(d, aes(x = Residuals, y = CDF, col = Distribution)) +
  geom_line() +
  xlim(quantile(r, c(0.01, 0.99))) +
  theme(legend.position = "bottom")