One way to approach this question is to look at it in reverse: how could we begin with normally distributed residuals and arrange them to be heteroscedastic? From this point of view the answer becomes obvious: associate the smaller residuals with the smaller predicted values.
To illustrate, here is an explicit construction.
The data at the left are clearly heteroscedastic relative to the linear fit (shown in red). This is driven home by the residuals vs predicted plot at the right. But--by construction--the unordered set of residuals is close to normally distributed, as their histogram in the middle shows. (The p-value in the Shapiro-Wilk test of normality is 0.60, obtained with the
shapiro.test(residuals(fit)) issued after running the code below.)
Real data can look like this, too. The moral is that heteroscedasticity characterizes a relationship between residual size and predictions whereas normality tells us nothing about how the residuals relate to anything else.
Here is the
R code for this construction.
n <- 256
x <- (1:n)/n # The set of x values
e <- rnorm(n, sd=1) # A set of *normally distributed* values
i <- order(runif(n, max=dnorm(e))) # Put the larger ones towards the end on average
y <- 1 + 5 * x + e[rev(i)] # Generate some y values plus "error" `e`.
fit <- lm(y ~ x) # Regress `y` against `x`.
par(mfrow=c(1,3)) # Set up the plots ...
plot(x,y, main="Data", cex=0.8)
plot(predict(fit), residuals(fit), cex=0.8, main="Residuals vs. Predicted")