Confirming the distribution of residuals in linear regression

Suppose we ran a simple linear regression $y=\beta_0+\beta_1x+u$, saved the residuals $\hat{u_i}$ and draw a histogram of distribution of residuals. If we get something which looks like a familiar distribution, can we assume that our error term has this distribution? Say, if we found out that residuals resemble normal distribution, does it make sense to assume normality of error term in population? I think it is sensible, but how can it be justified?

• Personally I find it rather difficult to assess normality from a histogram (or a kernel density plot). I would never rely on them as an "ultimate" evidence. QQ plots are much more powerful for this purpose. – user5644 Feb 9 '12 at 9:32

It all depends on how you estimate the parameters. Usually, the estimators are linear, which implies the residuals are linear functions of the data. When the errors $u_i$ have a Normal distribution, then so do the data, whence so do the residuals $\hat{u}_i$ ($i$ indexes the data cases, of course).

It's conceivable (and logically possible) that when the residuals appear to have approximately a Normal (univariate) distribution, that this arises from non-Normal distributions of errors. However, with least squares (or maximum likelihood) techniques of estimation, the linear transformation to compute the residuals is "mild" in the sense that the characteristic function of the (multivariate) distribution of the residuals cannot differ much from the cf of the errors.

In practice, we never need that the errors be exactly Normally distributed, so this is an unimportant issue. Of much greater import for the errors is that (1) their expectations should all be close to zero; (2) their correlations should be low; and (3) there should be an acceptably small number of outlying values. To check these, we apply various goodness-of-fit tests, correlation tests, and tests of outliers (respectively) to the residuals. Careful regression modeling always includes running such tests (which include various graphical visualizations of the residuals, such as supplied automatically by R's plot method when applied to an lm class).

Another way to get at this question is by simulating from the hypothesized model. Here is some (minimal, one-off) R code to do the job:

# Simulate y = b0 + b1*x + u and draw a normal probability plot of the residuals.
# (b0=1, b1=2, u ~ Normal(0,1) are hard-coded for this example.)
f<-function(n) { # n is the amount of data to simulate
x <- 1:n; y <- 1 + 2*x + rnorm(n);
model<-lm(y ~ x);
lines(qnorm(((1:n) - 1/2)/n), y=sort(model$residuals), col="gray") } # # Apply the simulation repeatedly to see what's happening in the long run. # n <- 6 # Specify the number of points to be in each simulated dataset plot(qnorm(((1:n) - 1/2)/n), seq(from=-3,to=3, length.out=n), type="n", xlab="x", ylab="Residual") # Create an empty plot out <- replicate(99, f(n)) # Overlay lots of probability plots abline(a=0, b=1, col="blue") # Draw the reference line y=x  For the case n=32, this overlaid probability plot of 99 sets of residuals shows they tend to be close to the error distribution (which is standard normal), because they uniformly cleave to the reference line$y=x$: For the case n=6, the smaller median slope in the probability plots hints that the residuals have a slightly smaller variance than the errors, but overall they tend to be normally distributed, because most of them track the reference line sufficiently well (given the small value of$n$): • things would get more interesting if you added say rexp(n) in place of rnorm(n) when generating your data. The distribution of the residuals would get WAY closer to the normal than you would think. – StasK Jul 6 '12 at 0:42 • But if we don't assume the residuals to be normal, how are the p-value of the resulting estimated coefficients calculated? What is the test statistics? – Ant Jan 15 '16 at 13:43 Let us recall the geometry of the least squares: we have the basic equation $$y_i = x_i'\beta + \epsilon_i$$ written in the matrix form as $$\mathbf{y} = \mathbf{X}\beta + \mathbf{\epsilon}$$ from which we derive the residuals $$\mathbf{e} = (I-H) \mathbf{y}$$ where $$H = X(X'X)^{-1} X'$$ is the projection matrix, or hat-matrix. We see that each individual residual$e_i$is a combination of potentially a large diagonal value$(1-h_{ii})$times its own residual$\epsilon_i$, and a bunch of small magnitude off-diagonal values$h_{ij}$times their residuals$\epsilon_j, j\neq i$. (The reason I am saying that the off-diagonal values are small is because$\sum_{j\neq i} h_{ij}^2 + h_{ii}^2 = h_{ii}$, and in fact either the diagonal or off-diagonal entries are roughly of order$O(1/n)$although this is not a very strict statement that is easily thrown off by the high leverage points.) So what happens if you sum up a lot of i.i.d. pieces with small weights? Right, you get the normal distribution by the central limit theorem. So the contribution of the off-diagonal terms to the residual will produce an essentially normal component in large samples, smoothing out the non-normality that the original distribution of the errors$\epsilon_i$may have featured. It is true of course that the major part of the residual$e_i$still comes from the own error,$(1-h_{ii})\epsilon_i$, but the interplay of all these terms may produce distributions that are much closer to the normal than the original distribution of errors. If we get something which looks like a familiar distribution, can we assume that our error term has this distribution? I would argue that you can't, since the model you have just fit is invalid if the normality assumption about the errors does not hold. (in the sense that the shape of the distribution is distinctly non-normal such as Cauchy etc.) The usual approach instead of assuming f.e. Poisson distributed errors, is to perform some form of data transformation such as log y, or 1/y in order to normalize the residuals. (also the true model might not be linear which would make the plotted residuals appear weirdly distributed even though they are in fact normal) Say, if we found out that residuals resemble normal distribution, does it make sense to assume normality of error term in population? You assumed the normality of errors once you have fit an OLS regression. Whether you have to provide arguments for that claim, depends on the type and level of your work. (it is often useful to look at what is the accepted practice in the field) Now, if the residuals do in fact appear to be normally distributed, you can pet yourself on the back, since you can use it as an empirical proof of your previous assumptions. :) Yes it is sensible. The residuals are the errors. You can also look at a normal Q-Q plot. • Yes, correct, "but how can it be justified?" What assures us that the empirical distribution of$\hat{u}_i$will approximate that of$u$? – whuber Feb 8 '12 at 19:29 • This is pedantic but the residuals are not the errors. The residuals are the observed differences from the estimated model,$y_{i} - x_{i} \hat{\beta}$. The errors are the differences from the true model$y_{i} - x_{i} \beta$. – Macro Feb 8 '12 at 19:43 • @whuber: I don't know what it is, but I'm assuming it's the same thing that justifies the sample$x$approximating the population$X\$, right? – Wayne Feb 8 '12 at 20:04
• @Wayne, I believe "it" refers to the procedure "if we found out that residuals resemble normal distribution, ... to assume normality of error terms in the population." I think you're basically right, but the subtlety is that the residuals are a product of both the sample and the method used to estimate the parameters. I find this to be a thoughtful and interesting question. – whuber Feb 8 '12 at 20:30
• @whuber I would be interested in your take on studentized versus standardized versus raw residuals. – Michelle Feb 9 '12 at 4:33