Regression Diagnostics

the exact \(t\) and \(F\) theory on estimation & inference is bought with assumptions: a linear mean, uncorrelated errors of constant variance, gaussian noise. diagnostics are how those assumptions are audited after the fit, and nearly all of them run through one matrix — the hat matrix \(\mathbf{H}\), whose diagonal measures how much each observation can pull the fit and whose complement \(\mathbf{I} - \mathbf{H}\) carries the residuals. 𐃏

the page builds the residual algebra first, then the influence measures it licenses (leverage, cook’s distance), then \(R^2\), then the catalogue of things that go wrong. the glm analogues (pearson and deviance residuals) live on generalised linear models and logistic regression.

the hat matrix and the residuals#

for \(\mathbf{y} = \mathbf{X}\beta + \boldsymbol\varepsilon\) with \(\mathbf{X}\) of full rank \(p\), least squares gives \(\hat{\mathbf{y}} = \mathbf{X}\hat\beta\), and every fitted value is a fixed linear combination of the responses.

symmetric plus idempotent is precisely an orthogonal projection: \(\hat{\mathbf{y}}\) is the point of the column space of \(\mathbf{X}\) nearest \(\mathbf{y}\), and \(\mathbf{X}^\top\mathbf{e} = \mathbf{0}\). two consequences are used constantly: with an intercept the residuals sum to zero, and \(\|\mathbf{y}\|^2 = \|\hat{\mathbf{y}}\|^2 + \|\mathbf{e}\|^2\) by pythagoras.

the punchline for diagnostics: raw residuals are neither independent nor equally variable, even when the errors are. the covariance is singular (rank \(N-p\)) because \(p\) linear combinations of the residuals vanish identically. the square root \(\hat\sigma\) is what R prints as “residual standard error”. 𐃏

proofs#

standardised and studentised residuals#

since \(\operatorname{Var}(e_i) = \sigma^2(1-h_{ii})\), the fair way to judge the size of a residual is to divide by its own standard error.

the name studentised is used two ways. the lecture uses it for exactly this plug-in quantity and flags \(|e_i’| > 3\) as a possible outlier; most texts reserve it for the version whose \(\hat\sigma_{(i)}\) is computed with observation \(i\) deleted, so that an outlier cannot inflate its own yardstick. in R, rstandard() is the first and rstudent() the second.

the policy matters as much as the formula: points are not deleted just because they fit badly. an outlier is a signal — a data error, or more often a feature of the problem the model has not represented — and the right response is frequently a different model in which the point is no longer extreme.

leverage and influence#

leverage#

high leverage is dangerous precisely because it is quiet: \(\operatorname{Var}(e_i) = \sigma^2(1-h_{ii})\) shrinks as \(h_{ii} \to 1\), so the fit is dragged toward a high-leverage point and its residual looks small. at \(h_{ii} = 1\) the line passes through the point whatever its response. leverage says a point could move the fit; whether it does depends on its residual.

cook’s distance#

so influence is leverage times standardised residual squared, and one fit suffices for all \(N\) deletions. the factor is the standardised residual, not the raw one (the lecture slide prints \(e_i^2\); with the raw residual the formula is off by \(\hat\sigma^2(1-h_{ii})\) and is not even dimensionless). 𐃏

cut-offs are rules of thumb: weisberg — a largest \(D_i\) well below 1 means no single deletion changes \(\hat{\mathbf{y}}\) much; fox — \(4/(N-2)\) flags noteworthy values in simple regression. better than any threshold: look for a gap in the sorted \(D_i\).

proofs#

the coefficient of determination#

the intercept-only model \(Y_i = \beta_0 + \varepsilon_i\) has \(\hat\beta_0 = \bar y\) and residual sum of squares \(\mathrm{TSS} = \sum_i(y_i - \bar y)^2\), the worst any model with an intercept can do.

\(R^2\) never decreases when a column is added — the saturated model scores 1 — so it cannot choose between models of different size; a small increase is evidence the new variable contributes little. the adjusted \(\bar R^2\), the nested \(F\)-test, aic and cross-validation are the tools for that job (model selection).

