How do residuals relate to the underlying disturbances?

Question

In the least squares method we want to estimate the unknown parameters in the model:

$$Y_j = \alpha + \beta x_j + \varepsilon_j \enspace (j=1...n)$$

Once we have done that (for some observed values), we get the fitted regression line:

$$Y_j = \hat{\alpha} + \hat{\beta}x +e_j \enspace (j =1,...n)$$

Now obviously we want to check some plots to ensure that the assumptions are fulfilled. Suppose you want to check for homoscedasticity, however, to do this we are actually checking the residuals $e_j$. Let say you examine the residual vs predicted values plot, if that shows us that heteroscedasticity is apparent, then how does that relate to the disturbance term $\varepsilon_j$? Does heteroscedasticity in the residuals imply heteroscedasticity in disturbance terms?

Answer 1

The simplest way to think about it is that your raw residuals ($e_j = y_j-\hat y_j$) are estimates of the corresponding disturbances ($\hat\varepsilon_j = e_j$). However, there are some extra complexities. For example, although we are assuming in the standard OLS model that the errors / disturbances are independent, the residuals cannot all be. In general, only $N-p-1$ residuals can be independent since you have used $p-1$ degrees of freedom in estimating the mean model and the residuals are constrained to sum to $0$. In addition, the standard deviation of the raw residuals is not actually constant. In general, the regression line is fitted such that it will be closer on average to those points with greater leverage. As a result, the standard deviation of the residuals for those points is smaller than that of low leverage points. (For more on this, it may help to read may answers here: Interpreting plot.lm(), and/or here: How to perform residual analysis for binary/dichotomous independent predictors in linear regression?)

Answer 2

The relationship between $\hat{\varepsilon}$ and $\varepsilon$ is:

$$\hat{\varepsilon} = (I-H) \varepsilon$$

where $H$, the hat matrix, is $X(X^TX)^{-1}X^T$.

Which is to say that $\hat{\varepsilon}_i$ is a linear combination of all the errors, but typically most of the weight falls on the $i$-th one.

Here's an example, using the cars data set in R. Consider the point marked in purple:

Let's call it point $i$. The residual, $\hat{\varepsilon}_i\approx 0.98\varepsilon_i +\sum_{j\neq i} w_j \varepsilon_j$, where the $w_j$ for the other errors are in the region of -0.02:

We can rewrite that as:

$\hat{\varepsilon}_i\approx 0.98\varepsilon_i +\eta_i$

or more generally

$\hat{\varepsilon}_i= (1-h_{ii})\varepsilon_i +\eta_i$

where $h_{ii}$ is the $i$-th diagonal element of $H$. Similarly, the $w_j$'s above are $h_{ij}$.

If the errors are iid $N(0,\sigma^2)$ then in this example, the weighted sum of those other errors will have a standard deviation corresponding to about 1/7th the effect of the error of the $i$th observation on its residual.

Which is to say, in well-behaved regressions, residuals can mostly be treated like a moderately noisy estimate of unobservable the error term. As we consider points further from the center, things work somewhat less nicely (the residual becomes less weighted on the error and the weights on the other errors become less even).

With many parameters, or with $X$'s not so nicely distributed, the residuals may be much less like the errors. You may like to try some examples.