Regression, Regularisation and Leakage

Week 4 of 16 · Supervised learning · 7 days

Full curriculum

Week 04 · Supervised learning

Regression, Regularisation and Leakage

Week 04 · Day 1 of 7

Linear Regression on a Real Target

Fitting, reading the score against the right baseline, and what residuals reveal

By Pawan Rama Mali1187 words

Week 3 fitted a line by hand. This week fits real regressions, and the first job is choosing a target worth predicting.

The problem

Predict monthly_charges from everything else. It is a genuine business question, what should this customer be paying, given who they are, and unlike week 3's toy version it has real signal in it.

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')
df = df.dropna(subset=['monthly_charges'])

NUM = ['tenure_months', 'support_calls']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
TARGET = 'monthly_charges'
print('rows %d' % len(df))
print('target: mean %.2f  std %.2f  min %.2f  max %.2f'
      % (df[TARGET].mean(), df[TARGET].std(), df[TARGET].min(), df[TARGET].max()))
print('\nby internet service:')
print(df.groupby('internet_service')[TARGET].agg(['size', 'mean', 'std']).round(2))

rows 2835
target: mean 59.41  std 23.47  min 15.00  max 105.96

by internet service:
                  size   mean   std
internet_service                   
DSL                965  55.19  8.99
Fibre optic       1282  79.56  9.18
No internet        588  22.38  7.34

Three tight clusters. Internet service almost determines the price, which is a promising sign: there is real structure for a model to find.

Fit it

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')
df = df.dropna(subset=['monthly_charges'])

NUM = ['tenure_months', 'support_calls']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
TARGET = 'monthly_charges'
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, OneHotEncoder

def preprocessor(num=NUM, cat=CAT):
    return ColumnTransformer([
        ('num', Pipeline([('i', SimpleImputer(strategy='median')),
                          ('s', StandardScaler())]), num),
        ('cat', Pipeline([('i', SimpleImputer(strategy='most_frequent')),
                          ('o', OneHotEncoder(handle_unknown='ignore'))]), cat),
    ])
from sklearn.model_selection import train_test_split

X = df[NUM + CAT]
y = df[TARGET]
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.25, random_state=42)
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_error, root_mean_squared_error, r2_score

model = Pipeline([('prep', preprocessor()),
                  ('reg', LinearRegression())]).fit(X_tr, y_tr)
pred = model.predict(X_te)

print('MAE  %6.3f' % mean_absolute_error(y_te, pred))
print('RMSE %6.3f' % root_mean_squared_error(y_te, pred))
print('R2   %6.4f' % r2_score(y_te, pred))
print('\nstd of the target was %.2f, so predicting the mean'
      ' would give RMSE about that.' % y_te.std())

MAE   7.200
RMSE  9.033
R2   0.8498

std of the target was 23.32, so predicting the mean would give RMSE about that.

Compare RMSE against the target's standard deviation

Predicting the mean for everybody gives an RMSE equal to the standard deviation, about 23.3 here. Our model gets just over 9. That ratio is the honest way to read a regression score, and it is what R² formalises: 0.85 means the model removed 85 percent of the variance the mean-only model left behind.

Residuals are where the truth is

Residual: Actual minus predicted, for one row. A well-fitted linear model leaves residuals that are centred on zero, roughly symmetric, and show no pattern against the prediction. Any pattern that remains is signal the model failed to use.

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')
df = df.dropna(subset=['monthly_charges'])

NUM = ['tenure_months', 'support_calls']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
TARGET = 'monthly_charges'
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, OneHotEncoder

def preprocessor(num=NUM, cat=CAT):
    return ColumnTransformer([
        ('num', Pipeline([('i', SimpleImputer(strategy='median')),
                          ('s', StandardScaler())]), num),
        ('cat', Pipeline([('i', SimpleImputer(strategy='most_frequent')),
                          ('o', OneHotEncoder(handle_unknown='ignore'))]), cat),
    ])
from sklearn.model_selection import train_test_split

X = df[NUM + CAT]
y = df[TARGET]
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.25, random_state=42)
from sklearn.linear_model import LinearRegression

model = Pipeline([('prep', preprocessor()),
                  ('reg', LinearRegression())]).fit(X_tr, y_tr)
resid = y_te - model.predict(X_te)

print('mean   %8.4f   <- should be near zero' % resid.mean())
print('median %8.4f' % resid.median())
print('std    %8.4f' % resid.std())
print('skew   %8.4f   <- should be near zero' % resid.skew())
print('\nlargest under-predictions:')
print(resid.nlargest(3).round(2).to_string())
print('largest over-predictions:')
print(resid.nsmallest(3).round(2).to_string())

mean    -0.4344   <- should be near zero
median  -0.6659
std      9.0292
skew     0.0188   <- should be near zero

largest under-predictions:
132     26.92
1074    25.72
2735    25.56
largest over-predictions:
3076   -30.66
767    -27.13
2637   -26.33

The four assumptions, and which ones matter

Assumption	What breaks if violated	How much you should care
Linearity	Systematic curved pattern in residuals	A lot. The model is simply wrong
Independent errors	Standard errors understated	A lot for time series, rarely otherwise
Constant variance	Predictions unreliable at the extremes	Moderate, affects intervals more than point estimates
Normal errors	Confidence intervals are off	Least of the four; irrelevant for pure prediction

These are assumptions about residuals, not about your columns

A common misreading: people check whether their features are normally distributed. Linear regression makes no such requirement. The normality assumption is about the errors, and it only matters when you want confidence intervals on the coefficients. For prediction you can ignore it.

Check linearity by plotting residuals against predictions

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')
df = df.dropna(subset=['monthly_charges'])

NUM = ['tenure_months', 'support_calls']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
TARGET = 'monthly_charges'
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, OneHotEncoder

def preprocessor(num=NUM, cat=CAT):
    return ColumnTransformer([
        ('num', Pipeline([('i', SimpleImputer(strategy='median')),
                          ('s', StandardScaler())]), num),
        ('cat', Pipeline([('i', SimpleImputer(strategy='most_frequent')),
                          ('o', OneHotEncoder(handle_unknown='ignore'))]), cat),
    ])
