In simple linear regression, where does the formula for the variance of the residuals come from?

Question

According to a text that I'm using, the formula for the variance of the $i^{th}$ residual is given by:

$\sigma^2\left ( 1-\frac{1}{n}-\frac{(x_{i}-\overline{x})^2}{S_{xx}} \right )$

I find this hard to believe since the $i^{th}$ residual is the difference between the $i^{th}$ observed value and the $i^{th}$ fitted value; if one were to compute the variance of the difference, at the very least I would expect some "pluses" in the resulting expression. Any help in understanding the derivation would be appreciated.

Answer 1

The intuition about the "plus" signs related to the variance (from the fact that even when we calculate the variance of a difference of independent random variables, we add their variances) is correct but fatally incomplete: if the random variables involved are not independent, then covariances are also involved -and covariances may be negative. There exists an expression that is almost like the expression in the question was thought that it "should" be by the OP (and me), and it is the variance of the prediction error, denote it $e^0 = y^0 - \hat y^0$, where $y^0 = \beta_0+\beta_1x^0+u^0$:

$$\text{Var}(e^0) = \sigma^2\cdot \left(1 + \frac 1n + \frac {(x^0-\bar x)^2}{S_{xx}}\right)$$

The critical difference between the variance of the prediction error and the variance of the estimation error (i.e. of the residual), is that the error term of the predicted observation is not correlated with the estimator, since the value $y^0$ was not used in constructing the estimator and calculating the estimates, being an out-of-sample value.

The algebra for both proceeds in exactly the same way up to a point (using $^0$ instead of $_i$), but then diverges. Specifically:

In the simple linear regression $y_i = \beta_0 + \beta_1x_i + u_i$, $\text{Var}(u_i)=\sigma^2$, the variance of the estimator $\hat \beta = (\hat \beta_0, \hat \beta_1)'$ is still

