Predicting Medical Charges: a regularization case study¶
The final capstone, and the one about many predictors. With 15 features, some redundant, some pure noise, the star technique is regularization (ridge, lasso, elastic net), which selects the signal and tames multicollinearity automatically. Same 12-step method, on a messy insurance file, library-first with pandas, seaborn, statsmodels, and scikit-learn.
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="#059669"; DEEP="#047857"; LIGHT="#6ee7b7"; EM="#059669"; EMDEEP="#047857"; 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/"
import warnings; warnings.filterwarnings('ignore')
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression, RidgeCV, LassoCV, ElasticNetCV, lasso_path
from sklearn.model_selection import cross_val_score
pd.set_option('display.max_columns', 40)
Objective. Predict a patient's annual medical charges and identify which factors drive cost, for an insurer setting premiums or a health program targeting interventions.
The challenge. We have 15 candidate features, several of which overlap (height, weight, and BMI all measure body size) or are pure noise. Throwing them all into ordinary regression gives unstable, hard-to-read coefficients. The job is to let the data tell us which features matter, that is what regularization does.
try: raw = pd.read_csv('../../data/medical_charges.csv')
except FileNotFoundError: raw = pd.read_csv(BASE + 'medical_charges.csv')
print(raw.shape); raw.head(4)
(1018, 17)
| patient_id | age | sex | height_cm | weight_kg | bmi | children | smoker | exercise_hours_wk | diet_score | blood_pressure | cholesterol | region | income_k | noise_a | noise_b | charges | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | P5201 | 24 | Female | 174.5 | 88.9 | 29.2 | 4 | no | 1.7 | 6 | 105 | 205 | southwest | 36 | 43.6 | 37 | 8288.11 |
| 1 | P5160 | 43 | M | 173.7 | 63.8 | 21.1 | 3 | Yes | 7.0 | 9 | 120 | 195 | SOUTHEAST | 54 | 43.7 | 31 | 46593.44 |
| 2 | P5843 | 42 | M | 178.5 | 130.0 | 40.8 | 1 | yes | 0.9 | 10 | 146 | 242 | southwest | 51 | 51.4 | 91 | 84317.58 |
| 3 | P5356 | 44 | Male | 169.8 | 77.5 | 26.9 | 5 | No | 0.8 | 4 | 126 | 190 | Southwest | 109 | 39.7 | 67 | 11342.60 |
raw.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 1018 entries, 0 to 1017 Data columns (total 17 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 patient_id 1018 non-null object 1 age 1018 non-null int64 2 sex 1018 non-null object 3 height_cm 1018 non-null float64 4 weight_kg 1018 non-null float64 5 bmi 994 non-null float64 6 children 1018 non-null int64 7 smoker 1018 non-null object 8 exercise_hours_wk 1018 non-null float64 9 diet_score 1018 non-null int64 10 blood_pressure 1018 non-null int64 11 cholesterol 1018 non-null int64 12 region 1018 non-null object 13 income_k 1018 non-null int64 14 noise_a 1018 non-null float64 15 noise_b 1018 non-null int64 16 charges 1012 non-null float64 dtypes: float64(6), int64(7), object(4) memory usage: 135.3+ KB
print('smoker :', raw.smoker.value_counts().to_dict())
print('sex :', raw.sex.value_counts().to_dict())
print('region :', raw.region.unique()[:8])
print('bmi dtype:', raw.bmi.dtype, '| duplicate ids:', raw.patient_id.duplicated().sum())
print('missing :', raw.isna().sum()[lambda s: s>0].to_dict())
smoker : {'no': 575, 'No': 134, 'yes': 133, 'N': 70, 'Yes': 37, '0': 35, 'Y': 22, '1': 12}
sex : {'Female': 483, 'Male': 420, 'M': 64, 'F': 51}
region : ['southwest' 'SOUTHEAST' 'Southwest' 'southeast' 'northwest' 'northeast'
'NORTHWEST' 'SOUTHWEST']
bmi dtype: float64 | duplicate ids: 18
missing : {'bmi': 24, 'charges': 6}
The familiar mess: smoker, sex, and region in mixed spellings, 24 blank bmi values, 18 duplicate patients, and a few rows with no charges. And notice height_cm, weight_kg, and bmi are three views of the same thing, a multicollinearity trap we will confirm in step 7.
df = raw.dropna(subset=['charges']).drop_duplicates(subset='patient_id').copy()
df['charges'] = df['charges'].astype(float)
df['sex'] = df['sex'].str.strip().str.upper().map({'MALE':'Male','M':'Male','FEMALE':'Female','F':'Female'})
df['smoker'] = df['smoker'].astype(str).str.strip().str.upper().map({'YES':1,'Y':1,'1':1,'NO':0,'N':0,'0':0})
df['region'] = df['region'].astype(str).str.strip().str.lower()
df['bmi'] = pd.to_numeric(df['bmi'], errors='coerce')
df['bmi'] = df['bmi'].fillna(df['weight_kg'] / (df['height_cm']/100)**2) # BMI = weight / height^2
print(f'rows {len(raw)} -> {len(df)} (dropped 6 missing-target + 18 duplicates); smoker rate {df.smoker.mean():.0%}')
rows 1018 -> 994 (dropped 6 missing-target + 18 duplicates); smoker rate 20%
fig, ax = plt.subplots(1, 2, figsize=(12,4.2))
ax[0].hist(df.charges, bins=40, color=EM, alpha=0.8); ax[0].set(title=f'Charges are very right-skewed (skew {df.charges.skew():.1f})', xlabel='charges')
ax[0].xaxis.set_major_formatter(plt.FuncFormatter(lambda v,_: f'${v/1000:.0f}k'))
for s,c,lab in [(0,EM,'non-smoker'),(1,RED,'smoker')]:
ax[1].hist(df[df.smoker==s].charges, bins=30, color=c, alpha=0.6, label=lab)
ax[1].set(title='Smoking splits the charges', xlabel='charges'); ax[1].legend()
ax[1].xaxis.set_major_formatter(plt.FuncFormatter(lambda v,_: f'${v/1000:.0f}k'))
plt.tight_layout(); plt.show()
Charges are dramatically right-skewed, a few very expensive patients, so we will model log(charges). And one feature stands out already: smokers form a separate, much more expensive cluster. Expect it to dominate the model.
df['log_charges'] = np.log(df['charges'])
feats = ['age','bmi','height_cm','weight_kg','children','smoker','exercise_hours_wk','diet_score',
'blood_pressure','cholesterol','income_k','noise_a','noise_b']
X = pd.get_dummies(df[feats].join(df['region']), columns=['region'], drop_first=True).astype(float)
Xs = StandardScaler().fit_transform(X); y = df['log_charges'].values
print(f'log fixed the skew: {df.charges.skew():.1f} -> {df.log_charges.skew():.2f}')
print(f'{X.shape[1]} model features after encoding')
log fixed the skew: 3.7 -> 0.44 16 model features after encoding
We model log(charges) (the skew drops from about 3.7 to about 0.4), one-hot encode region, and, crucially, standardize every feature to the same scale. Standardizing matters here because a penalty that shrinks coefficients must treat a feature measured in dollars and one measured in years even-handedly.
Xc = sm.add_constant(df[['weight_kg','bmi','height_cm','age','blood_pressure','cholesterol']].astype(float))
vif = pd.Series({c: variance_inflation_factor(Xc.values, i) for i,c in enumerate(Xc.columns)}).drop('const')
print('VIF:'); print(vif.round(1).to_string())
print('\nweight/height/BMI have enormous VIF -> they are redundant (BMI is built from the other two).')
VIF: weight_kg 97.2 bmi 70.8 height_cm 23.6 age 1.1 blood_pressure 1.0 cholesterol 1.0 weight/height/BMI have enormous VIF -> they are redundant (BMI is built from the other two).
The variance inflation factors are damning: weight, BMI, and height all have VIF far above 10, because BMI is literally computed from weight and height. In ordinary regression this makes their coefficients wild and untrustworthy. Regularization will resolve it automatically, no manual dropping required.
ols = LinearRegression().fit(Xs, y)
ridge = RidgeCV(alphas=np.logspace(-3,3,50)).fit(Xs, y)
lasso = LassoCV(cv=5, random_state=0).fit(Xs, y)
enet = ElasticNetCV(l1_ratio=[.2,.5,.8,1.0], cv=5, random_state=0).fit(Xs, y)
print(f'R2 OLS={ols.score(Xs,y):.3f} Ridge={ridge.score(Xs,y):.3f} Lasso={lasso.score(Xs,y):.3f} ElasticNet={enet.score(Xs,y):.3f}')
print(f'lasso alpha={lasso.alpha_:.4f} elastic-net l1_ratio={enet.l1_ratio_}')
R2 OLS=0.866 Ridge=0.866 Lasso=0.865 ElasticNet=0.865 lasso alpha=0.0087 elastic-net l1_ratio=1.0
All four fit about equally well (R-squared ~0.87), the penalty costs almost no accuracy. The difference is in which features survive, which we read next.
coef = pd.Series(lasso.coef_, index=X.columns)
print('KEPT :', [c for c in X.columns if abs(coef[c])>1e-6])
print('ZEROED:', [c for c in X.columns if abs(coef[c])<1e-6])
cv = cross_val_score(LassoCV(cv=5, random_state=0), Xs, y, cv=5)
print(f'\n5-fold cross-validated R2 = {cv.mean():.3f} (vs in-sample {lasso.score(Xs,y):.3f}) -> not overfit')
KEPT : ['age', 'bmi', 'children', 'smoker', 'exercise_hours_wk', 'blood_pressure', 'cholesterol', 'income_k', 'region_northwest', 'region_southwest'] ZEROED: ['height_cm', 'weight_kg', 'diet_score', 'noise_a', 'noise_b', 'region_southeast'] 5-fold cross-validated R2 = 0.860 (vs in-sample 0.865) -> not overfit
# the regularization path + standardized coefficients
alphas, coefs, _ = lasso_path(Xs, y - y.mean(), n_alphas=60)
fig, ax = plt.subplots(1, 2, figsize=(13,4.6))
for i, f in enumerate(X.columns):
real = abs(coef[f])>1e-6
ax[0].plot(np.log10(alphas), coefs[i], lw=2 if real else 1, ls='-' if real else '--', color=None if real else GREY, alpha=1 if real else 0.5)
ax[0].axvline(np.log10(lasso.alpha_), color=RED, ls='--'); ax[0].set(title='Lasso path (dashed gray = zeroed)', xlabel='log10(alpha)', ylabel='coefficient')
cs = coef.reindex(coef.abs().sort_values(ascending=False).index).head(8)[::-1]
ax[1].barh(range(len(cs)), cs.values, color=[DEEP if v>0 else RED for v in cs.values])
ax[1].set_yticks(range(len(cs))); ax[1].set_yticklabels(cs.index); ax[1].axvline(0,color=INK,lw=1)
ax[1].set(title='Lasso standardized coefficients (top 8)', xlabel='effect on log-charges')
plt.tight_layout(); plt.show()
Validation verdict. The lasso drove height, weight, the two noise columns, and weak features to exactly zero, keeping BMI to represent body size and resolving the multicollinearity on its own. Cross-validated R-squared (~0.86) sits right beside the in-sample fit, so the model is not overfit. The path (left) shows every coefficient shrinking as the penalty grows; the bar (right) ranks what survived, with smoking far and away the largest driver.
# predicted vs actual, back on the dollar scale
pred = np.exp(ols.predict(Xs))
fig, ax = plt.subplots(figsize=(6,5.4))
colors = [RED if s==1 else EM for s in df.smoker]
ax.scatter(df.charges, pred, s=14, c=colors, alpha=0.4)
lims=[df.charges.min(), df.charges.max()]; ax.plot(lims, lims, ls='--', color=INK)
ax.set(title='Predicted vs actual charges (red = smoker)', xlabel='actual', ylabel='predicted')
ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda v,_: f'${v/1000:.0f}k')); ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda v,_: f'${v/1000:.0f}k'))
plt.tight_layout(); plt.show()
kept = coef[coef.abs()>1e-6].reindex(coef.abs().sort_values(ascending=False).index)
print('Cost drivers, strongest first (standardized effect on log-charges):')
for f, v in kept.head(6).items(): print(f' {f:16s} {v:+.3f}')
print(f"\nSmoking multiplies charges by roughly e^(raw smoker coef). Its standardized effect ({kept['smoker']:.2f}) dwarfs the rest.")
Cost drivers, strongest first (standardized effect on log-charges): smoker +0.571 age +0.417 blood_pressure +0.139 cholesterol +0.117 bmi +0.109 children +0.068 Smoking multiplies charges by roughly e^(raw smoker coef). Its standardized effect (0.57) dwarfs the rest.
What the model says. Smoking is the dominant cost driver by a wide margin, followed by age, then blood pressure, cholesterol, and BMI. Diet score and the noise features carry no weight; income barely registers. In plain terms: the biggest, most actionable lever on medical cost is smoking, and the body-size signal is fully captured by BMI alone, no need for height and weight too.
Actionable: smoking is the strongest predictor of cost, so it is the first place to test whether a smoking-cessation program lowers charges; blood pressure and cholesterol point to preventive-care opportunities; and the insurer can price on a short, non-redundant feature set with no loss of accuracy.
- Persist the scaler + lasso as one pipeline (
joblib), so raw patient data is standardized and scored identically. - Score new applicants for expected cost, always returning a range, and refresh premiums on a schedule.
- Target prevention: rank members by the modifiable drivers (smoking, blood pressure, cholesterol) for outreach.
- Monitor and retrain as medical costs inflate and the population changes.
- Fairness first: a cost model must be checked so it does not proxy protected attributes, and used within legal and ethical limits on what may drive premiums. Regularization aids transparency by keeping the model short and explainable.
What we did. We took an insurance file with 15 fields per patient, cleaned it (fixed duplicates and inconsistent labels, rebuilt a few missing BMI values from height and weight, standardized the charge amounts), and built a formula that estimates a patient's yearly medical cost, letting the method decide which fields actually matter.
How good is it? It explains about 87% of the variation in log cost between patients and holds up on patients it has not seen.
What drives cost, in order: smoking (by far the biggest), then age, blood pressure, cholesterol, and body-mass index. Diet score and two junk columns added nothing and were automatically dropped, income kept only a negligible weight, and height and weight were dropped too, because BMI already captures body size.
Why that dropping matters. Three of the fields (height, weight, BMI) measured the same thing, which normally makes a model unstable. The method (a 'lasso') sorted this out on its own, keeping one and discarding the rest, and it did the same for the useless fields, leaving a short, trustworthy formula.
Bottom line: smoking is the number-one cost driver by a wide margin, and a lean model built on a handful of health measures predicts cost as well as one using everything.