from sklearn.model_selection import train_test_split

X = df[NUM + CAT]
y = df[TARGET]
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.25, random_state=42)
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

model = Pipeline([('prep', preprocessor()),
                  ('reg', LinearRegression())]).fit(X_tr, y_tr)
pred = model.predict(X_te)
resid = y_te - pred

fig, ax = plt.subplots(figsize=(7, 4))
ax.scatter(pred, resid, s=8, alpha=0.4, color='#2563eb')
ax.axhline(0, color='#dc2626', linewidth=1)
ax.set_xlabel('Predicted')
ax.set_ylabel('Residual')
fig.tight_layout()
fig.savefig('residuals.png', dpi=120)
plt.close(fig)

# The same check, numerically: is the mean residual flat across the range?
bins = pd.qcut(pred, 5, labels=False)
print(pd.DataFrame({'pred_bin': bins, 'resid': resid.to_numpy()})
      .groupby('pred_bin')['resid'].agg(['size', 'mean', 'std']).round(3))

          size   mean    std
pred_bin                    
        142  0.577  7.503
        142  0.261  9.132
        141 -0.570  9.695
        142 -0.045  9.648
        142 -2.397  8.816

The bin means sit between +0.58 and −2.40 against a residual standard deviation of 9, so there is no strong trend, a slight tendency to over-predict the most expensive customers, and little else. That is roughly what a correctly specified linear model looks like. A clear U-shape here would mean the relationship is curved and you need the polynomial features from day 3.

Day 1 takeaway

Judge a regression by comparing its RMSE against the standard deviation of the target, which is what a mean-only model would achieve. Then look at residuals: centred on zero, no pattern against the prediction. The normality assumption applies to errors, not features, and matters least of the four.

Week 04 · Day 2 of 7

Regression Metrics and Their Traps

MAE, RMSE, R-squared and MAPE, what each hides

By Pawan Rama Mali754 words

Four metrics, each answering a different question, and each with a way of misleading you.

What each one actually measures

Metric	Units	Reads as	Weakness
MAE	Same as target	Typical error size	Treats a huge miss as merely several small ones
RMSE	Same as target	Error, weighted toward big misses	One bad outlier dominates it
R²	None	Share of variance explained	Depends on the spread of your test set
MAPE	Percent	Relative error	Explodes when the target is near zero

import numpy as np
from sklearn.metrics import (mean_absolute_error, root_mean_squared_error,
                             r2_score)

truth = np.array([20.0, 40.0, 60.0, 80.0, 100.0])

spread = np.array([25.0, 35.0, 65.0, 75.0, 105.0])   # five misses of 5
single = np.array([20.0, 40.0, 60.0, 80.0, 75.0])    # one miss of 25

for name, p in [('five small misses', spread), ('one big miss', single)]:
    print('%-18s MAE %5.2f   RMSE %5.2f   R2 %6.3f'
          % (name, mean_absolute_error(truth, p),
             root_mean_squared_error(truth, p), r2_score(truth, p)))

five small misses  MAE  5.00   RMSE  5.00   R2  0.969
one big miss       MAE  5.00   RMSE 11.18   R2  0.844

Identical total absolute error, very different RMSE. If one large failure is much worse for your business than several small ones, a delivery estimate, a capacity forecast, use RMSE. If every unit of error costs the same, use MAE.

R² is not comparable across datasets

R² measures improvement over predicting the mean, so it depends on how variable the test set happens to be. The same model scores higher on a diverse test set than on a homogeneous one. Comparing your R² to a number from a paper on different data tells you nothing.

import numpy as np
from sklearn.metrics import r2_score, root_mean_squared_error

rng = np.random.default_rng(0)

for spread in [5, 20, 60]:
    truth = rng.normal(50, spread, 500)
    pred = truth + rng.normal(0, 4, 500)      # identical error every time
    print('target std %2d   RMSE %.2f   R2 %.4f'
          % (spread, root_mean_squared_error(truth, pred), r2_score(truth, pred)))

target std  5   RMSE 3.76   R2 0.4489
target std 20   RMSE 4.14   R2 0.9581
target std 60   RMSE 3.92   R2 0.9956

The model makes exactly the same size of error in all three cases. R² swings from 0.45 to 0.996. RMSE, which is in the units of the thing you care about, does not move.

MAPE and the zero problem

import numpy as np

def mape(truth, pred):
    return np.mean(np.abs((truth - pred) / truth)) * 100

truth = np.array([100.0, 50.0, 10.0, 0.5])
pred = np.array([110.0, 55.0, 11.0, 0.55])  # every one is 10% high
print('all errors are 10%%: MAPE %.1f%%' % mape(truth, pred))

truth2 = np.array([100.0, 50.0, 10.0, 0.01])
pred2 = np.array([110.0, 55.0, 11.0, 0.5])
print('one near-zero actual: MAPE %.1f%%' % mape(truth2, pred2))

all errors are 10%: MAPE 10.0%
one near-zero actual: MAPE 1232.5%

An error of 0.49 on a true value of 0.01 registers as 4,900 percent and swamps everything else. Never use MAPE on a target that can approach zero, and never on one that can be negative.

Cross-validation, not one split

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')
df = df.dropna(subset=['monthly_charges'])

NUM = ['tenure_months', 'support_calls']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
TARGET = 'monthly_charges'
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, OneHotEncoder

def preprocessor(num=NUM, cat=CAT):
    return ColumnTransformer([
        ('num', Pipeline([('i', SimpleImputer(strategy='median')),
                          ('s', StandardScaler())]), num),
        ('cat', Pipeline([('i', SimpleImputer(strategy='most_frequent')),
                          ('o', OneHotEncoder(handle_unknown='ignore'))]), cat),
    ])
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score, KFold
import numpy as np

X, y = df[NUM + CAT], df[TARGET]
model = Pipeline([('prep', preprocessor()), ('reg', LinearRegression())])

