←Back to the bookChapter 95 · Regression Assumptions & Diagnostics

Chapter 95 · Regression Assumptions & Diagnostics

Regression Assumptions & Diagnostics

Build the residual dashboard, run Breusch-Pagan / Jarque-Bera / Durbin-Watson, measure influence with Cook's distance, and fix a real model with a log transform, all library-first with statsmodels.

Statistics, Data Science and AI: A Visual Handbook · John Fisher · 2026

Regression diagnostics: checking (and fixing) the LINE assumptions

Library-first: statsmodels supplies the residual dashboard, the formal tests, and the influence measures. We diagnose a naive insurance model, expose two violations, cure them with one log transform, and hunt for influential points.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy import stats
import seaborn as sns   # seaborn = high-level statistical plots (heatmaps, regplots, pairplots)
import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.nonparametric.smoothers_lowess import lowess
from sklearn.linear_model import LinearRegression, Ridge, Lasso, LogisticRegression
from sklearn.preprocessing import StandardScaler, PolynomialFeatures
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.metrics import r2_score, mean_squared_error
ORG="#ea580c"; DEEP="#c2410c"; LIGHT="#fdba74"; INK="#1a2138"; GRID="#e6e9f2"; GREEN="#059669"; RED="#ef4444"; AMBER="#d97706"; BLUE="#2563eb"; PUR="#9333ea"; GREY="#94a3b8"; SLATE="#475569"
plt.rcParams.update({"figure.facecolor":"white","axes.facecolor":"white","figure.dpi":110,"font.size":11,
   "axes.edgecolor":GRID,"axes.grid":True,"grid.color":GRID,"axes.axisbelow":True,"axes.spines.top":False,
   "axes.spines.right":False,"axes.titlesize":12,"axes.titleweight":"bold","legend.frameon":False})
sns.set_style("whitegrid")
BASE="https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"

DEMO 1 · FIT & DASHBOARD

The four-plot diagnostic dashboard

Fit charges on age, BMI, and smoking, then build the classic residual dashboard.

try:
    ins = pd.read_excel('../../data/regression-assumptions-and-diagnostics--insurance.xlsx', sheet_name='Policies')
except FileNotFoundError:
    ins = pd.read_excel(BASE + 'regression-assumptions-and-diagnostics--insurance.xlsx', sheet_name='Policies')
print(ins.shape); ins.head(3)

(300, 5)

	policy_id	age	bmi	smoker	charges
0	P2000	39	29.0	no	5445.32
1	P2001	43	17.1	no	4725.71
2	P2002	20	36.5	no	3894.55

m = ols('charges ~ age + bmi + C(smoker)', ins).fit()
resid, fitted = m.resid, m.fittedvalues
infl = m.get_influence(); stud = infl.resid_studentized_internal
fig, ax = plt.subplots(2, 2, figsize=(11,8))
ax[0,0].scatter(fitted, resid, s=16, color=ORG, alpha=0.5); ax[0,0].axhline(0,color=INK,lw=1)
lo = lowess(resid, fitted, frac=0.6); ax[0,0].plot(lo[:,0], lo[:,1], color=DEEP, lw=2)
ax[0,0].set(title='Residuals vs fitted (funnel?)', xlabel='fitted', ylabel='residual')
sm.qqplot(resid, line='s', ax=ax[0,1], markerfacecolor=ORG, markeredgecolor=ORG, alpha=0.5)
ax[0,1].set_title('Q-Q plot (normal?)')
ax[1,0].scatter(fitted, np.sqrt(np.abs(stud)), s=16, color=ORG, alpha=0.5)
ax[1,0].set(title='Scale-location (spread trend?)', xlabel='fitted', ylabel='sqrt(|std resid|)')
ax[1,1].scatter(infl.hat_matrix_diag, stud, s=16, color=ORG, alpha=0.5); ax[1,1].axhline(0,color=INK,lw=1)
ax[1,1].set(title='Residuals vs leverage', xlabel='leverage', ylabel='studentized residual')
plt.tight_layout(); plt.show()

DEMO 2 · STEP 4 TESTS

Turn the plots into pass/fail decisions

Breusch-Pagan, Jarque-Bera, and Durbin-Watson put numbers on what the dashboard hints.

from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.stats.stattools import jarque_bera, durbin_watson
def report(model, name):
    bp = het_breuschpagan(model.resid, model.model.exog)[1]
    jb = jarque_bera(model.resid)[1]
    dw = durbin_watson(model.resid)
    print(f'{name:14s} R2={model.rsquared:.3f} | BP p={bp:.4f} {"OK" if bp>.05 else "FAIL"} | JB p={jb:.4f} {"OK" if jb>.05 else "FAIL"} | DW={dw:.2f}')
report(m, 'charges')

charges        R2=0.701 | BP p=0.0000 FAIL | JB p=0.0000 FAIL | DW=2.17

Reading the three tests. Each turns a residual plot into a yes/no decision:

• Breusch-Pagan tests equal variance. A p-value below 0.05 means the spread of the residuals changes with the fitted value (the funnel), so the assumption FAILS.
• Jarque-Bera tests normality of the residuals. Below 0.05 means they are too skewed or heavy-tailed to be Normal.
• Durbin-Watson tests independence. It runs 0 to 4; near 2 is good, far from 2 signals residuals that track each other (autocorrelation).

For the raw-dollar model both Breusch-Pagan and Jarque-Bera fail, confirming the funnel and skew from the dashboard. Durbin-Watson is fine because the rows are independent policies.

