What if residuals are normally distributed, but y is not?

Question

I've got a weird question. Assume that you have a small sample where the dependent variable that you're going to analyze with a simple linear model is highly left skewed. Thus you assume that $u$ is not normally distributed, because this would result in normally distributed $y$. But when you compute the QQ-Normal plot there is evidence, that the residuals are normally distributed. Thus anyone can assume that the error term is normally distributed, although $y$ is not. So what does it mean, when the error term seems to be normally distributed, but $y$ does not?

Answer 1

It is reasonable for the residuals in a regression problem to be normally distributed, even though the response variable is not. Consider a univariate regression problem where $y \sim \mathcal{N}(\beta x, \sigma^2)$. so that the regression model is appropriate, and further assume that the true value of $\beta=1$. In this case, while the residuals of the true regression model are normal, the distribution of $y$ depends on the distribution of $x$, as the conditional mean of $y$ is a function of $x$. If the dataset has a lot of values of $x$ that are close to zero and progressively fewer the higher the value of $x$, then the distribution of $y$ will be skewed to the left. If values of $x$ are distributed symmetrically, then $y$ will be distributed symmetrically, and so forth. For a regression problem, we only assume that the response is normal conditioned on the value of x.

Answer 2

@DikranMarsupial is exactly right, of course, but it occurred to me that it might be nice to illustrate his point, especially since this concern seems to come up frequently. Specifically, the residuals of a regression model should be normally distributed for the p-values to be correct. However, even if the residuals are normally distributed, that doesn't guarantee that $Y$ will be (not that it matters... ); it depends on the distribution of $X$.

Let's take a simple example (which I am making up). Let's say we're testing a drug for isolated systolic hypertension (i.e., the top blood pressure number is too high). Let's further stipulate that systolic bp is normally distributed within our patient population, with a mean of 160 & SD of 3, and that for each mg of the drug that patients take each day, systolic bp goes down by 1mmHg. In other words, the true value of $\beta_0$ is 160, and $\beta_1$ is -1, and the true data generating function is: $$ BP_{sys}=160-1*\text{daily drug dosage}+\varepsilon \\ \text{where }\varepsilon\sim\mathcal N(0, 9) $$ In our fictitious study, 300 patients are randomly assigned to take 0mg (a placebo), 20mg, or 40mg of this new medicine per day. (Notice that $X$ is not normally distributed.) Then, after an adequate period of time for the drug to take effect, our data might look like this:
(I jittered the dosages so that the points wouldn't overlap so much that they were hard to distinguish.) Now, let's check out the distributions of $Y$ (i.e., it's marginal / original distribution), and the residuals:
The qq-plots show us that $Y$ is not remotely normal, but that the residuals are reasonably normal. The kernel density plots give us a more intuitively accessible picture of the distributions. It is clear that $Y$ is tri-modal, whereas the residuals look much like a normal distribution is supposed to look.

But what about the fitted regression model, what is the effect of the non-normal $Y$ & $X$ (but normal residuals)? To answer this question, we need to specify what we might be worried about regarding the typical performance of a regression model in situations like this. The first issue is, are the betas, on average, right? (Of course, they'll bounce around some, but in the long run, are the sampling distributions of the betas centered on the true values?) This is the question of bias. Another issue is, can we trust the p-values we get? That is, when the null hypothesis true, is $p<.05$ only 5% of the time? To determine these things, we can simulate data from the above data generating process and a parallel case where the drug has no effect, a large number of times. Then we can plot the sampling distributions of $\beta_1$ and check to see if they're centered on the true value, and also check how often the relationship was 'significant' in the null case:

set.seed(123456789)                       # this make the simulation repeatable

b0 = 160;   b1 = -1;   b1_null = 0        # these are the true beta values
x  = rep(c(0, 20, 40), each=100)          # the (non-normal) drug dosages patients get

estimated.b1s  = vector(length=10000)     # these will store the simulation's results
estimated.b1ns = vector(length=10000)
null.p.values  = vector(length=10000)

for(i in 1:10000){
  residuals = rnorm(300, mean=0, sd=3)
  y.works = b0 + b1*x      + residuals
  y.null  = b0 + b1_null*x + residuals    # everything is identical except b1

  model.works = lm(y.works~x)
  model.null  = lm(y.null~x)
  estimated.b1s[i]  = coef(model.works)[2]
  estimated.b1ns[i] = coef(model.null)[2]
  null.p.values[i]  = summary(model.null)$coefficients[2,4]
}
mean(estimated.b1s)       # the sampling distributions are centered on the true values
[1] -1.000084                  
mean(estimated.b1ns)
[1] -8.43504e-05               
mean(null.p.values<.05)   # when the null is true, p<.05 5% of the time
[1] 0.0532                   

These results show that everything works out fine.

I won't go through the motions, but if $X$ had been normally distributed, with otherwise the same setup, the original / marginal distribution of $Y$ would have been normally distributed just as the residuals (albeit with a larger SD). I also didn't illustrate the effects of a skewed distribution of $X$ (which is was the impetus behind this question), but @DikranMarsupial's point is just as valid in that case, and it could be illustrated similarly.