cv = KFold(n_splits=5, shuffle=True, random_state=42)
for metric in ['r2', 'neg_mean_absolute_error', 'neg_root_mean_squared_error']:
    scores = cross_val_score(model, X, y, cv=cv, scoring=metric)
    if metric.startswith('neg_'):
        scores = -scores
    print('%-32s %.4f  +/- %.4f' % (metric, scores.mean(), scores.std()))

r2                               0.8593  +/- 0.0110
neg_mean_absolute_error          7.0785  +/- 0.2229
neg_root_mean_squared_error      8.7868  +/- 0.2839

sklearn negates losses so higher is always better

Every scorer in scikit-learn follows the rule “greater is better”, so errors come back negative. neg_mean_absolute_error of −3.6 means an MAE of 3.6. Forgetting to flip the sign is a common way to conclude your model got worse when it improved.

Day 2 takeaway

MAE treats all errors equally, RMSE punishes large ones, and which you want is a business question. R² is not comparable across datasets because it depends on the test set's spread. MAPE breaks near zero. Cross-validate rather than trusting one split, and remember sklearn negates its error scorers.

Week 04 · Day 3 of 7

Polynomials, Overfitting and Learning Curves

Adding flexibility, watching it backfire, and diagnosing which fix you need

By Pawan Rama Mali1026 words

A straight line cannot fit a curve. The obvious fix, give the model more flexibility, is also the fastest way to make it worse, and watching that happen is the most useful hour of this week.

Add curvature

import numpy as np
from sklearn.preprocessing import PolynomialFeatures

x = np.array([[2.0, 3.0]])
for degree in [1, 2, 3]:
    pf = PolynomialFeatures(degree=degree, include_bias=False)
    out = pf.fit_transform(x)
    print('degree %d -> %d features: %s'
          % (degree, out.shape[1], list(pf.get_feature_names_out())))

degree 1 -> 2 features: ['x0', 'x1']
degree 2 -> 5 features: ['x0', 'x1', 'x0^2', 'x0 x1', 'x1^2']
degree 3 -> 9 features: ['x0', 'x1', 'x0^2', 'x0 x1', 'x1^2', 'x0^3', 'x0^2 x1', 'x0 x1^2', 'x1^3']

Degree 2 adds the squares and the interaction x0 x1. That interaction term is often the valuable one: it lets the effect of one feature depend on another, which a plain linear model cannot express.

The feature count explodes

Ten features at degree 2 gives 65 columns; at degree 3, 285. At degree 5 you have 3,003 columns from ten. The model gains enough freedom to fit the noise exactly, and the number of rows you would need to constrain it grows just as fast.

Watch it overfit

import numpy as np
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline
from sklearn.metrics import root_mean_squared_error

rng = np.random.default_rng(0)
n = 40
x = rng.uniform(0, 6, n).reshape(-1, 1)
y = np.sin(1.5 * x.ravel()) + 0.3 * x.ravel() + rng.normal(0, 0.3, n)

x_val = rng.uniform(0, 6, 400).reshape(-1, 1)
y_val = np.sin(1.5 * x_val.ravel()) + 0.3 * x_val.ravel() + rng.normal(0, 0.3, 400)

print('%6s %10s %10s' % ('degree', 'train', 'validation'))
for degree in [1, 2, 3, 5, 9, 15, 25]:
    m = make_pipeline(PolynomialFeatures(degree), LinearRegression()).fit(x, y)
    tr = root_mean_squared_error(y, m.predict(x))
    va = root_mean_squared_error(y_val, m.predict(x_val))
    print('%6d %10.4f %10.4f' % (degree, tr, va))

degree      train validation
   0.7065     0.7743
   0.5418     0.6249
   0.5140     0.6155
   0.3146     0.3479
   0.3019     0.3262
   0.2953     0.3388
   0.4085     0.4709

Overfitting: The model has fitted the noise in the training rows rather than the pattern behind them. Recognised by training error continuing to fall while validation error rises. The gap between the two is the diagnostic, not the level of either.

Training error falls as flexibility rises, in principle it must, since a richer model can always fit the given points better. Validation error bottoms out at degree 9 and then climbs. The best degree is where validation is lowest, and nothing in the training column would have told you where that was.

Look at degree 25, though: training error goes up, from 0.2953 to 0.4085. That is not overfitting, it is arithmetic failing. Raising a value near 6 to the 25th power produces columns spanning twenty orders of magnitude: the matrix becomes so ill-conditioned that the solver cannot work accurately, and the fit degrades. A training error that rises with complexity is a numerical warning, not a statistical one.

The learning curve: more data, or a better model?

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')
df = df.dropna(subset=['monthly_charges'])

NUM = ['tenure_months', 'support_calls']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
TARGET = 'monthly_charges'
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, OneHotEncoder

def preprocessor(num=NUM, cat=CAT):
    return ColumnTransformer([
        ('num', Pipeline([('i', SimpleImputer(strategy='median')),
                          ('s', StandardScaler())]), num),
        ('cat', Pipeline([('i', SimpleImputer(strategy='most_frequent')),
                          ('o', OneHotEncoder(handle_unknown='ignore'))]), cat),
    ])
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import learning_curve
import numpy as np

X, y = df[NUM + CAT], df[TARGET]
model = Pipeline([('prep', preprocessor()), ('reg', LinearRegression())])

sizes, train_scores, val_scores = learning_curve(
    model, X, y, cv=5, scoring='neg_root_mean_squared_error',
    train_sizes=np.linspace(0.1, 1.0, 6), random_state=42)

print('%8s %10s %10s %8s' % ('rows', 'train', 'val', 'gap'))
for n, tr, va in zip(sizes, -train_scores.mean(axis=1), -val_scores.mean(axis=1)):
    print('%8d %10.4f %10.4f %8.4f' % (n, tr, va, va - tr))

    rows      train        val      gap
   8.3919     9.3860   0.9941
   8.6457     8.8709   0.2252
   8.7090     8.8098   0.1008
   8.7100     8.8003   0.0903
   8.7191     8.7845   0.0654
   8.7448     8.7775   0.0327

How to read a learning curve

