Feature Engineering · Notebook 10, the final articleCommon Pitfalls and Best Practices
Notebook 9 settled the leakage rule: every transformer sits inside a pipeline, each one refits on the training fold, and predict can never quietly relearn anything. This notebook exists because that's true and still doesn't prevent any of the failures below. The first three pass straight through a correct pipeline, and cross-validation can't see them. The rest are about what happens around the model: too many columns, a world that changes, and a threshold nobody chose.
Start hereIntroduction#
A pipeline controls when a transformer is allowed to learn. It says nothing about whether a column should have been collected at all, whether it will actually exist when you need a prediction, or whether the score you're reading has been used up by the process that produced it. Those questions are outside what a pipeline can enforce, and they're what this final notebook is about.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import (StandardScaler, OneHotEncoder, OrdinalEncoder,
TargetEncoder, FunctionTransformer)
from sklearn.impute import SimpleImputer
from sklearn.linear_model import LogisticRegression, Ridge
from sklearn.model_selection import cross_val_score, StratifiedKFold, train_test_split, KFold
from sklearn.metrics import roc_auc_score, r2_score, recall_score, f1_score, precision_score
SEED = 42
CV = StratifiedKFold(n_splits=5, shuffle=True, random_state=SEED)
def new_lr():
return LogisticRegression(max_iter=2000, random_state=SEED)
telco = pd.read_csv('datasets/telco_churn.csv')
telco['TotalCharges'] = pd.to_numeric(telco['TotalCharges'], errors='coerce')
telco = telco[telco['tenure'] > 0].reset_index(drop=True)
SERVICES = ['PhoneService', 'MultipleLines', 'InternetService', 'OnlineSecurity',
'OnlineBackup', 'DeviceProtection', 'TechSupport', 'StreamingTV',
'StreamingMovies']
telco['bundle'] = telco[SERVICES].apply(lambda row: '|'.join(row.astype(str)), axis=1)
y = (telco['Churn'] == 'Yes').astype(int)
X = telco.drop(columns=['Churn'])
NUMERIC = ['tenure', 'MonthlyCharges', 'TotalCharges']
NOMINAL = ['gender', 'SeniorCitizen', 'Partner', 'Dependents',
'PaperlessBilling', 'PaymentMethod', 'InternetService']
CONTRACT_ORDER = [['Month-to-month', 'One year', 'Two year']]
def count_services(frame):
return frame.isin(['Yes', 'DSL', 'Fiber optic']).sum(axis=1).to_frame('num_services')
def service_feature_names(transformer, input_features):
return ['num_services']
service_counter = FunctionTransformer(count_services, feature_names_out=service_feature_names)
def make_preprocessor(extra_numeric=()):
"""Notebook 9's preprocessor, with a hook for adding a column to the numeric group."""
return ColumnTransformer([
('num', Pipeline([('impute', SimpleImputer(strategy='median')),
('scale', StandardScaler())]), NUMERIC + list(extra_numeric)),
('svc', Pipeline([('count', service_counter), ('scale', StandardScaler())]), SERVICES),
('ord', OrdinalEncoder(categories=CONTRACT_ORDER), ['Contract']),
('nom', OneHotEncoder(handle_unknown='ignore', drop='if_binary',
sparse_output=False), NOMINAL),
('hi', TargetEncoder(random_state=SEED), ['bundle']),
], remainder='drop', verbose_feature_names_out=False)
def make_pipe(extra_numeric=()):
return Pipeline([('prep', make_preprocessor(extra_numeric)), ('model', new_lr())])
honest = cross_val_score(make_pipe(), X, y, cv=CV, scoring='roc_auc', n_jobs=1)
print('the honest pipeline from notebook 9: %.4f' % honest.mean())
the honest pipeline from notebook 9: 0.8416
0.8416, exactly where the series left off. Every failure in this notebook is measured against this one trusted number. That's why it's worth re-checking first: when something below reports 0.98, you can see right away how big the lie is.
Section 1A column that exists because of the answer#
Target leakage is a column recorded after the outcome happened, or because it happened. The model stops predicting the outcome and starts reading a record of it
Telco doesn't have a column like this, so this section builds one on purpose. Imagine a support system logs a cancellation_ticket when a customer calls to cancel. It's a real field with a real business reason, and 90% of customers who churn have one. But it can't help in production, because by the time it exists, the customer has already decided to leave.
rng = np.random.default_rng(SEED)
X_leaky = X.copy()
X_leaky['cancellation_ticket'] = np.where(y == 1,
rng.binomial(1, 0.90, len(y)), # 90% of churners phone in
rng.binomial(1, 0.02, len(y))) # 2% of stayers phone in and stay
print('ticket rate among churners %.3f, among stayers %.3f'
% (X_leaky.loc[y == 1, 'cancellation_ticket'].mean(),
X_leaky.loc[y == 0, 'cancellation_ticket'].mean()))
leaky = cross_val_score(make_pipe(['cancellation_ticket']), X_leaky, y,
cv=CV, scoring='roc_auc', n_jobs=1)
print()
print('honest pipeline %.4f folds %s' % (honest.mean(), honest.round(4)))
print('one column added %.4f folds %s' % (leaky.mean(), leaky.round(4)))
series_gain = 0.8416 - 0.8087 # what all the real feature work in notebooks 1 to 9 was worth
print('the column is worth %+.4f, which is %.1f times everything notebooks 1 to 9 found put together'
% (leaky.mean() - honest.mean(), (leaky.mean() - honest.mean()) / series_gain))
ticket rate among churners 0.900, among stayers 0.019 honest pipeline 0.8416 folds [0.8517 0.8438 0.8405 0.8372 0.8349] one column added 0.9809 folds [0.9821 0.9771 0.9792 0.9807 0.9853] the column is worth +0.1393, which is 4.2 times everything notebooks 1 to 9 found put together
0.9809, up from 0.8416. That one made-up column is worth +0.1393, which is 4.2 times the 0.0329 this series spent nine notebooks building with real features. In a real project, this is the moment to be suspicious, not proud. The useful question is which of your usual checks would have caught it.
print('check 1: is the cross-validation unstable?')
print(' honest pipeline fold spread %.4f' % honest.std())
print(' leaky pipeline fold spread %.4f' % leaky.std())
print()
print('check 2: does any single column predict the target implausibly well?')
num_services = X[SERVICES].isin(['Yes', 'DSL', 'Fiber optic']).sum(axis=1)
singles = {
'cancellation_ticket': X_leaky['cancellation_ticket'],
'tenure (negated)': -X['tenure'],
'MonthlyCharges': X['MonthlyCharges'],
'num_services': num_services,
}
for name, column in singles.items():
print(' univariate AUC, %-22s %.4f' % (name, roc_auc_score(y, column)))
check 1: is the cross-validation unstable? honest pipeline fold spread 0.0059 leaky pipeline fold spread 0.0028 check 2: does any single column predict the target implausibly well? univariate AUC, cancellation_ticket 0.9406 univariate AUC, tenure (negated) 0.7414 univariate AUC, MonthlyCharges 0.6203 univariate AUC, num_services 0.4967
The first check doesn't just miss it. It points the wrong way. The leaky model's five folds agree with each other more closely than the honest model's do, 0.0028 against 0.0059, because a column that contains the answer works equally well in every fold, and scores close to 1.0 have less room to vary. A fold-stability check would say the broken model is the more trustworthy one.
The second check works, up to a point. A single yes/no column scoring 0.9406 on its own, when the best real column tested reaches 0.7414, is definitely worth investigating. But look at the last row before turning that into a rule. num_services scores 0.4967, a coin flip, yet section 4 shows it adds +0.0082 to the model (notebook 7 measured +0.0109 on its own).
Its relationship with churn is an upside-down U (mostly because of contract mix, as notebook 7 found), and a single-column AUC can't see a relationship that doesn't move in one direction. The same screen that would flag the leak would also throw away one of the more useful real features in the series.
Section 2A feature that will not be there when you need it#
Section 1's column was ruled out because it simply can't exist yet. This one is subtler. A feature can be completely legitimate and honestly linked to the target, but still not be available when you need it
This happens most often with rolling windows, so this section switches to hourly bike rentals, where notebook 5 found strong daily patterns. A short rolling average of recent rentals is an obviously sensible feature. The only question is which hours go into it.
bike = pd.read_csv('datasets/bike_sharing.csv')
bike['dteday'] = pd.to_datetime(bike['dteday'])
bike = bike.sort_values(['dteday', 'hr']).reset_index(drop=True) # nb5: sort by both, the default sort is unstable
bike['year'] = bike['dteday'].dt.year
counts = bike['cnt']
for window in (3, 24):
bike['trailing_%d' % window] = counts.shift(1).rolling(window, min_periods=1).mean()
bike['centred_%d' % window] = counts.rolling(window, center=True, min_periods=1).mean()
print('correlation with the value we are trying to predict:')
for window in (3, 24):
print(' %2dh trailing (shifted first) %.3f %2dh centred %.3f'
% (window, bike['trailing_%d' % window].corr(counts),
window, bike['centred_%d' % window].corr(counts)))
correlation with the value we are trying to predict: 3h trailing (shifted first) 0.670 3h centred 0.977 24h trailing (shifted first) 0.401 24h centred 0.436
trailing shifts first, so the average for an hour uses only earlier hours. centred doesn't shift, and center=True makes it reach into the future on purpose.
BIKE_BASE = ['hr', 'workingday', 'temp', 'hum', 'windspeed', 'weathersit', 'season', 'holiday', 'weekday']
def bike_model(extra):
return Pipeline([
('prep', ColumnTransformer([
('onehot', OneHotEncoder(handle_unknown='ignore', sparse_output=False),
['hr', 'weekday', 'season', 'weathersit']),
('scale', StandardScaler(), ['temp', 'hum', 'windspeed'] + extra),
('raw', 'passthrough', ['workingday', 'holiday'])])),
('model', Ridge(alpha=1.0))])
usable = bike.dropna(subset=['trailing_3']).reset_index(drop=True)
train_2011 = usable[usable['year'] == 2011]
test_2012 = usable[usable['year'] == 2012]
def bike_r2(extra):
fitted = bike_model(extra).fit(train_2011[BIKE_BASE + extra], train_2011['cnt'])
return r2_score(test_2012['cnt'], fitted.predict(test_2012[BIKE_BASE + extra]))
print('train on 2011, predict 2012:')
print(' no window feature at all R2 %.4f' % bike_r2([]))
for window in (3, 24):
print(' %2dh trailing, shifted (legitimate) R2 %.4f' % (window, bike_r2(['trailing_%d' % window])))
print(' %2dh centred (uses future hours) R2 %.4f' % (window, bike_r2(['centred_%d' % window])))
train on 2011, predict 2012: no window feature at all R2 0.4294 3h trailing, shifted (legitimate) R2 0.6774 3h centred (uses future hours) R2 0.9716 24h trailing, shifted (legitimate) R2 0.5636 24h centred (uses future hours) R2 0.6230
Two window sizes, two lessons, and both matter. At three hours, the cheat is huge and obvious. The centred window averages the hour before, the hour itself and the hour after, so it contains the answer directly. It correlates 0.977 with the value being predicted and pushes R-squared from 0.6774 up to 0.9716. Nobody would ship that model. A 0.977 correlation with your own target makes anyone look twice.
At twenty-four hours, the cheat is almost invisible, but just as real. Correlation moves from 0.401 to 0.436, a difference you'd never notice in a correlation table, and R-squared moves from 0.5636 to 0.6230. That's a gain of 0.0594 that will disappear on the model's first real day in production, and nothing in a normal workflow flags it.
Be careful not to draw the wrong lesson. This doesn't show that subtle leaks are worse than obvious ones: the three-hour cheat is both louder and bigger. It shows that the dangerous leak is the one that survives review, and only one of these two would. A 0.977 correlation gets caught in the first code review. A move from 0.401 to 0.436 gets caught by nobody, and it's still 0.0594 of R-squared that won't exist on day one.
One more point about "legitimate." The trailing window uses the previous hour's rental count, so it's only available if you predict one hour ahead. If the model has to forecast tomorrow, even the trailing window needs a shift of at least 24 hours. Whether a feature is available depends on when the prediction is made, not just on the code.
A smaller detail: rolling(3) counts rows, not hours. The dataset has 17,379 of the 17,544 possible hours, so a few windows span a gap. A time-based window such as rolling('3h') on a datetime index avoids that. It doesn't change the lesson here.
Section 3The validation set you used up#
This one is different. There's no bad column, no future information and no fitting on the wrong rows. Every single step is correct. The failure comes from repeating a correct step
Notebook 9 measured the start of this: quoting the best of 81 grid-search combinations inflated the score by only 0.0001, because all 81 candidates agreed closely. This section looks at the case where they don't, which is what happens when a person, not a grid, does the choosing.
Here's the loop almost everyone runs without thinking of it as a loop. Hold out a validation set. Think of a feature. Add it. Did validation go up? Keep it. Think of another one.
To measure what this loop does on its own, the person is replaced with 200 columns of pure random noise. So any improvement it reports is an illusion. And since a validation set used for decisions is no longer a fair test, a third split is held back that nothing in this experiment can touch.
X_fit, X_rest, y_fit, y_rest = train_test_split(X, y, test_size=0.4, random_state=SEED, stratify=y)
X_val, X_test, y_val, y_test = train_test_split(X_rest, y_rest, test_size=0.5,
random_state=SEED, stratify=y_rest)
print('fit %d rows validation %d rows untouched test %d rows' % (len(X_fit), len(X_val), len(X_test)))
prep = make_preprocessor().fit(X_fit, y_fit)
A_fit, A_val, A_test = prep.transform(X_fit), prep.transform(X_val), prep.transform(X_test)
def fit_and_score(fit_cols, val_cols, test_cols):
model = new_lr().fit(fit_cols, y_fit)
return (roc_auc_score(y_val, model.predict_proba(val_cols)[:, 1]),
roc_auc_score(y_test, model.predict_proba(test_cols)[:, 1]))
start_val, start_test = fit_and_score(A_fit, A_val, A_test)
print('before we start: validation %.4f untouched test %.4f' % (start_val, start_test))
fit 4219 rows validation 1406 rows untouched test 1407 rows before we start: validation 0.8126 untouched test 0.8316
noise_rng = np.random.default_rng(0)
cur_fit, cur_val, cur_test = A_fit.copy(), A_val.copy(), A_test.copy()
best_val = start_val
history = [(0, 0, start_val, start_test)]
kept = 0
CANDIDATES = 200
for trial in range(1, CANDIDATES + 1):
col_fit = noise_rng.normal(size=(len(A_fit), 1))
col_val = noise_rng.normal(size=(len(A_val), 1))
col_test = noise_rng.normal(size=(len(A_test), 1))
val_score, test_score = fit_and_score(np.hstack([cur_fit, col_fit]),
np.hstack([cur_val, col_val]),
np.hstack([cur_test, col_test]))
if val_score > best_val: # the only rule: keep it if validation improves
best_val = val_score
kept += 1
cur_fit = np.hstack([cur_fit, col_fit])
cur_val = np.hstack([cur_val, col_val])
cur_test = np.hstack([cur_test, col_test])
history.append((trial, kept, val_score, test_score))
print('tried %d pure-noise features, kept %d of them' % (CANDIDATES, kept))
print()
print(' after keeping validation untouched test')
for trial, k, v, t in history[::6] + [history[-1]]:
print(' %2d features %.4f %.4f' % (k, v, t))
print()
print('validation improved %+.4f (%.4f -> %.4f)'
% (history[-1][2] - start_val, start_val, history[-1][2]))
print('the untouched test %+.4f (%.4f -> %.4f)'
% (history[-1][3] - start_test, start_test, history[-1][3]))
tried 200 pure-noise features, kept 43 of them
after keeping validation untouched test
0 features 0.8126 0.8316
6 features 0.8153 0.8310
12 features 0.8164 0.8304
18 features 0.8174 0.8301
24 features 0.8192 0.8295
30 features 0.8209 0.8305
36 features 0.8219 0.8301
42 features 0.8230 0.8296
43 features 0.8230 0.8298
validation improved +0.0104 (0.8126 -> 0.8230)
the untouched test -0.0019 (0.8316 -> 0.8298)200 pure-noise columns went in. 43 were kept, just because each one happened to nudge the validation score up when it was tried.
fig, ax = plt.subplots(figsize=(9, 4.2))
kept_axis = [h[1] for h in history]
ax.plot(kept_axis, [h[2] for h in history], 'o-', ms=4, color='#c44e52',
label='validation set, used to decide')
ax.plot(kept_axis, [h[3] for h in history], 's-', ms=4, color='#4c72b0',
label='test set, never looked at')
ax.set_xlabel('number of pure-noise features kept because validation went up')
ax.set_ylabel('ROC AUC')
ax.set_title('Every one of these features is random. Only one of these lines knows that.')
ax.legend(loc='center right', fontsize=9)
plt.tight_layout()
plt.show()
The red line rises by 0.0104, and every point on it is false. The blue line, measured on rows that played no part in any decision, moves by −0.0019, which is the right answer for adding 43 columns of noise.
Focus on that 0.0104. Section 4 lists what the real features in this series were worth, and 0.0104 is bigger than four of the five. The two sets of numbers come from different setups (one validation split here, 5-fold cross-validation on a smaller model there), so it's a rough comparison, but the point holds. Anyone running this loop by hand would think they'd found gains as big as most of the real work in notebooks 1 to 8, and they'd have nothing.
Notice what didn't happen. There was no leak in the usual sense. The preprocessor was fitted on the training rows only, and the validation rows never went into a fit call. Notebook 9's pipeline would accept every step of this loop, because the mistake isn't in any single step. A held-out set stops being truly held out as soon as you use it to choose between options, and it gets a little worse with every choice.
Section 4Over-engineering, measured on real data#
Section 3 added noise and kept only what the validation score seemed to reward. This section adds noise and keeps all of it, which is the more common version in practice: someone creates every ratio, product and polynomial term possible and gives the whole pile to the model
The belief being tested is that extra columns can't hurt, because the model can just give the useless ones a weight of zero. With unlimited data, that's close to true. With real, limited data it does hurt, because a model with more columns than it needs has more ways to fit its training rows exactly, and none of those extra ways work on new data.
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3, random_state=SEED, stratify=y)
prep_oe = make_preprocessor().fit(X_tr, y_tr)
base_tr, base_te = prep_oe.transform(X_tr), prep_oe.transform(X_te)
noise = np.random.default_rng(0)
print(' noise cols total train AUC test AUC gap')
curve = []
for n_noise in [0, 5, 20, 50, 100, 300, 1000, 3000]:
extra_tr = noise.normal(size=(len(base_tr), n_noise))
extra_te = noise.normal(size=(len(base_te), n_noise))
wide_tr = np.hstack([base_tr, extra_tr])
wide_te = np.hstack([base_te, extra_te])
model = LogisticRegression(max_iter=3000, random_state=SEED).fit(wide_tr, y_tr)
tr_auc = roc_auc_score(y_tr, model.predict_proba(wide_tr)[:, 1])
te_auc = roc_auc_score(y_te, model.predict_proba(wide_te)[:, 1])
curve.append((wide_tr.shape[1], tr_auc, te_auc))
print(' %10d %7d %.4f %.4f %+.4f' % (n_noise, wide_tr.shape[1], tr_auc, te_auc, tr_auc - te_auc))
noise cols total train AUC test AUC gap
0 18 0.8689 0.8220 +0.0469
5 23 0.8690 0.8218 +0.0472
20 38 0.8705 0.8204 +0.0501
50 68 0.8730 0.8149 +0.0581
100 118 0.8779 0.8124 +0.0655
300 318 0.8928 0.8026 +0.0902
1000 1018 0.9471 0.7388 +0.2083
3000 3018 1.0000 0.6537 +0.3463fig, ax = plt.subplots(figsize=(9, 4.2))
widths = [c[0] for c in curve]
ax.plot(widths, [c[1] for c in curve], 'o-', color='#c44e52', label='training rows')
ax.plot(widths, [c[2] for c in curve], 's-', color='#4c72b0', label='held-out rows')
ax.set_xscale('log')
ax.set_xlabel('total columns given to the model (18 real, the rest random)')
ax.set_ylabel('ROC AUC')
ax.set_title('Training AUC reaches %.4f at %d columns while held-out AUC falls to %.4f'
% (curve[-1][1], curve[-1][0], curve[-1][2]))
ax.legend(fontsize=9)
plt.tight_layout()
plt.show()
With 3,018 columns and 4,922 training rows, the model separates its training data perfectly, with an AUC of 1.0000, and scores 0.6537 on held-out rows, against 0.8220 with just the 18 real columns. It memorized its way to a worse model.
Notice the first row: the train-test gap is already +0.0469 with no noise added. A gap by itself doesn't mean anything is wrong. A model is fitted to do well on its training rows, so scoring a bit higher there is normal. The warning sign is the gap growing while the held-out score keeps falling.
The other side of over-engineering is quieter, with no dramatic number. It's the fifteenth feature that helps by 0.0005 and costs you a column to compute, check, document and keep working forever. Here's what the real features in this series were worth, added one at a time in the order the notebooks introduced them.
def logistic_cv(frame):
return cross_val_score(Pipeline([('scale', StandardScaler()), ('model', new_lr())]),
frame, y, cv=CV, scoring='roc_auc', n_jobs=1).mean()
chain = [('3 raw numeric columns', {})]
chain.append(('+ log_tenure and tenure_squared (nb1)',
dict(log_tenure=np.log1p(X['tenure']), tenure_squared=X['tenure'] ** 2)))
chain.append(('+ num_services (nb7)', dict(num_services=num_services)))
chain.append(('+ Contract as ordinal (nb4)',
dict(contract=X['Contract'].map({'Month-to-month': 0, 'One year': 1, 'Two year': 2}))))
chain.append(('+ PaperlessBilling', dict(paperless=(X['PaperlessBilling'] == 'Yes').astype(int))))
chain.append(('+ SeniorCitizen', dict(senior=X['SeniorCitizen'])))
frame = X[NUMERIC].copy()
previous = None
gains = []
for label, additions in chain:
frame = frame.assign(**additions)
score = logistic_cv(frame)
step = '' if previous is None else '%+.4f' % (score - previous)
if previous is not None:
gains.append((label, score - previous))
print(' %-40s %2d cols %.4f %s' % (label, frame.shape[1], score, step))
previous = score
print()
illusory = 0.0104
print('the noise loop in section 3 "gained" %+.4f on its validation set.' % illusory)
print('that is larger than %d of the %d real steps above.'
% (sum(1 for _, g in gains if g < illusory), len(gains)))
3 raw numeric columns 3 cols 0.8087 + log_tenure and tenure_squared (nb1) 5 cols 0.8159 +0.0072 + num_services (nb7) 6 cols 0.8241 +0.0082 + Contract as ordinal (nb4) 7 cols 0.8367 +0.0126 + PaperlessBilling 8 cols 0.8394 +0.0028 + SeniorCitizen 9 cols 0.8409 +0.0015 the noise loop in section 3 "gained" +0.0104 on its validation set. that is larger than 4 of the 5 real steps above.
Two things in that table matter more than the final number. The gains aren't in the order anyone thought of them. Contract came fourth and is the biggest gain in the list, at +0.0126, bigger than the log and square transforms that opened the series. Diminishing returns is a real pattern over a project, but you can't use it to decide when to stop, because the next idea might be the best one yet.
Each step's gain also depends on what came before it. Contract alone added +0.0191 in notebook 4, but only +0.0126 here, because num_services had already picked up part of its information (notebook 7 showed how tangled the two are). What a feature is worth is always relative to the features already in the model, as notebook 9 put it.
The last two features add +0.0028 and +0.0015. Both are positive under cross-validation, though a gain as small as 0.0015 could easily shrink on new data. Whether they're worth keeping is no longer a modeling question. It depends on whether PaperlessBilling is a stable field that will mean the same thing eighteen months from now, and that's exactly the kind of judgment a feature registry records, which section 7 covers.
Section 5The model was fine, the world moved#
Every failure so far existed on the day the model was trained. This one shows up later, and no amount of care during training prevents it
Notebook 5 already found this without naming it: a model validated with random KFold scored 0.6006, while the same model trained on 2011 and tested on 2012 scored 0.4096. The model here uses a few more calendar and weather columns than notebook 5's, so the numbers differ, but the direction and the reason are the same.
mean_2011 = bike.loc[bike['year'] == 2011, 'cnt'].mean()
mean_2012 = bike.loc[bike['year'] == 2012, 'cnt'].mean()
print('average hourly rentals, 2011 %.1f 2012 %.1f growth %.1f%%'
% (mean_2011, mean_2012, 100 * (mean_2012 / mean_2011 - 1)))
fixed = bike_model([]).fit(train_2011[BIKE_BASE], train_2011['cnt'])
print()
profile_2011 = bike[bike['year'] == 2011].groupby('hr')['cnt'].mean()
profile_2012 = bike[bike['year'] == 2012].groupby('hr')['cnt'].mean()
print('shape of the daily demand curve, 2011 against 2012: correlation %.4f'
% profile_2011.corr(profile_2012))
print()
print('random KFold over both years R2 %.4f'
% cross_val_score(bike_model([]), usable[BIKE_BASE], usable['cnt'],
cv=KFold(5, shuffle=True, random_state=SEED), scoring='r2', n_jobs=1).mean())
print('train on 2011, predict all of 2012 R2 %.4f'
% r2_score(test_2012['cnt'], fixed.predict(test_2012[BIKE_BASE])))
average hourly rentals, 2011 143.8 2012 234.7 growth 63.2% shape of the daily demand curve, 2011 against 2012: correlation 0.9992 random KFold over both years R2 0.6253 train on 2011, predict all of 2012 R2 0.4294
Average hourly rentals grew 63.2% from 2011 to 2012, and a random KFold estimate of 0.6253 sits well above the honest 2011-to-2012 estimate of 0.4294. The random split let the model learn from December to help predict June, which answers a question nobody needed answered.
monthly = []
for month in range(1, 13):
block = test_2012[test_2012['dteday'].dt.month == month]
predicted = fixed.predict(block[BIKE_BASE])
monthly.append((month, r2_score(block['cnt'], predicted),
(predicted - block['cnt']).mean(), block['cnt'].mean()))
print(' month R2 mean error actual mean')
for month, r2, bias, actual in monthly:
print(' 2012-%02d %+.4f %+8.1f %8.1f' % (month, r2, bias, actual))
print()
print('the model under-predicts in %d of 12 months' % sum(1 for _, _, b, _ in monthly if b < 0))
month_index = np.array([m[0] for m in monthly])
print('correlation between month and R2: %+.3f' % np.corrcoef(month_index, [m[1] for m in monthly])[0, 1])
print('correlation between month and mean error: %+.3f' % np.corrcoef(month_index, [m[2] for m in monthly])[0, 1])
month R2 mean error actual mean 2012-01 +0.3619 -62.4 130.6 2012-02 +0.3718 -71.9 149.0 2012-03 +0.3396 -103.3 221.9 2012-04 +0.4302 -88.6 242.7 2012-05 +0.4306 -90.6 263.3 2012-06 +0.4127 -89.0 281.7 2012-07 +0.4401 -73.9 273.7 2012-08 +0.3891 -98.0 288.3 2012-09 +0.2664 -125.9 303.6 2012-10 +0.3234 -110.1 280.8 2012-11 +0.4397 -68.3 212.6 2012-12 +0.4940 -55.8 166.7 the model under-predicts in 12 of 12 months correlation between month and R2: +0.209 correlation between month and mean error: -0.082
fig, axes = plt.subplots(1, 2, figsize=(11, 3.8))
months = [m[0] for m in monthly]
axes[0].plot(months, [m[1] for m in monthly], 'o-', color='#4c72b0')
axes[0].set_ylim(0, 0.6); axes[0].set_xlabel('month of 2012'); axes[0].set_ylabel('R2')
axes[0].set_title('R2 wobbles and never trends')
axes[1].bar(months, [m[2] for m in monthly], color='#c44e52')
axes[1].axhline(0, color='black', lw=1)
axes[1].set_xlabel('month of 2012'); axes[1].set_ylabel('mean prediction error')
axes[1].set_title('The bias is below zero every single month')
plt.tight_layout()
plt.show()
This is why drift is hard to catch, and why the two charts are worth keeping separate. If you only watched R-squared, you'd see it wander between 0.2664 and 0.4940 over the year, with only a weak +0.209 correlation with the month, and you'd reasonably conclude the model was noisy but fine.
If you watched the average error instead, you'd see it below zero in all twelve months, between 56 and 126 rentals short every time, and you'd know within weeks that something systematic was wrong.
The level of demand shifted, but the pattern didn't. The shape of the daily demand curve in 2012 correlates 0.9992 with 2011, so rush hour is still rush hour. There's just 63% more of everything. The model has no input that tracks this overall growth (no year or trend feature), so it can't represent the shift, and it ends up wrong in one consistent direction instead of randomly.
This matters for monitoring. Watching the inputs alone wouldn't have caught it, because 2012's weather and calendar looked much like 2011's. What changed was how many rentals the same inputs produce. When the incoming data itself changes, that's usually called data drift, and input monitoring catches it. When the link between inputs and outcome changes, as here, it's called concept drift, and only the errors show it, which needs the real outcomes to arrive.
print('retraining on everything available before each month:')
for month in (1, 4, 7, 10):
cutoff = pd.Timestamp(2012, month, 1)
history_rows = usable[usable['dteday'] < cutoff]
block = usable[(usable['year'] == 2012) & (usable['dteday'].dt.month == month)]
refreshed = bike_model([]).fit(history_rows[BIKE_BASE], history_rows['cnt'])
print(' 2012-%02d %5d rows of history refreshed R2 %+.4f frozen 2011 model %+.4f'
% (month, len(history_rows),
r2_score(block['cnt'], refreshed.predict(block[BIKE_BASE])),
monthly[month - 1][1]))
retraining on everything available before each month: 2012-01 8644 rows of history refreshed R2 +0.3619 frozen 2011 model +0.3619 2012-04 10820 rows of history refreshed R2 +0.4731 frozen 2011 model +0.4302 2012-07 13002 rows of history refreshed R2 +0.5148 frozen 2011 model +0.4401 2012-10 15210 rows of history refreshed R2 +0.4779 frozen 2011 model +0.3234
Retraining recovers a good part of the gap, and the gap grows the longer the frozen model stays in place. In January the two are identical, because "everything before now" is still exactly the 2011 training set. By October it's 0.4779 against 0.3234: the frozen model's training data now ends more than nine months back, and it has fallen well behind.
Section 6The metric is a business decision, not a modelling one#
Notebook 6 found that the same change to a spam model was worth +0.0033 of AUC and +0.0872 of recall, and which number you quote decides whether the change looks worth it. Notebook 6 also warned that recall depends on the cutoff. This section shows how much
ROC AUC sums up a model's performance across every possible decision threshold. But nobody deploys every threshold. Deploying a model means picking one, and that choice moves the number your colleagues care about far more than any feature in this notebook.
churn_model = make_pipe().fit(X_tr, y_tr)
scores = churn_model.predict_proba(X_te)[:, 1]
print('ROC AUC %.4f, and it is the same number for every row below' % roc_auc_score(y_te, scores))
print()
print(' threshold flagged recall precision f1')
for threshold in (0.60, 0.50, 0.40, 0.30, 0.25, 0.20):
flagged = (scores >= threshold).astype(int)
print(' %.2f %5d %.4f %.4f %.4f'
% (threshold, flagged.sum(), recall_score(y_te, flagged),
precision_score(y_te, flagged), f1_score(y_te, flagged)))
print()
print('%d of the %d held-out customers actually churned' % (y_te.sum(), len(y_te)))
ROC AUC 0.8360, and it is the same number for every row below
threshold flagged recall precision f1
0.60 296 0.3601 0.6824 0.4714
0.50 459 0.5330 0.6514 0.5863
0.40 638 0.6684 0.5878 0.6255
0.30 836 0.7754 0.5203 0.6228
0.25 951 0.8235 0.4858 0.6111
0.20 1040 0.8503 0.4587 0.5959
561 of the 2110 held-out customers actually churnedOne model, one AUC, and recall ranging from 0.3601 to 0.8503 depending on a number someone has to pick. At the default threshold of 0.50, the model catches just over half of the customers who leave. At 0.25 it catches 82%, but flags 951 customers instead of 459.
Neither threshold is more correct in a way the model can decide. The right answer depends on what a retention offer costs and what a leaving customer is worth, and that's a conversation with the business, not a setting to tune quietly.
This belongs in a pitfalls notebook because 0.5 is a default nobody chose. You get it automatically when you call predict() instead of predict_proba(). It's rarely the right choice for an unbalanced problem, and a useful model can be judged a failure just because it inherited it without anyone checking.
If you do tune the threshold, tune it on validation data, or with scikit-learn's TunedThresholdClassifierCV, never on the final test set. Otherwise it's section 3's used-up validation set all over again.
Section 7Working by hypothesis, and writing down why#
Section 3 showed what happens when features come from a random number generator and are judged only on a validation score. The defense isn't a cleverer statistical test. It's having a real reason before you start
The loop that works is simple. State what you believe about the domain and why. Build the feature that would show it. Measure it on data that played no part in forming the belief. Then look at where the model is still wrong to find the next idea to test.
This series has followed that loop from the start. num_services came from noticing nine service columns sitting unused. Contract as an ordinal came from knowing that a two-year contract is a commitment, not a random label. Neither came from searching through candidates.
The practical test: can you finish the sentence "I added this because…"? If you can't, the feature is probably a validation-set coincidence pretending to be an idea.
That's also what a feature registry is for. Every number in the one below was measured earlier in this series, and the last column is the one people usually forget to write down.
registry = pd.DataFrame([
dict(feature='log_tenure, tenure_squared',
source='tenure', hypothesis='churn falls steeply then flattens, so the effect is not linear',
measured='+0.0072 over raw numeric', notebook='1, 3',
watch_out='two columns for one idea; a spline would do the same job'),
dict(feature='num_services',
source='9 service columns', hypothesis='customers with more of the product are more entangled with it',
measured='+0.0082 in the chain', notebook='7',
watch_out='univariate AUC 0.4967; the relationship is an inverted U, and it is partly Contract in disguise'),
dict(feature='Contract (ordinal)',
source='Contract', hypothesis='month-to-month, one year, two year is a commitment ladder, not three labels',
measured='+0.0126, the largest single gain', notebook='4',
watch_out='the category order must be stated explicitly or it is arbitrary'),
dict(feature='bundle (target encoded)',
source='9 service columns joined', hypothesis='the combination of services says more than the count',
measured='+0.0163 alone, +0.0017 in the full model', notebook='4, 9',
watch_out='322 categories; must be cross-fitted or it leaks 0.0201'),
])
pd.set_option('display.max_colwidth', 52)
print(registry[['feature', 'hypothesis', 'measured']].to_string(index=False))
print()
print('the column that saves you later:')
for _, row in registry.iterrows():
print(' %-26s %s' % (row['feature'], row['watch_out']))
feature hypothesis measured
log_tenure, tenure_squared churn falls steeply then flattens, so the effect is not linear +0.0072 over raw numeric
num_services customers with more of the product are more entangled with it +0.0082 in the chain
Contract (ordinal) month-to-month, one year, two year is a commitment ladder, not three labels +0.0126, the largest single gain
bundle (target encoded) the combination of services says more than the count +0.0163 alone, +0.0017 in the full model
the column that saves you later:
log_tenure, tenure_squared two columns for one idea; a spline would do the same job
num_services univariate AUC 0.4967; the relationship is an inverted U, and it is partly Contract in disguise
Contract (ordinal) the category order must be stated explicitly or it is arbitrary
bundle (target encoded) 322 categories; must be cross-fitted or it leaks 0.0201The watch_out column pays for itself, usually much later than you'd expect. Six months from now, nobody will remember that bundle leaks two points of AUC if someone target-encodes it the obvious way, or that num_services looks useless on a single-column screen even though it helps the model. Someone will run a routine feature-selection pass, see 0.4967, drop the column, and never know what they lost.
Section 8What the series covered#
Ten notebooks, one dataset through most of them, and one rule that kept coming back in a different form until notebook 9 built it in
| Notebook | What it actually established |
|---|---|
| 1, Foundations | What a feature is, why the target isn't one, and the ID trap: customerID scores 100% on training rows and 61.3% on test, below the 73.5% from always guessing the majority class. |
| 2, Cleaning | Missing values that carry information, filling methods judged against the real values instead of a falling NaN count, and outliers that turned out to be a shared ticket, not a data error. |
| 3, Numeric transforms | Scaling changes units, never shape. Trees don't care either way. Raw polynomial features pushed R-squared to −11.16 before a log transform brought it back to 0.6505. |
| 4, Categorical encoding | An ordinal order is knowledge you have to supply. One-hot on 322 categories scored worse than dropping the column. Target encoding leaks about two points of AUC unless it's cross-fitted. |
| 5, Datetime | Cyclical encoding isn't automatically right. One-hot on the 24 hours beat sine and cosine, 0.6006 to 0.4236, because one wave can only make one peak and demand has two. |
| 6, Text | Nine hand-built style features beat a 3,802-word vocabulary on F1, and scikit-learn's own stopword list includes the word not. |
| 7, Construction | TotalCharges is tenure times MonthlyCharges with an R-squared of 0.999, so most math on those three columns just renames a column that already exists. |
| 8, Selection | Five methods, three different answers, and selection done outside the folds reporting 74.5% accuracy on a target that was a pure coin flip. |
| 9, Pipelines | The leakage rule built into the code, and each kind of leak measured side by side: 0.000001 AUC for a scaler, 0.27 of accuracy for a selector on pure noise. |
| 10, this one | The failures that get past a pipeline anyway, because a pipeline was never responsible for them. |
If there's one habit to take from all ten notebooks, it's this: measure what you're about to claim, on rows that played no part in producing the claim. Almost every mistake in this notebook, and many from the nine before it, comes from skipping that step.
Where to go next. The obvious gap in this series is time series: lags, rolling and expanding windows built correctly, and the validation methods forecasting needs. Section 2 is a small first taste of why that topic needs its own course.
Beyond that, gradient boosting handles missing values and raw categories natively, which changes which of these transformations are still worth doing. SHAP values show what a trained model is using for each individual prediction, not just on average. And text embeddings from pretrained language models are the next step beyond word counts. They don't need a huge dataset of your own, because the heavy training was done elsewhere.
Wrap-upKey things to remember from this notebook#
A correct pipeline doesn't make your project correct. Notebook 9 made preprocessing leakage impossible, and every failure in this notebook passed straight through that same pipeline.
The most expensive column is the one recorded after the answer. Adding one made-up cancellation_ticket flag took the honest 0.8416 to 0.9809. That's not a modeling result. It's a data-plumbing question, and the only reliable way to catch it is to ask, column by column, when each field actually gets filled in.
Your checks can point the wrong way, not just fail to help. The leaky model's folds agreed more closely than the honest model's, 0.0028 against 0.0059, so a fold-stability check would have called the broken model more trustworthy.
Single-column screening cuts both ways, and neither edge is sharp. The leaked column scored 0.9406 on its own, a real warning sign. But num_services scores 0.4967, a coin flip, and still adds +0.0082 to the model, because a single-column AUC can't see a relationship that goes up and then down.
Shift before you roll. A centred 3-hour window correlates 0.977 with the target and lifts R-squared from 0.6774 to 0.9716, which anyone would catch. A centred 24-hour window correlates 0.436 against a legitimate 0.401, which nobody would catch, and still adds 0.0594 of R-squared that won't exist in production. The leak that costs real money is the one small enough to get past review. And even a trailing window only works if its inputs exist when you predict.
A validation set gets used up as you use it. Keeping 43 pure-noise features because each nudged the validation score up raised that score by 0.0104, while the untouched test set moved −0.0019. That imaginary 0.0104 is bigger than four of the five real features this series added to the Telco model.
Train-test gaps are normal. Growing ones aren't. With 18 real columns the gap was already +0.0469. With 3,000 columns of added noise, training AUC hit exactly 1.0000 while held-out AUC fell to 0.6537.
Returns shrink over a project, but not in the order you find things. Contract as an ordinal column came fourth and was the biggest single gain, +0.0126, bigger than the transforms the series started with. And a feature's gain depends on what's already in the model: Contract was worth +0.0191 on its own in notebook 4.
Watch errors, not just scores. A 2011 model run through 2012 had an R-squared that wandered between 0.2664 and 0.4940, with only a +0.209 correlation with the month. Its average error, meanwhile, was negative in all twelve months. Ridership had grown 63%, and only the errors showed it. The inputs looked unchanged, so input monitoring alone would have missed this kind of drift.
The threshold is a business decision, and nobody actually chose 0.5. The same model, with the same AUC, catches 36% of churners at a threshold of 0.60 and 85% at 0.20. Choose the threshold with the business, and tune it on validation data, never on the final test set.
Write down why each feature exists and what to watch out for. In six months, nobody will remember that a column with a single-column AUC of 0.4967 was one of the model's more useful features.
In one sentence: measure what you're about to claim, on rows that played no part in the claim, and be able to say plainly why every column in the model is there.
That's the end of the series. Ten notebooks, and the common thread was never really a technique. It was the habit of checking, on data that couldn't have been used to produce the answer, whether what you believe is actually true.