proofs#

worked: cheddar cheese#

cheese.lm <- lm(taste ~ Acetic + H2S + Lactic, data = cheddar)  # faraway package
summary(cheese.lm)

            Estimate Std. Error t value Pr(>|t|)
(Intercept) -28.8768    19.7354  -1.463  0.15540
Acetic        0.3277     4.4598   0.073  0.94198
H2S           3.9118     1.2484   3.133  0.00425 **
Lactic       19.6705     8.6291   2.280  0.03108 *
Residual standard error: 10.13 on 26 degrees of freedom
Multiple R-squared:  0.6518,   Adjusted R-squared:  0.6116
F-statistic: 16.22 on 3 and 26 DF,  p-value: 3.81e-06

what can go wrong#

each gauss–markov assumption has a characteristic signature in a residual plot and a standard remedy.

non-linearity#

on the Auto data, mpg regressed on horsepower leaves a clear U in residuals against fitted values: positive at both ends, negative in the middle. adding \(\text{horsepower}^2\) flattens it. the quadratic is still a linear model — linear in \(\beta\) — so the cure is a new column, not a new estimator; polynomial, fractional-polynomial and interaction columns are developed on analysis of variance, and flexible ones on splines and smoothing. a funnel that survives the quadratic fit is a different problem.

correlated errors#

if \(\operatorname{Var}(\boldsymbol\varepsilon)\) has off-diagonal entries, \[ \operatorname{Var}(\hat\beta) = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\operatorname{Var}(\boldsymbol\varepsilon)\mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1} \] is not the printed \(\sigma^2(\mathbf{X}^\top\mathbf{X})^{-1}\). with positive correlation (the time-series case) the printed standard errors are too small, intervals too narrow, \(p\)-values too optimistic. the size of the damage is easy to compute in the simplest case:

plot residuals in time order: at \(\rho = 0\) the signs flip almost every step; at \(\rho = 0.5\) runs appear; at \(\rho = 0.9\) the residuals form slow waves — and their spread shrinks, because the fitted line absorbs part of the drift, which is exactly how the residual variance comes to understate the error variance.

heteroscedasticity#

non-constant \(\operatorname{Var}(\varepsilon_i)\) shows as a funnel widening with the fitted values. two remedies. a concave transform of the response (\(\log y\), \(\sqrt y\)) shrinks large responses more and damps a variance that grows with the mean. weighted least squares applies when variances are known up to a constant — if \(y_i\) averages \(n_i\) raw measurements, \(\sigma_i^2 = \sigma^2/n_i\) — and minimises \(\sum_i w_i(y_i - \mathbf{x}_i^\top\beta)^2\):

\[ \hat\beta_W = (\mathbf{X}^\top\mathbf{W}\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{W}\mathbf{y}, \qquad \mathbf{W} = \operatorname{diag}(w_i),\ w_i \propto 1/\operatorname{Var}(Y_i), \]

which is the same \(\mathbf{X}^\top\mathbf{W}\mathbf{X}\) system fisher scoring solves at each step of irls: weighting by precision is what maximum likelihood does whenever variances differ (lm(..., weights = n)).

collinearity#

when predictors are nearly linearly dependent, their separate effects are poorly identified and small perturbations of the data swing the estimates. the pairwise correlation matrix catches pairs; multicollinearity among three or more predictors can hide from it entirely.

transforming the response#

proofs#

see also#

• estimation & inference — the exact \(t\)/\(F\) results these checks protect
• analysis of variance — factor models, interactions and polynomial terms, the usual fixes for a wrong mean
• generalised linear models — pearson and deviance residuals, leverage from the final irls step
• logistic regression — why residual plots degenerate for binary responses
• model selection — where \(R^2\) stops being useful and penalised criteria take over
• cross-validation — the leave-one-out identity \(\hat\beta_{-i}\) above is what makes loocv free for least squares
• regularised regression — ridge as the other answer to collinearity