The two curves converging to a low error means you are done. Converging to a high error means high bias: more data will not help, you need a better model. A persistent gap means high variance: more data will help, or regularise. Here the gap closes to almost nothing at a low error, so this model is well matched to this problem.

Finding the right complexity

import numpy as np
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import cross_val_score

rng = np.random.default_rng(0)
x = rng.uniform(0, 6, 120).reshape(-1, 1)
y = np.sin(1.5 * x.ravel()) + 0.3 * x.ravel() + rng.normal(0, 0.3, 120)

best = None
for degree in range(1, 16):
    m = make_pipeline(PolynomialFeatures(degree), LinearRegression())
    rmse = -cross_val_score(m, x, y, cv=5,
                            scoring='neg_root_mean_squared_error').mean()
    flag = ''
    if best is None or rmse < best[1]:
        best, flag = (degree, rmse), '  <- best so far'
    if degree <= 8 or flag:
        print('degree %2d  cv RMSE %8.4f%s' % (degree, rmse, flag))

print('\nchosen degree: %d' % best[0])

degree  1  cv RMSE   0.7715  <- best so far
degree  2  cv RMSE   0.6255  <- best so far
degree  3  cv RMSE   0.6250  <- best so far
degree  4  cv RMSE   0.3330  <- best so far
degree  5  cv RMSE   0.3355
degree  6  cv RMSE   0.3052  <- best so far
degree  7  cv RMSE   0.3086
degree  8  cv RMSE   0.3116

chosen degree: 6

Cross-validation picks the complexity for you, and it does so without ever looking at the test set. That is the pattern for every hyperparameter in this course.

Day 3 takeaway

Polynomial features add curvature and interactions, and the column count explodes with degree. Overfitting is diagnosed by the gap between training and validation error, never by training error alone. A learning curve tells you whether more data or a better model is the answer. Choose complexity by cross-validation.

Week 04 · Day 4 of 7

Ridge Regression

Keeping a flexible model honest with an L2 penalty

By Pawan Rama Mali775 words

Yesterday's fix for overfitting was to use a simpler model. Regularisation is the better fix: keep the flexibility, but penalise the model for using it.

The idea

Ridge regression: Ordinary least squares plus a penalty on the sum of the squared coefficients. The model now minimises error plus alpha times that penalty, so a coefficient has to earn its size by reducing error more than it adds to the penalty.

import numpy as np
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.pipeline import make_pipeline
from sklearn.metrics import root_mean_squared_error

rng = np.random.default_rng(0)
x = rng.uniform(0, 6, 40).reshape(-1, 1)
y = np.sin(1.5 * x.ravel()) + 0.3 * x.ravel() + rng.normal(0, 0.3, 40)
x_val = rng.uniform(0, 6, 400).reshape(-1, 1)
y_val = np.sin(1.5 * x_val.ravel()) + 0.3 * x_val.ravel() + rng.normal(0, 0.3, 400)

plain = make_pipeline(PolynomialFeatures(15), StandardScaler(),
                      LinearRegression()).fit(x, y)
ridge = make_pipeline(PolynomialFeatures(15), StandardScaler(),
                      Ridge(alpha=1.0)).fit(x, y)

for name, m in [('plain deg 15', plain), ('ridge deg 15', ridge)]:
    coefs = m[-1].coef_
    print('%-14s val RMSE %7.4f   largest |coef| %12.2f'
          % (name, root_mean_squared_error(y_val, m.predict(x_val)),
             np.abs(coefs).max()))

plain deg 15   val RMSE  0.3480   largest |coef| 1094963376.11
ridge deg 15   val RMSE  0.5461   largest |coef|         0.41

Same fifteen-degree polynomial, same forty rows. The unpenalised version produces enormous coefficients that cancel each other out to thread through every training point. Ridge keeps them small, and generalises far better.

Ridge requires scaled features

The penalty is on the coefficients, and a coefficient's size depends on its feature's units. A feature measured in pounds gets a small coefficient and is barely penalised; the same feature in thousands of pounds gets a large one and is crushed. Without scaling, the penalty lands arbitrarily. Always put a scaler in front, and note that the intercept is never penalised, which is why it is fitted separately.

What alpha does

import numpy as np
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import Ridge
from sklearn.pipeline import make_pipeline
from sklearn.metrics import root_mean_squared_error

rng = np.random.default_rng(0)
x = rng.uniform(0, 6, 40).reshape(-1, 1)
y = np.sin(1.5 * x.ravel()) + 0.3 * x.ravel() + rng.normal(0, 0.3, 40)
x_val = rng.uniform(0, 6, 400).reshape(-1, 1)
y_val = np.sin(1.5 * x_val.ravel()) + 0.3 * x_val.ravel() + rng.normal(0, 0.3, 400)

print('%10s %10s %10s %14s' % ('alpha', 'train', 'val', 'largest |coef|'))
for alpha in [0.0001, 0.01, 1, 100, 10000]:
    m = make_pipeline(PolynomialFeatures(15), StandardScaler(),
                      Ridge(alpha=alpha)).fit(x, y)
    print('%10s %10.4f %10.4f %14.2f'
          % (alpha, root_mean_squared_error(y, m.predict(x)),
             root_mean_squared_error(y_val, m.predict(x_val)),
             np.abs(m[-1].coef_).max()))

     alpha      train        val largest |coef|
    0.0001     0.3061     0.3248          18.87
      0.01     0.3421     0.3808           6.58
         1     0.4653     0.5461           0.41
       100     0.5946     0.6756           0.07
     10000     0.9418     0.9272           0.00

At tiny alpha you have plain least squares and it overfits. At enormous alpha every coefficient is crushed toward zero and the model underfits, predicting close to a constant. The useful values are in between, and you find them by cross-validation.

RidgeCV finds alpha for you

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')
df = df.dropna(subset=['monthly_charges'])

NUM = ['tenure_months', 'support_calls']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
TARGET = 'monthly_charges'
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, OneHotEncoder