$$\text{Var}(\hat \beta) = \sigma^2 \left(\mathbf X' \mathbf X\right)^{-1}$$

We have

$$\mathbf X' \mathbf X= \left[ \begin{matrix} n & \sum x_i\\ \sum x_i & \sum x_i^2 \end{matrix}\right]$$

and so

$$\left(\mathbf X' \mathbf X\right)^{-1}= \left[ \begin{matrix} \sum x_i^2 & -\sum x_i\\ -\sum x_i & n \end{matrix}\right]\cdot \left[n\sum x_i^2-\left(\sum x_i\right)^2\right]^{-1}$$

We have

$$\left[n\sum x_i^2-\left(\sum x_i\right)^2\right] = \left[n\sum x_i^2-n^2\bar x^2\right] = n\left[\sum x_i^2-n\bar x^2\right] \\= n\sum (x_i^2-\bar x^2) \equiv nS_{xx}$$

So

$$\left(\mathbf X' \mathbf X\right)^{-1}= \left[ \begin{matrix} (1/n)\sum x_i^2 & -\bar x\\ -\bar x & 1 \end{matrix}\right]\cdot (1/S_{xx})$$

which means that

$$\text{Var}(\hat \beta_0) = \sigma^2\left(\frac 1n\sum x_i^2\right)\cdot \ (1/S_{xx}) = \frac {\sigma^2}{n}\frac{S_{xx}+n\bar x^2} {S_{xx}} = \sigma^2\left(\frac 1n + \frac{\bar x^2} {S_{xx}}\right) $$

$$\text{Var}(\hat \beta_1) = \sigma^2(1/S_{xx}) $$

$$\text{Cov}(\hat \beta_0,\hat \beta_1) = -\sigma^2(\bar x/S_{xx}) $$

The $i$-th residual is defined as

$$\hat u_i = y_i - \hat y_i = (\beta_0 - \hat \beta_0) + (\beta_1 - \hat \beta_1)x_i +u_i$$

The actual coefficients are treated as constants, the regressor is fixed (or conditional on it), and has zero covariance with the error term, but the estimators are correlated with the error term, because the estimators contain the dependent variable, and the dependent variable contains the error term. So we have

$$\text{Var}(\hat u_i) = \Big[\text{Var}(u_i)+\text{Var}(\hat \beta_0)+x_i^2\text{Var}(\hat \beta_1)+2x_i\text{Cov}(\hat \beta_0,\hat \beta_1)\Big] + 2\text{Cov}([(\beta_0 - \hat \beta_0) + (\beta_1 - \hat \beta_1)x_i],u_i) $$

$$=\Big[\sigma^2 + \sigma^2\left(\frac 1n + \frac{\bar x^2} {S_{xx}}\right) + x_i^2\sigma^2(1/S_{xx}) +2\text{Cov}([(\beta_0 - \hat \beta_0) + (\beta_1 - \hat \beta_1)x_i],u_i)$$

Pack it up a bit to obtain

$$\text{Var}(\hat u_i)=\left[\sigma^2\cdot \left(1 + \frac 1n + \frac {(x_i-\bar x)^2}{S_{xx}}\right)\right]+ 2\text{Cov}([(\beta_0 - \hat \beta_0) + (\beta_1 - \hat \beta_1)x_i],u_i)$$

The term in the big parenthesis has exactly the same structure with the variance of the prediction error, with the only change being that instead of $x_i$ we will have $x^0$ (and the variance will be that of $e^0$ and not of $\hat u_i$). The last covariance term is zero for the prediction error because $y^0$ and hence $u^0$ is not included in the estimators, but not zero for the estimation error because $y_i$ and hence $u_i$ is part of the sample and so it is included in the estimator. We have

$$2\text{Cov}([(\beta_0 - \hat \beta_0) + (\beta_1 - \hat \beta_1)x_i],u_i) = 2E\left([(\beta_0 - \hat \beta_0) + (\beta_1 - \hat \beta_1)x_i]u_i\right)$$

$$=-2E\left(\hat \beta_0u_i\right)-2x_iE\left(\hat \beta_1u_i\right) = -2E\left([\bar y -\hat \beta_1 \bar x]u_i\right)-2x_iE\left(\hat \beta_1u_i\right)$$

the last substitution from how $\hat \beta_0$ is calculated. Continuing,

$$...=-2E(\bar yu_i) -2(x_i-\bar x)E\left(\hat \beta_1u_i\right) = -2\frac {\sigma^2}{n} -2(x_i-\bar x)E\left[\frac {\sum(x_i-\bar x)(y_i-\bar y)}{S_{xx}}u_i\right]$$

$$=-2\frac {\sigma^2}{n} -2\frac {(x_i-\bar x)}{S_{xx}}\left[ \sum(x_i-\bar x)E(y_iu_i-\bar yu_i)\right]$$

$$=-2\frac {\sigma^2}{n} -2\frac {(x_i-\bar x)}{S_{xx}}\left[ -\frac {\sigma^2}{n}\sum_{j\neq i}(x_j-\bar x) + (x_i-\bar x)\sigma^2(1-\frac 1n)\right]$$

$$=-2\frac {\sigma^2}{n}-2\frac {(x_i-\bar x)}{S_{xx}}\left[ -\frac {\sigma^2}{n}\sum(x_i-\bar x) + (x_i-\bar x)\sigma^2\right]$$

$$=-2\frac {\sigma^2}{n}-2\frac {(x_i-\bar x)}{S_{xx}}\left[ 0 + (x_i-\bar x)\sigma^2\right] = -2\frac {\sigma^2}{n}-2\sigma^2\frac {(x_i-\bar x)^2}{S_{xx}}$$

Inserting this into the expression for the variance of the residual, we obtain

$$\text{Var}(\hat u_i)=\sigma^2\cdot \left(1 - \frac 1n - \frac {(x_i-\bar x)^2}{S_{xx}}\right)$$

So hats off to the text the OP is using.

(I have skipped some algebraic manipulations, no wonder OLS algebra is taught less and less these days...)

SOME INTUITION

So it appears that what works "against" us (larger variance) when predicting, works "for us" (lower variance) when estimating. This is a good starting point for one to ponder why an excellent fit may be a bad sign for the prediction abilities of the model (however counter-intuitive this may sound...).
The fact that we are estimating the expected value of the regressor, decreases the variance by $1/n$. Why? because by estimating, we "close our eyes" to some error-variability existing in the sample,since we essentially estimating an expected value. Moreover, the larger the deviation of an observation of a regressor from the regressor's sample mean, the smaller the variance of the residual associated with this observation will be... the more deviant the observation, the less deviant its residual... It is variability of the regressors that works for us, by "taking the place" of the unknown error-variability.

But that's good for estimation. For prediction, the same things turn against us: now, by not taking into account, however imperfectly, the variability in $y^0$ (since we want to predict it), our imperfect estimators obtained from the sample show their weaknesses: we estimated the sample mean, we don't know the true expected value -the variance increases. We have an $x^0$ that is far away from the sample mean as calculated from the other observations -too bad, our prediction error variance gets another boost, because the predicted $\hat y^0$ will tend to go astray... in more scientific language "optimal predictors in the sense of reduced prediction error variance, represent a shrinkage towards the mean of the variable under prediction". We do not try to replicate the dependent variable's variability -we just try to stay "close to the average".

Answer 2

I find this hard to believe since the ith residual is the difference between the ith observed value and the ith fitted value; if one were to compute the variance of the difference, at the very least I would expect some "pluses" in the resulting expression

(i) The two things are dependent (positively correlated), and while there is a + term when you expand it, there's also two - terms that come up, and some cancellation

(ii) The variance of a residual should be smaller than $\sigma^2$, since the fitted line will "pick up" any little linear component that by chance happens to occur in the errors (there's always some). There's a reduction due to the intercept and a reduction due to the slope around the center of the data whose effect is strongest at the ends of the data. You can see both terms in your question.

Algebraically:

Given $H=X(X^TX)^{-1}X^T$,

\begin{eqnarray} \text{Var}(y-\hat{y})&=&\text{Var}((I-H)y)\\ &=&(I-H)\text{Var}(y)(I-H)^T\\ &=&\sigma^2(I-H)^2\\ &=&\sigma^2(I-H) \end{eqnarray}

Hence

$$\text{Var}(y_i-\hat{y}_i)=\sigma^2(1-h_{ii})$$

In the case of simple linear regression ... this gives the answer in your question.

This answer also makes sense: since $\hat{y}_i$ is positively correlated with $y_i$, the variance of the difference should be smaller than the sum of the variances.

--

Edit: Explanation of why $(I-H)$ is idempotent.

(i) $H$ is idempotent:

$H^2=X(X^TX)^{-1}X^TX(X^TX)^{-1}X^T $ $= X\ [(X^TX)^{-1}X^TX]\ (X^TX)^{-1}X^T=X(X^TX)^{-1}X^T=H$

(ii) $(I-H)^2= I^2-IH-HI+H^2=I-2H+H=I-H$

Answer 3

Here's a hybrid of the two previous solutions. The variance of the $i$th residual, by @Glen_b's answer, is $$\operatorname{Var}(y_i-\hat y_i)=\sigma^2(1-h_{ii})$$ where $h_{ii}$ is the $(i,i)$ entry of the hat matrix $H:=X(X^TX)^{-1}X^T$. This entry can be computed as the multiplication $$h_{ii}=(X)_{i\bullet}\ (X^TX)^{-1}\ (X^T)_{\bullet i}$$ where $(X)_{i\bullet}$ denotes row $i$ of matrix $X$, and $(X^T)_{\bullet i}$ denotes column $i$ of matrix $X^T$ (i.e., row $i$ of $X$ as a column vector). To simplify @AlecosPapadopoulos's calculation of $(X^TX)^{-1}$, the trick is to observe that the residual vector remains unchanged if we switch the regression model to: $$y_i=\beta_0+\beta_1(x_i-\bar x)+u_i.$$ (In essence we've shifted the $x$ values so that they average to zero.) The design matrix for this model is $$ X=\begin{pmatrix}1&x_1-\bar x\\ 1&x_2-\bar x\\ \vdots&\vdots\\ 1&x_n-\bar x\end{pmatrix},$$ and it is straightforward to calculate: $$ X^TX=\begin{pmatrix}n&0\\ 0&S_{xx}\end{pmatrix},\qquad (X^TX)^{-1}=\begin{pmatrix}\frac1n&0\\0&\frac1{S_{xx}}\end{pmatrix}$$ and finally $$ h_{ii}= \begin{pmatrix}1 &x_i-\bar x\end{pmatrix} \begin{pmatrix}\frac1n&0\\0&\frac1{S_{xx}}\end{pmatrix} \begin{pmatrix}1 \\x_i-\bar x\end{pmatrix} =\frac1n +\frac{(x_i-\bar x)^2}{S_{xx}}. $$