DEMO 3 · THE FIX

One log transform repairs both violations

Model log(charges) and re-run the same checks, the funnel flattens and the Q-Q straightens.

mlog = ols('np.log(charges) ~ age + bmi + C(smoker)', ins).fit()
report(m,    'charges')
report(mlog, 'log(charges)')
fig, ax = plt.subplots(1, 2, figsize=(11,4.3))
ax[0].scatter(m.fittedvalues, m.resid, s=14, color=RED, alpha=0.5); ax[0].axhline(0,color=INK,lw=1)
ax[0].set(title='BEFORE: charges (funnel)', xlabel='fitted', ylabel='residual')
ax[1].scatter(mlog.fittedvalues, mlog.resid, s=14, color=GREEN, alpha=0.5); ax[1].axhline(0,color=INK,lw=1)
ax[1].set(title='AFTER: log(charges) (flat band)', xlabel='fitted', ylabel='residual')
plt.tight_layout(); plt.show()

charges        R2=0.701 | BP p=0.0000 FAIL | JB p=0.0000 FAIL | DW=2.17
log(charges)   R2=0.781 | BP p=0.8794 OK | JB p=0.4984 OK | DW=1.94

# Interpret the log-scale coefficients as percentage / multiplicative effects
print('log(charges) coefficients:')
for k,v in mlog.params.items():
    print(f'  {k:20s} {v:+.4f}   -> x{np.exp(v):.2f} per unit')
print(f"\nSmokers pay about {np.exp(mlog.params['C(smoker)[T.yes]']):.1f}x more, all else equal.")

log(charges) coefficients:
  Intercept            +6.6207   -> x750.48 per unit
  C(smoker)[T.yes]     +1.4650   -> x4.33 per unit
  age                  +0.0324   -> x1.03 per unit
  bmi                  +0.0308   -> x1.03 per unit

Smokers pay about 4.3x more, all else equal.

DEMO 4 · INFLUENCE

Leverage, Cook's distance, and a robustness check

Find the points that most move the fit, then confirm the conclusions survive without them.

infl = mlog.get_influence(); cook = infl.cooks_distance[0]
fig, ax = plt.subplots(figsize=(8,3.6))
ax.stem(np.arange(len(cook)), cook, linefmt=ORG, markerfmt=' ', basefmt=' ')
thr = 4/len(ins)
ax.axhline(thr, color=RED, ls='--', lw=1, label=f'4/n rule = {thr:.3f}')
ax.set(title="Cook's distance by observation", xlabel='observation', ylabel="Cook's D"); ax.legend()
plt.tight_layout(); plt.show()
flag = np.where(cook > thr)[0]
print(f'{len(flag)} points exceed the 4/n rule; max Cook D = {cook.max():.3f}')

15 points exceed the 4/n rule; max Cook D = 0.032

# Robustness: drop the flagged points and refit; do the coefficients hold?
keep = ~ins.index.isin(flag)
refit = ols('np.log(charges) ~ age + bmi + C(smoker)', ins[keep]).fit()
comp = pd.DataFrame({'all': mlog.params.round(4), 'without_influential': refit.params.round(4)})
print(comp)
print('\nCoefficients barely move: the conclusions do not hinge on a few points.')

                     all  without_influential
Intercept         6.6207               6.6073
C(smoker)[T.yes]  1.4650               1.4694
age               0.0324               0.0337
bmi               0.0308               0.0296

Coefficients barely move: the conclusions do not hinge on a few points.

Putting it all together: the diagnostic workflow

This chapter is steps 1 and 4 of the workflow up close. Here is the full run and how each tool was used.

Step 1, explore. The charges distribution is heavily right-skewed, a long tail of expensive policies. A skewed, positive response is a red flag for the equal-variance and normality assumptions before we even fit.

Step 2, fit. We fit charges ~ age + bmi + C(smoker) in raw dollars, which posted a respectable R-squared of 0.70.

Steps 3-4, evaluate and check the conditions. The four-plot dashboard (residuals-vs-fitted, Q-Q, scale-location, residuals-vs-leverage) is the analyst's first look. We then made the decision formal: Breusch-Pagan (equal variance) and Jarque-Bera (normality) both failed for the raw model, while Durbin-Watson (independence) was fine.

The remedy. Because the response was skewed and positive, we modeled log(charges). Re-running the same checks, the funnel flattened, the Q-Q straightened, both tests now passed, and R-squared even rose from 0.70 to 0.78. That is the feedback loop: diagnose, apply the matching fix, re-check.

Influence. Leverage measures an unusual x, a large residual an unusual y; Cook's distance combines them to flag points that actually move the fit. We flagged the highest-Cook's-distance policies with the 4/n rule, dropped them, and re-fit, the coefficients barely moved, so the conclusions do not hinge on a few rows. (Investigate influential points; do not blindly delete them.)

Interpretation. On the log scale a coefficient b means the response multiplies by e^b per unit: age adds about 3.2% per year, BMI about 3.1% per point, and being a smoker multiplies charges by about 4.3x, the biggest driver. To predict a dollar amount, exponentiate the log-scale prediction.

Takeaway. A high R-squared does not certify a model. The residual diagnostics are what let you trust the standard errors, p-values, and intervals; when they fail, a transform or a robust method is the cure.

---

Statistics, Data Science and AI: A Visual Handbook · © 2026 John Fisher

←Back to the bookChapter 95 · Regression Assumptions & Diagnostics