def preprocessor(num=NUM, cat=CAT):
    return ColumnTransformer([
        ('num', Pipeline([('i', SimpleImputer(strategy='median')),
                          ('s', StandardScaler())]), num),
        ('cat', Pipeline([('i', SimpleImputer(strategy='most_frequent')),
                          ('o', OneHotEncoder(handle_unknown='ignore'))]), cat),
    ])
from sklearn.linear_model import RidgeCV
from sklearn.model_selection import cross_val_score
import numpy as np

X, y = df[NUM + CAT], df[TARGET]
alphas = np.logspace(-3, 3, 25)

model = Pipeline([('prep', preprocessor()),
                  ('reg', RidgeCV(alphas=alphas))]).fit(X, y)
print('chosen alpha: %.4f' % model.named_steps['reg'].alpha_)

rmse = -cross_val_score(model, X, y, cv=5,
                        scoring='neg_root_mean_squared_error').mean()
print('cv RMSE %.4f' % rmse)

chosen alpha: 1.7783
cv RMSE 8.7776

Search alpha on a log scale

The useful range spans orders of magnitude, so np.logspace(-3, 3, 25) is the right shape of grid and np.linspace is not. If the chosen value lands at either end of your grid, the grid was too narrow. Widen it and search again.

Day 4 takeaway

Ridge adds a penalty on squared coefficients, which lets you keep a flexible model without letting it produce wild parameters. It requires scaled features, because the penalty is unit-dependent. Alpha trades fit against simplicity; search it on a log scale and widen the grid if the winner is at an edge.

Week 04 · Day 5 of 7

Lasso and ElasticNet

L1 as feature selection, its instability, and the mix that fixes it

By Pawan Rama Mali827 words

Ridge shrinks every coefficient toward zero but never quite to it. Lasso takes them all the way, which turns regularisation into feature selection.

L1 against L2

	Ridge (L2)	Lasso (L1)
Penalty	Sum of squared coefficients	Sum of absolute coefficients
Effect	Shrinks all, eliminates none	Sets some exactly to zero
Correlated features	Splits the weight between them	Picks one arbitrarily, drops the rest
Use when	All features plausibly matter	You expect most to be irrelevant

import numpy as np
from sklearn.linear_model import Ridge, Lasso
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline

rng = np.random.default_rng(0)
n, p = 200, 20
X = rng.normal(size=(n, p))
# Only the first three columns matter; the other seventeen are noise.
true_coef = np.zeros(p)
true_coef[:3] = [4.0, -3.0, 2.0]
y = X @ true_coef + rng.normal(0, 1.0, n)

for name, model in [('ridge', Ridge(alpha=1.0)), ('lasso', Lasso(alpha=0.1))]:
    m = make_pipeline(StandardScaler(), model).fit(X, y)
    coef = m[-1].coef_
    print('%-6s first three %s   exactly zero: %d of %d'
          % (name, np.round(coef[:3], 2), (coef == 0).sum(), p))

ridge  first three [ 3.84 -2.94  1.92]   exactly zero: 0 of 20
lasso  first three [ 3.72 -2.87  1.77]   exactly zero: 15 of 20

Ridge leaves all twenty coefficients non-zero, including seventeen that should not be there. Lasso zeroes most of the noise columns outright, giving you a model you can read.

The path from full model to empty one

import numpy as np
from sklearn.linear_model import Lasso
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline

rng = np.random.default_rng(0)
X = rng.normal(size=(200, 20))
true_coef = np.zeros(20)
true_coef[:3] = [4.0, -3.0, 2.0]
y = X @ true_coef + rng.normal(0, 1.0, 200)

print('%8s %10s %s' % ('alpha', 'kept', 'which'))
for alpha in [0.001, 0.05, 0.2, 0.5, 1.0, 3.0]:
    m = make_pipeline(StandardScaler(), Lasso(alpha=alpha, max_iter=5000)).fit(X, y)
    kept = np.flatnonzero(m[-1].coef_)
    print('%8s %10d %s' % (alpha, len(kept), list(kept[:6])))

   alpha       kept which
001         20 [np.int64(0), np.int64(1), np.int64(2), np.int64(3), np.int64(4), np.int64(5)]
05          9 [np.int64(0), np.int64(1), np.int64(2), np.int64(5), np.int64(7), np.int64(9)]
2          3 [np.int64(0), np.int64(1), np.int64(2)]
5          3 [np.int64(0), np.int64(1), np.int64(2)]
0          3 [np.int64(0), np.int64(1), np.int64(2)]
0          1 [np.int64(0)]

As alpha rises the model discards features, and the last three standing are exactly the three that generated the data. That is the best possible outcome and it will not always happen, but it shows what lasso is for.

Lasso is unstable when features are correlated

Given two nearly identical columns, lasso keeps one and zeroes the other, and which one it keeps can flip on a slightly different sample. Do not read “lasso dropped this feature” as “this feature does not matter”. It may simply be standing next to a twin.

import numpy as np
from sklearn.linear_model import Lasso
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline

rng = np.random.default_rng(1)
n = 150
a = rng.normal(size=n)
twin = a + rng.normal(0, 0.05, n)          # almost the same column
noise = rng.normal(size=(n, 5))
y = 3 * a + rng.normal(0, 0.5, n)

for seed in range(4):
    idx = np.random.default_rng(seed).choice(n, 120, replace=False)
    X = np.column_stack([a, twin, noise])[idx]
    m = make_pipeline(StandardScaler(), Lasso(alpha=0.15)).fit(X, y[idx])
    print('sample %d: coef on a %6.3f   coef on its twin %6.3f'
          % (seed, m[-1].coef_[0], m[-1].coef_[1]))

sample 0: coef on a  1.966   coef on its twin  0.229
sample 1: coef on a  2.102   coef on its twin  0.468
sample 2: coef on a  2.379   coef on its twin  0.000
sample 3: coef on a  2.278   coef on its twin  0.064

Same underlying truth, four resamples, and the weight jumps between the two twins. Ridge would have split it evenly and stably between them.

ElasticNet: both penalties

import numpy as np
from sklearn.linear_model import ElasticNetCV
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline

