Therefore, if we assume that the error term is Normally distributed, doesn't that imply that the response is also Normally distributed?
Not even remotely. The way I remember this is that the residuals are normal conditional on the deterministic portion of the model. Here's a demonstration of what that looks like in practice.
I start by randomly generating some data. Then I define an outcome which is a linear function of the predictors and estimate a model.
N <- 100
x1 <- rbeta(N, shape1=2, shape2=10)
x2 <- rbeta(N, shape1=10, shape2=2)
x <- c(x1,x2)
plot(density(x, from=0, to=1))
y <- 1+10*x+rnorm(2*N, sd=1)
Let's take a look at what these residuals look like. I suspect that they should be normally distributed, since the outcome
y had iid normal noise added to it. And indeed that is the case.
plot(density(model$residuals), main="Model residuals", lwd=2)
s <- seq(-5,20, len=1000)
lines(s, dnorm(s), col="red")
plot(density(y), main="KDE of y", lwd=2)
lines(s, dnorm(s, mean=mean(y), sd=sd(y)), col="red")
Checking the distribution of y, however, we can see that it's definitely not normal! I've overlaid the density function with the same mean and variance as
y, but it's obviously a terrible fit!
The reason that this happened in this case is that the input data is not even remotely normal. Nothing about this regression model requires normality except in the residuals -- not in the independent variable, and not in the dependent variable.