4.2 Residuals

There are many types of residuals. The model residuals are simply \(\boldsymbol{e}=\mathbf{M}_{\mathbf{X}}\boldsymbol{y}\), which can be obtained through resid for lm objects. We can verify numerically that \(\hat{{\boldsymbol{y}}} \perp \boldsymbol{e}\) and verify that \(\mathbf{X}^\top\boldsymbol{e}=\boldsymbol{0}_p\).

#Fitted values    
yhat <- fitted(ols)
#Residuals
e <- resid(ols)

#Orthogonality (by construction)
isTRUE(all.equal(c(e %*% yhat), 0))

## [1] TRUE

isTRUE(all.equal(c(e %*% X), rep(0, 2)))

## [1] TRUE

In the sequel, we will look at calculation of various variants of the residuals. The first are the standardized residuals, also internally studentized residuals. These are defined as \(r_i = e_i/\{s(1-h_{ii})^{1/2}\}\), i.e. each residual \(e_i\) is scaled by its individual variance to create homoscedastic residuals \(r_i\).

#we divide so they have the same variance - but not independent
r <- e/sqrt(s2*(1-leverage)) 
#also obtainable via rstandard(ols)
isTRUE(all.equal(rstandard(ols), r))

## [1] TRUE

Because the \(i\)th residual is used in both the numerator and in the denominator (in the calculation of \(s^2\)), the standardized (internally studentized) residual follows marginally an approximate scaled Student distribution. However, because of the use of \(s^2\) in the denominator, the entries of \(\boldsymbol{r}\) are bounded by \(\pm n-p\). They are also not independent, even if this fact is often omitted in practice. While they will be approximately centered (with mean zero and variance one), they can (and should) be recentered before undertaking visual diagnostics.

The externally studentized residuals \(t_i\) are obtained by excluding the \(i\)th observation from the calculation of the variance. The advantage of doing this is that \(\{t_i\}_{i=1}^n\) are marginally Student distributed with \(n-p-1\) degrees of freedom (but they are again not independent). These are typically the residuals that are displayed in Q-Q plots. The externally studentized residuals can be obtained with the function rstudent.

We will derive formulas for \(\hat{\boldsymbol{\beta}}_{-i}\), \(s^2_{-i}\), Cook distance and \(t_i\) later in the exercises. Two of these are used below, namely \[t_i = \frac{e_i}{[s^2_{-i}(1-h_{ii})]^{1/2}}, \qquad s^2_{-i} = \frac{(n-p)s^2 -e_i^2/(1-h_{ii})}{n-p-1}. \]

#Externally studentized residuals
smi <- influence(ols)$sigma
s2mi <-((n-2)*s2-e^2/(1-leverage))/(n-3)
isTRUE(all.equal(smi^2, s2mi))

## [1] TRUE

esr <- e/sqrt(s2mi*(1-leverage))
isTRUE(all.equal(rstudent(ols), esr))

## [1] TRUE

The last type of residual is the leave-one-out cross validation residual. These are the residuals obtained by fitting the linear model to all observations, but the \(i\)th, i.e., \(\boldsymbol{y}_{-i}= \mathbf{X}_{-i}\boldsymbol{\beta}+ \boldsymbol{\varepsilon}\). Let \(\hat{\boldsymbol{\beta}}_{-i}\) denote the OLS coefficients from this regression and \(\hat{y}_{i,-i}=\mathbf{x}_i\hat{\boldsymbol{\beta}}_{-i}\) the predicted value for the left-out \(\mathbf{x}_i\) regressor. The \(i\)th leave-one-out cross validation residual is \(e_{i,-i}=y_i - \hat{y}_{i,-i}=e_i/(1-h_{ii})\). We can use these to calculate the PRESS statistic, \(\mathsf{PRESS}=\sum_{i=1}^n e_{i, -i}^2\)

# LOOCV residuals e/(1-leverage)
loocv <- e/(1-leverage)
loocv_err <- rstandard(ols, type = "pred")
PRESS <- crossprod(loocv_err)