rng = np.random.default_rng(1)
n = 150
a = rng.normal(size=n)
twin = a + rng.normal(0, 0.05, n)
noise = rng.normal(size=(n, 5))
X = np.column_stack([a, twin, noise])
y = 3 * a + rng.normal(0, 0.5, n)

m = make_pipeline(StandardScaler(),
                  ElasticNetCV(l1_ratio=[0.1, 0.5, 0.9, 1.0], cv=5,
                               random_state=0, max_iter=5000)).fit(X, y)
net = m[-1]
print('chosen l1_ratio %.2f   alpha %.4f' % (net.l1_ratio_, net.alpha_))
print('coef on a %.3f, on its twin %.3f' % (net.coef_[0], net.coef_[1]))
print('noise coefficients:', np.round(net.coef_[2:], 3))

chosen l1_ratio 0.90   alpha 0.0336
coef on a 1.589, on its twin 1.007
noise coefficients: [0.    0.    0.047 0.015 0.003]

l1_ratio mixes the two: 1.0 is pure lasso, 0.0 pure ridge. ElasticNet keeps lasso's ability to discard irrelevant columns while sharing weight between correlated ones rather than choosing arbitrarily.

Day 5 takeaway

L2 shrinks, L1 selects. Lasso produces readable models by zeroing coefficients, but distributes weight arbitrarily among correlated features and is unstable across samples as a result. ElasticNet mixes both and is the safer default when you have correlated columns and still want selection.

Week 04 · Day 6 of 7

Multicollinearity and Target Leakage

The failure that corrupts coefficients, and the one that fakes success

By Pawan Rama Mali1084 words

Two failures that both look like success. One inflates your test score to something impossible; the other makes your coefficients meaningless while the score stays fine.

Multicollinearity

Multicollinearity: Two or more features carrying nearly the same information. The model can trade weight between them freely without changing its predictions, so the individual coefficients become large, unstable and uninterpretable, while the overall fit stays perfectly good.

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')
df = df.dropna(subset=['monthly_charges'])

NUM = ['tenure_months', 'support_calls']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
TARGET = 'monthly_charges'
d = df.dropna(subset=['total_charges'])
print(d[['tenure_months', 'monthly_charges', 'total_charges']].corr().round(3))

                 tenure_months  monthly_charges  total_charges
tenure_months            1.000            0.001          0.828
monthly_charges          0.001            1.000          0.450
total_charges            0.828            0.450          1.000

import numpy as np
from sklearn.linear_model import LinearRegression, Ridge

rng = np.random.default_rng(0)
n = 300
x1 = rng.normal(0, 1, n)
x2 = x1 + rng.normal(0, 0.02, n)     # correlation about 0.9998
y = 2 * x1 + rng.normal(0, 0.3, n)

print('correlation between the two features: %.4f' % np.corrcoef(x1, x2)[0, 1])
print('\nrefitting on four different subsamples:')
print('%12s %10s %10s %10s' % ('', 'coef x1', 'coef x2', 'sum'))
for seed in range(4):
    idx = np.random.default_rng(seed).choice(n, 250, replace=False)
    X = np.column_stack([x1, x2])[idx]
    m = LinearRegression().fit(X, y[idx])
    print('%12s %10.3f %10.3f %10.3f'
          % ('sample %d' % seed, m.coef_[0], m.coef_[1], m.coef_.sum()))

correlation between the two features: 0.9998

refitting on four different subsamples:
                coef x1    coef x2        sum
    sample 0      1.421      0.594      2.015
    sample 1      0.799      1.212      2.011
    sample 2      0.950      1.060      2.010
    sample 3      1.364      0.640      2.004

The individual coefficients swing violently, and their sum stays near 2, the truth. The model knows the total effect perfectly well; it simply cannot tell which of two identical columns deserves the credit. Reporting either coefficient on its own would be meaningless.

Measure it with the variance inflation factor

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')
df = df.dropna(subset=['monthly_charges'])

NUM = ['tenure_months', 'support_calls']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
TARGET = 'monthly_charges'
import numpy as np
from sklearn.linear_model import LinearRegression

d = df.dropna(subset=['total_charges'])
cols = ['tenure_months', 'support_calls', 'monthly_charges', 'total_charges']
M = d[cols].to_numpy(dtype=float)

print('%-18s %8s' % ('feature', 'VIF'))
for i, name in enumerate(cols):
    others = np.delete(M, i, axis=1)
    r2 = LinearRegression().fit(others, M[:, i]).score(others, M[:, i])
    print('%-18s %8.2f' % (name, 1 / (1 - r2)))

feature                 VIF
tenure_months          7.03
support_calls          1.02
monthly_charges        2.79
total_charges          8.82

How to read a VIF

VIF is how much the variance of a coefficient is inflated by its correlation with the others. Below 5 is fine, above 10 is a problem worth acting on. The fixes, in order of preference: drop one of the pair, combine them into a single meaningful feature, or use ridge, which handles collinearity gracefully by design.

Leakage, which is the more dangerous of the two

Target leakage: A feature that contains information about the target that would not be available at prediction time. It produces excellent test scores and a model that fails completely in production, because the leaking column is the answer in disguise.

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')
df = df.dropna(subset=['monthly_charges'])

NUM = ['tenure_months', 'support_calls']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
TARGET = 'monthly_charges'
d = df.dropna(subset=['total_charges']).copy()
d = d[d['tenure_months'] > 0]
d['implied_monthly'] = d['total_charges'] / d['tenure_months']

print('correlation with the target: %.4f'
      % d['implied_monthly'].corr(d[TARGET]))
print('\nfirst five rows:')
print(d[['tenure_months', 'total_charges', 'implied_monthly', TARGET]]
      .head().round(2).to_string(index=False))

correlation with the target: 0.9956

first five rows:
 tenure_months  total_charges  implied_monthly  monthly_charges
            10         904.14            90.41            90.21
             8         621.20            77.65            81.45
             3         181.14            60.38            57.29
             3          44.38            14.79            15.00
            21         564.82            26.90            25.85

total_charges divided by tenure_months is, almost exactly, monthly_charges, because that is how the bill was calculated. A perfectly sensible-looking engineered feature reconstructs the answer.

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')
df = df.dropna(subset=['monthly_charges'])

