Regression Assumptions & Diagnostics
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/"
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()
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 only compares neighbors in row order, and these policies have no time order, so independence rests on the design (one row per policy) and cannot be read from this test.
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.
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 says little here, because the rows have no time order.
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, and neither test flagged a problem any more (R-squared of 0.78 on the log scale is not comparable with the 0.70 in dollars). 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.3% 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.