Chapter 98 · Solutions
Chapter 98 · Solutions
Five challenges, each verified in code.
Solutions to the five challenges from Chapter 98, using the statsmodels diagnostic tools.
In [1]:
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/"
In [2]:
from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.stats.stattools import jarque_bera, durbin_watson
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')
m = ols('charges ~ age + bmi + C(smoker)', ins).fit()
CHALLENGE 1
Build the diagnostic dashboard
Residual-vs-fitted, Q-Q, scale-location.
In [3]:
fig, ax = plt.subplots(1,3, figsize=(13,3.8))
ax[0].scatter(m.fittedvalues, m.resid, s=14, color=ORG, alpha=0.5); ax[0].axhline(0,color=INK,lw=1); ax[0].set_title('resid vs fitted')
sm.qqplot(m.resid, line='s', ax=ax[1], markerfacecolor=ORG, markeredgecolor=ORG, alpha=0.5); ax[1].set_title('Q-Q')
stud = m.get_influence().resid_studentized_internal
ax[2].scatter(m.fittedvalues, np.sqrt(np.abs(stud)), s=14, color=ORG, alpha=0.5); ax[2].set_title('scale-location')
plt.tight_layout(); plt.show()
CHALLENGE 2
Test every assumption
BP, JB, DW.
In [4]:
print(f'Breusch-Pagan p = {het_breuschpagan(m.resid,m.model.exog)[1]:.4f} (equal variance)')
print(f'Jarque-Bera p = {jarque_bera(m.resid)[1]:.4f} (normality)')
print(f'Durbin-Watson = {durbin_watson(m.resid):.2f} (independence)')
print('Equal variance and normality FAIL in levels; DW says little here (rows have no time order).')
Breusch-Pagan p = 0.0000 (equal variance) Jarque-Bera p = 0.0000 (normality) Durbin-Watson = 2.17 (independence) Equal variance and normality FAIL in levels; DW says little here (rows have no time order).
CHALLENGE 3
Fix it with a log
Refit on log(charges).
In [5]:
mlog = ols('np.log(charges) ~ age + bmi + C(smoker)', ins).fit()
print(f'levels BP p = {het_breuschpagan(m.resid,m.model.exog)[1]:.4f} -> log BP p = {het_breuschpagan(mlog.resid,mlog.model.exog)[1]:.3f}')
print(f'levels JB p = {jarque_bera(m.resid)[1]:.4f} -> log JB p = {jarque_bera(mlog.resid)[1]:.3f}')
print(f'R2: {m.rsquared:.3f} -> {mlog.rsquared:.3f}')
levels BP p = 0.0000 -> log BP p = 0.879 levels JB p = 0.0000 -> log JB p = 0.498 R2: 0.701 -> 0.781
CHALLENGE 4
Find influential points
Leverage, studentized residuals, Cook's distance.
In [6]:
infl = m.get_influence(); cook = infl.cooks_distance[0]
top = np.argsort(cook)[::-1][:5]
print(pd.DataFrame({'obs':top, 'cooks_D':cook[top].round(3), 'leverage':infl.hat_matrix_diag[top].round(3)}))
obs cooks_D leverage 0 174 0.261 0.028 1 34 0.134 0.020 2 245 0.091 0.029 3 269 0.083 0.029 4 236 0.072 0.039
CHALLENGE 5
Refit without the influencers
Robustness check.
In [7]:
thr = 4/len(ins); flag = np.where(cook>thr)[0]
refit = ols('charges ~ age + bmi + C(smoker)', ins[~ins.index.isin(flag)]).fit()
print(pd.DataFrame({'all':m.params.round(2), 'without':refit.params.round(2)}))
print(f'dropped {len(flag)} points; conclusions stable if coefficients barely move.')
all without Intercept -17636.37 -11146.34 C(smoker)[T.yes] 27434.75 25421.54 age 381.10 266.09 bmi 352.57 291.94 dropped 30 points; conclusions stable if coefficients barely move.
Statistics, Data Science and AI: A Visual Handbook · © 2026 John Fisher