NUM = ['tenure_months', 'support_calls']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
TARGET = 'monthly_charges'
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, OneHotEncoder

def preprocessor(num=NUM, cat=CAT):
    return ColumnTransformer([
        ('num', Pipeline([('i', SimpleImputer(strategy='median')),
                          ('s', StandardScaler())]), num),
        ('cat', Pipeline([('i', SimpleImputer(strategy='most_frequent')),
                          ('o', OneHotEncoder(handle_unknown='ignore'))]), cat),
    ])
from sklearn.model_selection import train_test_split

X = df[NUM + CAT]
y = df[TARGET]
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.25, random_state=42)
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score

honest = Pipeline([('prep', preprocessor()),
                   ('reg', LinearRegression())]).fit(X_tr, y_tr)
print('honest model      R2 %.4f' % r2_score(y_te, honest.predict(X_te)))

leaky_df = df.dropna(subset=['total_charges']).copy()
leaky_df = leaky_df[leaky_df['tenure_months'] > 0]
leaky_df['implied_monthly'] = (leaky_df['total_charges']
                               / leaky_df['tenure_months'])
LNUM = NUM + ['implied_monthly']
Xl, yl = leaky_df[LNUM + CAT], leaky_df[TARGET]
Xl_tr, Xl_te, yl_tr, yl_te = train_test_split(Xl, yl, test_size=0.25,
                                              random_state=42)
leaky = Pipeline([('prep', preprocessor(num=LNUM)),
                  ('reg', LinearRegression())]).fit(Xl_tr, yl_tr)
print('with the leak     R2 %.4f' % r2_score(yl_te, leaky.predict(Xl_te)))

honest model      R2 0.8498
with the leak     R2 0.9920

An R² that jumps to 0.99 is not good news

It is the strongest available signal that a feature knows the answer. The instinct on seeing a score like that should be suspicion, not celebration. Ask of every feature: would I have this value, for a new customer, at the moment I need the prediction? If the answer is no, or “only after the outcome is known”, it leaks.

Where leakage hides

Aggregates computed over the whole dataset: a per-category mean that included the test rows.
Anything recorded after the event: a cancellation reason code, when predicting cancellation.
Identifiers correlated with time: sequential IDs encode when a row was created.
Duplicated rows across the split: exactly the 138 you removed in week 1.
Preprocessing fitted before the split: the reason every scaler in this course sits inside a pipeline.

Day 6 takeaway

Multicollinearity destabilises coefficients while leaving predictions intact, measure it with VIF and prefer ridge if you keep the correlated columns. Leakage inflates the score itself and is far more dangerous. Test every feature against one question: would I actually have this value at prediction time?

Week 04 · Day 7 of 7

A Complete Regression Project

Every decision of the week, applied end to end

By Pawan Rama Mali1282 words

A complete regression, start to finish, with every decision this week justified.

1. Frame and split

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')
df = df.dropna(subset=['monthly_charges'])

NUM = ['tenure_months', 'support_calls']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
TARGET = 'monthly_charges'
from sklearn.model_selection import train_test_split

X = df[NUM + CAT]
y = df[TARGET]
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.25, random_state=42)
print('train %d rows, test %d rows' % (len(X_tr), len(X_te)))
print('target mean: train %.2f, test %.2f' % (y_tr.mean(), y_te.mean()))

train 2126 rows, test 709 rows
target mean: train 59.48, test 59.20

2. Establish the floor

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')
df = df.dropna(subset=['monthly_charges'])

NUM = ['tenure_months', 'support_calls']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
TARGET = 'monthly_charges'
from sklearn.model_selection import train_test_split

X = df[NUM + CAT]
y = df[TARGET]
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.25, random_state=42)
from sklearn.dummy import DummyRegressor
from sklearn.metrics import root_mean_squared_error, r2_score

base = DummyRegressor(strategy='mean').fit(X_tr, y_tr)
pred = base.predict(X_te)
print('predicting the mean: RMSE %.3f   R2 %.4f'
      % (root_mean_squared_error(y_te, pred), r2_score(y_te, pred)))

predicting the mean: RMSE 23.310 R2 -0.0001

R² of essentially zero, by construction. Every model from here has to beat an RMSE of about 23.4.

3. Compare candidates by cross-validation

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')
df = df.dropna(subset=['monthly_charges'])

NUM = ['tenure_months', 'support_calls']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
TARGET = 'monthly_charges'
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, OneHotEncoder

def preprocessor(num=NUM, cat=CAT):
    return ColumnTransformer([
        ('num', Pipeline([('i', SimpleImputer(strategy='median')),
                          ('s', StandardScaler())]), num),
        ('cat', Pipeline([('i', SimpleImputer(strategy='most_frequent')),
                          ('o', OneHotEncoder(handle_unknown='ignore'))]), cat),
    ])
from sklearn.model_selection import train_test_split

X = df[NUM + CAT]
y = df[TARGET]
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.25, random_state=42)
from sklearn.linear_model import LinearRegression, RidgeCV, LassoCV
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import cross_val_score, KFold
import numpy as np

cv = KFold(n_splits=5, shuffle=True, random_state=42)
candidates = {
    'linear    ': LinearRegression(),
    'ridge     ': RidgeCV(alphas=np.logspace(-3, 3, 25)),
    'lasso     ': LassoCV(cv=5, random_state=0, max_iter=5000),
    'forest    ': RandomForestRegressor(n_estimators=200, random_state=42),
}
for name, reg in candidates.items():
    pipe = Pipeline([('prep', preprocessor()), ('reg', reg)])
    s = -cross_val_score(pipe, X_tr, y_tr, cv=cv,
                         scoring='neg_root_mean_squared_error')
    print('%s cv RMSE %.4f  +/- %.4f' % (name, s.mean(), s.std()))

linear     cv RMSE 8.7237  +/- 0.1761
ridge      cv RMSE 8.7234  +/- 0.1735
lasso      cv RMSE 8.7068  +/- 0.1714
forest     cv RMSE 9.9159  +/- 0.2897

