Regression Assumptions and Diagnostics

A good score can hide a wrong model

A model with an R squared of 0.957 sounds finished. Run regression diagnostics on its residuals and you may find it is systematically too high at both ends of the range and too low in the middle, with errors three times larger for big predictions than small ones. The score averages all that away.

Residuals are what the model got wrong, one per row. Their pattern tells you something the score cannot: whether the errors look like noise, or like structure the model failed to capture.

Which assumptions matter depends on what you are doing. For pure prediction, you need the relationship to be captured correctly and little else. For confidence intervals, p-values or prediction intervals, you need more, which is where linear regression and its assumptions goes into the inferential side properly. This chapter takes the diagnostic angle: what to compute, what the numbers mean, and which problems are worth fixing.

Regression diagnostics without plotting anything

The standard advice is to plot residuals against fitted values. You can get the same information from a table, which also makes it checkable in a script.

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(23)
n = 500
x1 = rng.uniform(5, 60, n)
x2 = rng.normal(0, 1, n)
truth = 40 + 1.2 * x1 ** 1.6 / 10 + 3 * x2           # genuinely curved in x1
y = truth + rng.normal(0, 1, n) * (1 + 0.09 * x1)    # spread grows with x1

data = pd.DataFrame({"x1": x1, "x2": x2})
model = LinearRegression().fit(data, y)
fitted = model.predict(data)
resid = y - fitted
print(round(model.score(data, y), 4))     # 0.9566

summary = (pd.DataFrame({"fitted": fitted, "resid": resid,
                         "bucket": pd.qcut(fitted, 5, labels=False)})
           .groupby("bucket")
           .agg(fitted_mid=("fitted", "mean"),
                mean_resid=("resid", "mean"),
                sd_resid=("resid", "std")))
print(summary.round(2).to_string())
#         fitted_mid  mean_resid  sd_resid
# bucket
# 0            41.79        3.21      3.33
# 1            60.49       -2.50      3.18
# 2            77.20       -3.09      4.10
# 3            92.65       -0.78      5.03
# 4           108.78        3.16      5.70

Two diagnoses in one table.

The mean_resid column should be near zero everywhere. Instead it runs positive, negative, negative, negative, positive: a U-shape. The model underpredicts at both ends and overpredicts in the middle, which is what fitting a straight line to a curve looks like. This is the one problem you must fix regardless of purpose, because the model is wrong in a patterned way across the whole range.

The sd_resid column should be roughly constant. It climbs from 3.33 to 5.70, so errors are about 70% larger for the biggest predictions. That is heteroscedasticity, and it does not bias the predictions; it means any single prediction interval you quote is too wide at one end and too narrow at the other.

Fix the curvature first, since it often removes the apparent spread problem too.

from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline

quad = make_pipeline(PolynomialFeatures(2, include_bias=False),
                     LinearRegression()).fit(data, y)
resid2 = y - quad.predict(data)
print(pd.DataFrame({"r": resid2, "b": pd.qcut(quad.predict(data), 5, labels=False)})
      .groupby("b")["r"].mean().round(3).to_dict())
# {0: 0.224, 1: -0.686, 2: -0.153, 3: 0.561, 4: 0.054}

The U-shape is gone. Adding a squared term, as covered in polynomial regression, was the whole fix.

Normality, and when it actually matters

from scipy import stats

print(round(float(stats.skew(resid)), 3))        # 0.383
print(round(float(stats.kurtosis(resid)), 3))    # 2.194
print(f"{stats.shapiro(resid).pvalue:.2e}")      # 7.99e-07

Shapiro-Wilk rejects normality decisively, and the model still predicts well. Both things are true at once, which confuses people.

Normality of residuals is needed for exact p-values and confidence intervals on coefficients. It is not needed for the predictions themselves. With 500 rows the central limit theorem covers most inference anyway, and the test becomes so sensitive at large sample sizes that it rejects almost everything.

My position: check skew and kurtosis to spot a badly misspecified model, and ignore the formal normality test unless you are reporting inferential statistics. Skew of 0.38 and kurtosis of 2.19 here indicate heavier tails than a normal distribution, which is worth noticing and not worth acting on.

Autocorrelation and influential rows

Independence of errors is the assumption most often violated silently, because it only shows up when rows have an order. The Durbin-Watson statistic measures whether each residual resembles the one before it, running from 0 to 4 with 2 meaning no correlation.

dw = float((np.diff(resid) ** 2).sum() / (resid ** 2).sum())
print(round(dw, 3))     # 1.912

series = np.zeros(n)
e = rng.normal(0, 1, n)
for i in range(1, n):
    series[i] = 0.85 * series[i - 1] + e[i]
print(round(float((np.diff(series) ** 2).sum() / (series ** 2).sum()), 3))   # 0.345

1.912 is fine. A strongly autocorrelated series scores 0.345, and anything below about 1.5 means your rows are not independent. On time-ordered data this is the norm rather than the exception, and it means your error estimates are optimistic.

Leverage and influence identify single rows that are steering the fit.

Xd = np.hstack([np.ones((n, 1)), data.to_numpy()])
H = Xd @ np.linalg.pinv(Xd.T @ Xd) @ Xd.T
leverage = np.diag(H)
p = Xd.shape[1]
cooks = (resid ** 2 / (p * (resid ** 2).sum() / (n - p))) * (leverage / (1 - leverage) ** 2)

print(round(float(leverage.mean()), 4), round(p / n, 4))        # 0.006 0.006
print(int((leverage > 2 * p / n).sum()))                        # 25
print(round(float(cooks.max()), 4), int((cooks > 4 / n).sum())) # 0.1448 35

Average leverage always equals the number of coefficients divided by the row count. Rows above twice that have unusual input values; rows with Cook’s distance above 4/n are changing the fitted coefficients noticeably. The maximum here is 0.145, comfortably below the threshold of 1 that marks a genuinely dominant row.

Fix curvature with transformations or a more flexible model. Fix autocorrelation by splitting on time rather than at random. Investigate influential rows individually rather than deleting them. Heteroscedasticity and non-normality are usually worth recording and leaving alone, unless you are quoting intervals. Reading what the fitted model says about each feature is covered in reading a model’s output, and the normality test is documented under scipy.stats.shapiro.

Frequently Asked Questions

What do residual plots tell you about a regression model?

Whether the errors look like noise or structure. A pattern in the average residual across the range of predictions means the model is missing a relationship, usually curvature. Changing spread means heteroscedasticity, which affects prediction intervals but not the predictions themselves.

Do residuals need to be normally distributed?

Not for prediction. Normality matters for exact p-values and confidence intervals on coefficients, and with a few hundred rows the central limit theorem covers most of that anyway. The model above failed a Shapiro-Wilk test decisively while achieving an R squared of 0.957.

What is a good Durbin-Watson value?

Around 2 means no autocorrelation between consecutive residuals. Below about 1.5 indicates positive autocorrelation, which is common in time-ordered data and means your error estimates are too optimistic. A strongly autocorrelated series scored 0.345 in the example above.

Key Takeaways

• Group residuals into buckets by fitted value and check both the mean and the standard deviation, since one table diagnoses curvature and changing spread at once.
• Fix any pattern in mean residuals before anything else, because that is the model being systematically wrong rather than merely imprecise.
• Treat heteroscedasticity and non-normal residuals as notes rather than blockers unless you are quoting prediction intervals or p-values.
• Compute Durbin-Watson on any time-ordered data and split by time when it falls below 1.5, as correlated errors make held-out scores optimistic.
• Investigate rows with Cook’s distance above 4/n individually, rather than deleting them because they are inconvenient.