Linear Models and Regularisation
Regression diagnostics
A regression can score well and still be lying — diagnostics interrogate the residuals to find missed curves, unequal spread and influential points.
- 8 min read
- 3 reading levels
- Published
Read these first
On this page 5
One lesson, three depths. Pick the one that fits you today — you can switch any time.
Beginner — No maths. Plain English.
Regression diagnostics are the health checks you run on a fitted model's mistakes, because the mistakes confess what the fit conceals.
A tailor stitches you a kurta and proudly announces "96% perfect". Before paying, you try it on. Tight at one shoulder, loose at the wrist. One sleeve is slightly short, and that repeats on every kurta he makes. The overall score hid a pattern in the errors — and patterned errors mean a fixable flaw in the method, not bad luck.
A regression's errors — called residuals, each one being "actual minus predicted" — deserve the same try-on.
Why it exists
A fitted line always produces a score, and scores can flatter. The model may have missed a curve in the data. It overpredicts in the middle and underpredicts at both ends. The errors cancel politely into a good average. Its errors may grow with the prediction — small mistakes on cheap flats, giant ones on penthouses — so the single quoted "typical error" describes neither. Or three unusual rows may be steering the whole line, as the robust regression lesson showed.
None of these show up in the headline score. All of them show up in the residuals, if you look. Diagnostics are the discipline of looking. If the model has truly captured the pattern, what remains should be patternless leftovers. Shapeless, evenly spread, belonging to no one point in particular. Any structure in the leftovers is unfinished work.
How it works
residual │ ∙ ∙ ∙ ∙ ∙ residual │ ∙ ∙ ∙
0 ─────┼─∙───∙───∙──∙─── 0 ───────┼─∙──∙──∙──∙──∙──
│∙ ∙ ∙ ∙ ∙ │ ∙ ∙ ∙ ∙
└── prediction → └── prediction →
healthy: shapeless cloud sick: fan shape —
errors grow to the rightA real example you have seen
Weather forecasts are diagnosed exactly this way. If the forecast runs 2 degrees hot every monsoon afternoon, the average error over the year still looks tiny — the bias hides in a pattern. Forecasters hunt patterned errors by season and hour; that hunting is residual diagnostics under another name.
Remember this
- The score can flatter; the residuals confess.
- A healthy model leaves a shapeless cloud — no curve, no fan, no bosses.
- Curve in the leftovers: missing feature. Fan shape: unequal spread. Lone bosses: influential rows.
What to learn next
- Multicollinearity and VIF — the design-side diagnostic, completing the checklist.
- Robust regression — what to fit when influence checks fail.
- Quantile regression — embracing unequal spread instead of testing for it.
Developer — Code and libraries.
Setup
pip install statsmodelsOutputs verified with statsmodels 0.14.6 on CPU.
Two checks on one suspicious fit
import numpy as np
import statsmodels.api as sm
from statsmodels.stats.diagnostic import het_breuschpagan
rng = np.random.default_rng(11)
size = np.sort(rng.uniform(400, 2000, 80)) # flat size in sq ft
price = 10 + 0.05 * size + rng.normal(0, 0.01 * size) # noise grows with size
X = sm.add_constant(size)
fit = sm.OLS(price, X).fit()
print(f"R^2: {fit.rsquared:.3f} coefficients: {fit.params.round(3)}")
# check 1: does the error spread grow with the prediction?
_, p_bp, _, _ = het_breuschpagan(fit.resid, X)
print(f"Breusch-Pagan p-value: {p_bp:.4f} (small = spread is not constant)")
# check 2: is there leftover curvature the line missed?
curve = np.corrcoef(fit.resid, fit.fittedvalues ** 2)[0, 1]
print(f"residual vs fitted^2 correlation: {curve:+.2f} (far from 0 = missed curve)")R^2: 0.811 coefficients: [7.048 0.054] Breusch-Pagan p-value: 0.0000 (small = spread is not constant) residual vs fitted^2 correlation: +0.03 (far from 0 = missed curve)
The walkthrough
The headline said 0.811 — respectable. The coefficients are even close to truth (10 and 0.05). A careless analyst ships it here.
Check 1 failed loudly. The Breusch-Pagan test asks whether residual size depends on the features; a p-value near zero says yes. We built the disease in — rng.normal(0, 0.01 * size) makes noise proportional to size — and the test caught it. The name for the disease: heteroscedasticity, unequal spread. Its price: predictions stay usable, but every standard error and confidence interval the model quotes is wrong — typically overconfident for the big flats.
Check 2 passed. Residuals do not track curvature (+0.03 ≈ 0): the straight-line shape is fine. Diagnostics localise problems — this fit needs a spread fix, not more polynomial terms.
The fixes for a failed spread check, in practical order: model the log of the target if effects are multiplicative; use weighted least squares (sm.WLS) with weights from the spread pattern; or keep OLS coefficients and switch to robust standard errors — fit.get_robustcov_results(cov_type="HC3") — which repair the intervals without re-modelling.
The third check — influence — takes one more line:
influence = fit.get_influence().cooks_distance[0]
print("largest Cook's distance:", influence.max().round(3))largest Cook's distance: 0.285
Cook's distance measures how much the whole fit moves if one row is deleted. A common alarm level is 1; our worst row scores 0.285 — noticeable, worth a glance, but nobody is secretly steering this model.
Common mistakes
Diagnosing only with the score. R² cannot see fan shapes, curves or bosses; that is the whole lesson. fit.summary() prints several diagnostics free — read them.
Fixing heteroscedasticity by deleting the big flats. The spread pattern is information about the world, not dirt. Model it (WLS, log target) or accommodate it (robust errors); do not amputate it.
Running diagnostics on training residuals of a flexible model. These classical checks assume the model class is rigid (a line). A boosted tree's training residuals are near zero by construction — diagnose flexible models on held-out data instead.
Treating p-values as effect sizes. With 80,000 rows instead of 80, Breusch-Pagan flags even tiny spread differences. Pair every test with a look at magnitude: plot residuals, or compare residual spread across prediction deciles.
Try it yourself
Change the noise to a constant rng.normal(0, 15, 80) and rerun — both checks should go quiet. Then instead corrupt the shape: generate price = 10 + 0.00003 * size**2 + ... and watch which check fires now. You built each disease, so you know each test told the truth.
What to learn next
- Multicollinearity and VIF — the design-side diagnostic, completing the checklist.
- Robust regression — what to fit when influence checks fail.
- Quantile regression — embracing unequal spread instead of testing for it.
Researcher — Mathematics and papers.
What OLS assumes, and which check probes what
The Gauss–Markov conditions for $y = X\beta + \varepsilon$: errors have mean zero (linearity/exogeneity), constant variance $\sigma^2$ (homoscedasticity), and no correlation; add normality for exact finite-sample inference. Violations map to diagnostics:
| Assumption | Symptom in residuals | Standard tests |
|---|---|---|
| Linearity | curvature vs fitted values | RESET (Ramsey, 1969) |
| Constant variance | fan/funnel shape | Breusch-Pagan, White |
| Independence | serial correlation | Durbin-Watson, Ljung-Box |
| Normality | heavy tails, skew in Q-Q plot | Jarque-Bera, Shapiro-Wilk |
| No dominant rows | leverage/influence spikes | hat values, Cook's distance |
The tests, briefly formalised
Breusch-Pagan (1979): regress squared residuals $\hat\varepsilon_i^2$ on the features; under homoscedasticity the explained share is negligible, and $n R^2_{aux} \sim \chi^2_{k}$ with $k$ regressors in the auxiliary fit. White (1980) generalises by including squares and cross-products, testing against any smooth variance function — and the same paper's heteroscedasticity-consistent covariance $(X^\top X)^{-1} X^\top \hat\Omega X (X^\top X)^{-1}$ founds the "robust standard errors" that let you tolerate the disease rather than cure it (HC0–HC3; MacKinnon and White, 1985 — HC3 preferred at small $n$).
Leverage and influence: the hat matrix $H = X(X^\top X)^{-1} X^\top$ gives self-sensitivity $h_{ii} = \partial \hat y_i / \partial y_i$. Cook's distance for row $i$:
$$ D_i = \frac{\hat\varepsilon_i^2}{p\,\hat\sigma^2} \cdot \frac{h_{ii}}{(1 - h_{ii})^2} $$
with $p$ the parameter count — large when a row is both surprising (big residual) and remote in feature space (big $h_{ii}$). Reference treatment: Belsley, Kuh and Welsch (1980); Cook (1977).
Durbin-Watson: $DW \approx 2(1 - \hat\rho_1)$ for lag-1 residual autocorrelation $\hat\rho_1$; 2 is healthy, values near 0 or 4 signal serial dependence — vital for time-series regressions, where autocorrelated errors fake precision.
The four-plot canon and its limits
R's plot(lm) institutionalised: residuals vs fitted (shape), Q-Q (normality), scale-location (spread), residuals vs leverage (influence). Anscombe (1973) built his famous quartet — four datasets, identical fitted lines and R², wildly different pathologies — precisely to argue that these plots are not optional.
Caveats worth teaching: single-deletion diagnostics miss groups of jointly influential points (masking again — the robust-regression connection); formal normality tests are near-meaningless at large $n$ (always rejected) and weak at small $n$; and all classical results condition on the model being chosen before seeing the data — post-selection inference is its own field (Berk et al., 2013). For ML practice, the transferable core is: inspect errors conditionally — by segment, by feature slice, by prediction magnitude — never only on average. That habit, under the name "error analysis", is exactly what modern model debugging inherited from this classical toolkit.
What to learn next
- Multicollinearity and VIF — the design-side diagnostic, completing the checklist.
- Robust regression — what to fit when influence checks fail.
- Quantile regression — embracing unequal spread instead of testing for it.