Cross-validate on the training set only

The test set has not been touched yet and must not be, until the model is chosen. Every comparison, every hyperparameter, every feature decision comes from cross-validation within the training data. The test set is spent the first time you look at it.

4. Fit the winner and score once

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')
df = df.dropna(subset=['monthly_charges'])

NUM = ['tenure_months', 'support_calls']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
TARGET = 'monthly_charges'
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, OneHotEncoder

def preprocessor(num=NUM, cat=CAT):
    return ColumnTransformer([
        ('num', Pipeline([('i', SimpleImputer(strategy='median')),
                          ('s', StandardScaler())]), num),
        ('cat', Pipeline([('i', SimpleImputer(strategy='most_frequent')),
                          ('o', OneHotEncoder(handle_unknown='ignore'))]), cat),
    ])
from sklearn.model_selection import train_test_split

X = df[NUM + CAT]
y = df[TARGET]
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.25, random_state=42)
from sklearn.linear_model import RidgeCV
from sklearn.metrics import (root_mean_squared_error, mean_absolute_error,
                             r2_score)
import numpy as np

final = Pipeline([('prep', preprocessor()),
                  ('reg', RidgeCV(alphas=np.logspace(-3, 3, 25)))])
final.fit(X_tr, y_tr)
pred = final.predict(X_te)

print('TEST SET')
print('  RMSE %.4f' % root_mean_squared_error(y_te, pred))
print('  MAE  %.4f' % mean_absolute_error(y_te, pred))
print('  R2   %.4f' % r2_score(y_te, pred))
print('\nagainst a mean-only RMSE of %.3f' % y_te.std())

TEST SET
  RMSE 9.0293
  MAE  7.2018
  R2   0.8499

against a mean-only RMSE of 23.325

5. Explain it

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')
df = df.dropna(subset=['monthly_charges'])

NUM = ['tenure_months', 'support_calls']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
TARGET = 'monthly_charges'
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, OneHotEncoder

def preprocessor(num=NUM, cat=CAT):
    return ColumnTransformer([
        ('num', Pipeline([('i', SimpleImputer(strategy='median')),
                          ('s', StandardScaler())]), num),
        ('cat', Pipeline([('i', SimpleImputer(strategy='most_frequent')),
                          ('o', OneHotEncoder(handle_unknown='ignore'))]), cat),
    ])
from sklearn.model_selection import train_test_split

X = df[NUM + CAT]
y = df[TARGET]
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.25, random_state=42)
from sklearn.linear_model import RidgeCV
import numpy as np
import pandas as pd

final = Pipeline([('prep', preprocessor()),
                  ('reg', RidgeCV(alphas=np.logspace(-3, 3, 25)))]).fit(X_tr, y_tr)

names = final.named_steps['prep'].get_feature_names_out()
coefs = pd.Series(final.named_steps['reg'].coef_, index=names)
print(coefs.reindex(coefs.abs().sort_values(ascending=False).index)
      .head(8).round(3).to_string())

cat__internet_service_No internet      -30.107
cat__internet_service_Fibre optic       27.346
cat__internet_service_DSL                2.762
cat__contract_Month-To-Month             0.322
cat__payment_method_Electronic check    -0.304
num__tenure_months                       0.299
cat__contract_One Year                  -0.290
cat__payment_method_Bank transfer        0.215

Fibre optic adds about 27 against the baseline and having no internet subtracts about 30, while contract type and payment method barely move the price at all. The model has recovered the pricing structure, which is exactly what the generator in week 1 put there, and note that the features which drive churn are almost irrelevant to price. Different target, different important features.

6. Check the residuals one last time

import numpy as np
import pandas as pd

df = pd.read_csv('customers.csv').drop_duplicates().copy()
df['contract'] = df['contract'].str.strip().str.title()
df['total_charges'] = pd.to_numeric(df['total_charges'], errors='coerce')
df = df.dropna(subset=['monthly_charges'])

NUM = ['tenure_months', 'support_calls']
CAT = ['contract', 'internet_service', 'payment_method', 'has_dependents']
TARGET = 'monthly_charges'
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, OneHotEncoder

def preprocessor(num=NUM, cat=CAT):
    return ColumnTransformer([
        ('num', Pipeline([('i', SimpleImputer(strategy='median')),
                          ('s', StandardScaler())]), num),
        ('cat', Pipeline([('i', SimpleImputer(strategy='most_frequent')),
                          ('o', OneHotEncoder(handle_unknown='ignore'))]), cat),
    ])
from sklearn.model_selection import train_test_split

X = df[NUM + CAT]
y = df[TARGET]
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.25, random_state=42)
from sklearn.linear_model import RidgeCV
import numpy as np
import pandas as pd

final = Pipeline([('prep', preprocessor()),
                  ('reg', RidgeCV(alphas=np.logspace(-3, 3, 25)))]).fit(X_tr, y_tr)
pred = final.predict(X_te)
resid = y_te - pred

print('residual mean %.4f, std %.4f, skew %.4f'
      % (resid.mean(), resid.std(), resid.skew()))
worst = resid.abs().nlargest(3).index
print('\nworst three predictions:')
print(pd.DataFrame({'actual': y_te.loc[worst].round(2),
                    'predicted': pred[[list(y_te.index).index(i) for i in worst]].round(2),
                    'service': df.loc[worst, 'internet_service']}).to_string())

residual mean -0.4340, std 9.0253, skew 0.0221

worst three predictions:
      actual  predicted      service
3076   48.61      79.22  Fibre optic
767    52.42      79.49  Fibre optic
132   105.96      78.98  Fibre optic

Your assignment

Rebuild this predicting total_charges instead. You will find an R² above 0.99, because tenure multiplied by monthly charges is the total. Then remove both of those features and try again. Write down which version you would deploy, and why the first one is worthless despite scoring better.

Day 7 takeaway

Split first, establish the floor a mean-only model sets, compare candidates by cross-validation inside the training set, then touch the test set exactly once. Read the coefficients to check the model learned something sensible, and read the residuals to see what it still cannot explain.