Feature Engineering · Notebook 9 companion articlePutting It All Together in Pipelines

Start hereIntroduction#

Notebook 3 stated the rule in full, and it has come back in every notebook since: anything with a fit method learns from data, and everything it learns must come from the training set. Followed correctly by hand, eight times, the rule works. But it depends on nobody ever forgetting it, in the ninth notebook or the ninetieth. Pipeline and ColumnTransformer take the rule out of your memory by making the wrong order of steps hard to write by accident.

That's the honest reason to use them. The less honest reason, repeated so often it's become folklore, is that skipping a pipeline will wreck your score. Section 2 measures what these leaks are actually worth, for four of the six transformer types. The answer ranges from a millionth of an AUC point to a huge 0.27 of accuracy. The whole case for a pipeline is that you never know in advance which end of that range you're at.

Section 7 then picks up where notebook 8 ended: if you try many options and keep the best score, how much does that choice inflate the number you report?

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import joblib, os, sklearn

from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import (StandardScaler, MinMaxScaler, RobustScaler,
                                   OneHotEncoder, OrdinalEncoder, TargetEncoder,
                                   FunctionTransformer)
from sklearn.impute import SimpleImputer
from sklearn.linear_model import LogisticRegression
from sklearn.feature_selection import SelectKBest, f_classif
from sklearn.model_selection import (cross_val_score, StratifiedKFold, GridSearchCV,
                                     train_test_split)
from sklearn.metrics import roc_auc_score

SEED = 42
CV = StratifiedKFold(n_splits=5, shuffle=True, random_state=SEED)

def new_lr():
    return LogisticRegression(max_iter=2000, random_state=SEED)

print('sklearn', sklearn.__version__)

sklearn 1.8.0

telco = pd.read_csv('datasets/telco_churn.csv')
telco['TotalCharges'] = pd.to_numeric(telco['TotalCharges'], errors='coerce')

print('missing values before the tenure filter:', int(telco.isna().sum().sum()))
telco = telco[telco['tenure'] > 0].reset_index(drop=True)
print('missing values after the tenure filter: ', int(telco.isna().sum().sum()))

SERVICES = ['PhoneService', 'MultipleLines', 'InternetService', 'OnlineSecurity',
            'OnlineBackup', 'DeviceProtection', 'TechSupport', 'StreamingTV',
            'StreamingMovies']

# nb4 built this: the nine service columns pasted into one high-cardinality label
telco['bundle'] = telco[SERVICES].apply(lambda row: '|'.join(row.astype(str)), axis=1)

y = (telco['Churn'] == 'Yes').astype(int)
X = telco.drop(columns=['Churn'])

print()
print('rows', len(X), ' churn rate %.4f' % y.mean())
print('bundle categories:', X['bundle'].nunique())

missing values before the tenure filter: 11
missing values after the tenure filter:  0

rows 7032  churn rate 0.2658
bundle categories: 322

These are the same 7,032 rows and 26.6% churn rate used since the first notebook, so every number stays comparable. Two details in that output come back later. The eleven missing values are all in rows where tenure == 0, so after filtering those rows there are no missing values left. An imputer still goes into the pipeline anyway, for a reason section 8 makes clear.

And bundle, notebook 4's column that joins nine service columns into one 322-category label, is the real high-cardinality decision here. customerID has even more values, 7,032, but it's just an ID to throw away.

Section 1A pipeline is one object with one fit#

The whole safety idea behind this notebook fits in one sentence. Here's exactly what that sentence rules out

A Pipeline is a list of (name, step) pairs. Every step except the last must have both fit and transform. The last one only needs fit, which is why it's usually the model. Calling .fit() goes through the list once, calling fit_transform on each step and passing the result to the next. Calling .predict() goes through the same list calling only transform, and never calls fit again.

That last point is the whole safety idea: there's no way for predict to refit anything.

NUMERIC = ['tenure', 'MonthlyCharges', 'TotalCharges']

simple = Pipeline([
    ('impute', SimpleImputer(strategy='median')),
    ('scale',  StandardScaler()),
    ('model',  new_lr()),
])

simple.fit(X[NUMERIC], y)

print('steps:', list(simple.named_steps))
print()
print('what each fitted step learned, all of it from the rows passed to fit:')
print('  imputer medians:', simple.named_steps['impute'].statistics_.round(2))
print('  scaler means:   ', simple.named_steps['scale'].mean_.round(2))
print('  model coefs:    ', simple.named_steps['model'].coef_.round(3))
print()
print('5-fold AUC %.4f  (notebook 1 baseline: 0.8087)'
      % cross_val_score(simple, X[NUMERIC], y, cv=CV, scoring='roc_auc', n_jobs=1).mean())

steps: ['impute', 'scale', 'model']

what each fitted step learned, all of it from the rows passed to fit:
  imputer medians: [  29.     70.35 1397.48]
  scaler means:    [  32.42   64.8  2283.3 ]
  model coefs:     [[-1.616  0.913  0.299]]
5-fold AUC 0.8087  (notebook 1 baseline: 0.8087)

0.8087, the same baseline as notebook 1 and every notebook since. The pipeline didn't change the model at all. It changed where the medians and means are allowed to come from. And the fitted pieces are still easy to reach through named_steps. A pipeline isn't a black box that trades visibility for safety. It's the same objects, with the plumbing handled for you.

Section 2How much is each kind of leak actually worth#

The leakage rule has come up for every type of transformer in this series. Here, four of the six are measured side by side, with the same model and the same seed

Every row in that table is a fit, and every one of them can be run on the whole dataset before splitting. The usual advice stops at "don't do that." This section measures it instead for the scaler, the imputer, the target encoder and the selector (the vectorizer and group-aggregate leaks were measured in notebooks 6 and 7). The sizes turn out to be very different, and that difference is the useful part.

The recipe is the same for each row. The leaky version fits the transformer on all the data, then cross-validates the model on the result. The honest version puts the same transformer inside a pipeline, so each fold refits it on its own training rows. Same folds, same model, same seed every time.

leaks = {}

# --- scaler, on Telco ---
leaky_X = StandardScaler().fit_transform(X[NUMERIC])          # fitted on all 7,032 rows
a = cross_val_score(new_lr(), leaky_X, y, cv=CV, scoring='roc_auc', n_jobs=1).mean()
b = cross_val_score(Pipeline([('scale', StandardScaler()), ('model', new_lr())]),
                    X[NUMERIC], y, cv=CV, scoring='roc_auc', n_jobs=1).mean()
leaks['scaler (Telco)'] = (a, b)
print('scaler   leaky %.6f   honest %.6f   difference %+.6f' % (a, b, a - b))
first_fold = next(iter(CV.split(X, y)))[0]
all_rows = StandardScaler().fit(X[NUMERIC]).mean_
one_fold = StandardScaler().fit(X[NUMERIC].iloc[first_fold]).mean_
print('  column means, all 7,032 rows vs the %d rows of one training fold:' % len(first_fold))
print('   ', all_rows.round(3))
print('   ', one_fold.round(3))
print('  they differ by %s, at most %.2f%%'
      % (np.round(np.abs(one_fold - all_rows), 3), (np.abs(one_fold - all_rows) / all_rows).max() * 100))

scaler   leaky 0.808677   honest 0.808676   difference +0.000001
  column means, all 7,032 rows vs the 5625 rows of one training fold:
    [  32.422   64.798 2283.3  ]
    [  32.522   64.878 2297.903]
  they differ by [ 0.1    0.079 14.602], at most 0.64%

One part in a million. That's a real effect, but a tiny one: the slightly different means and standard deviations shift the regularized model's fit by a hair.

The reason it's so small is right below. A mean over all 7,032 rows and a mean over one training fold's 5,625 rows differ by less than one percent on every column, because a mean is stable, and removing a random fifth of the rows barely moves it. So the leak tutorials use to scare people is, here at least, worth nothing.

# --- imputer, on Titanic, because Telco has nothing left to impute ---
titanic = pd.read_csv('datasets/titanic.csv')
X_tt = titanic[['Age', 'Fare', 'Pclass', 'SibSp', 'Parch']]
y_tt = titanic['Survived']
print('Titanic Age missing: %d of %d rows' % (X_tt['Age'].isna().sum(), len(X_tt)))

leaky_X = SimpleImputer(strategy='median').fit_transform(X_tt)   # median over all 891 rows
a = cross_val_score(Pipeline([('scale', StandardScaler()), ('model', new_lr())]),
                    leaky_X, y_tt, cv=CV, scoring='roc_auc', n_jobs=1).mean()
b = cross_val_score(Pipeline([('impute', SimpleImputer(strategy='median')),
                              ('scale', StandardScaler()), ('model', new_lr())]),
                    X_tt, y_tt, cv=CV, scoring='roc_auc', n_jobs=1).mean()
leaks['imputer (Titanic)'] = (a, b)
print('imputer  leaky %.6f   honest %.6f   difference %+.6f' % (a, b, a - b))
print('  Age median over everything: %.2f' % X_tt['Age'].median())
print('  Age median per training fold:', [X_tt.iloc[tr]['Age'].median() for tr, _ in CV.split(X_tt, y_tt)])

Titanic Age missing: 177 of 891 rows
imputer  leaky 0.725968   honest 0.725834   difference +0.000134
  Age median over everything: 28.00
  Age median per training fold: [np.float64(28.0), np.float64(29.0), np.float64(28.0), np.float64(28.0), np.float64(28.0)]

A little bigger, but still nothing you could act on. The medians for each fold are 28, 29, 28, 28 and 28, against 28.00 for the whole file, so four of the five folds fill in exactly the value the leaky version used. Both of these cases stay small for the same reason: both transformers learn a couple of summary numbers per column, and removing 20% of the rows barely changes a summary. Next, the transformer learns one number per category instead.

# --- target encoding, on the 322-category bundle ---
category_means = y.groupby(X['bundle']).mean()          # computed from every row, labels included
leaky_frame = X[NUMERIC].assign(bundle_encoded=X['bundle'].map(category_means))
a = cross_val_score(Pipeline([('scale', StandardScaler()), ('model', new_lr())]),
                    leaky_frame, y, cv=CV, scoring='roc_auc', n_jobs=1).mean()

honest = Pipeline([
    ('prep', ColumnTransformer([('num', StandardScaler(), NUMERIC),
                                ('cat', TargetEncoder(random_state=SEED), ['bundle'])])),
    ('model', new_lr()),
])
b = cross_val_score(honest, X[NUMERIC + ['bundle']], y, cv=CV, scoring='roc_auc', n_jobs=1).mean()
leaks['target encoder (bundle)'] = (a, b)
print('encoder  leaky %.6f   honest %.6f   difference %+.6f' % (a, b, a - b))

sizes = X['bundle'].value_counts()
print()
print('  bundle categories: %d, median size %d rows, %d of them have 5 rows or fewer'
      % (len(sizes), int(sizes.median()), int((sizes <= 5).sum())))

encoder  leaky 0.845031   honest 0.824909   difference +0.020121

  bundle categories: 322, median size 9 rows, 84 of them have 5 rows or fewer

Two points of AUC. Notebook 4 measured this same encoder on this same column and got 0.8450 leaked against 0.8250 honest. The leaked figure matches exactly here, and the honest one is a ten-thousandth away. The reason is category size. The median bundle has only 9 rows, and 84 bundles have 5 or fewer. In a category that small, each row's own label is a big part of the category's average, so the encoder is handing the model a smoothed copy of the answer.

Notebook 7 showed that a row's own share of its group's average is about 1 / group_size, so the leak shrinks as groups grow. That's why a 322-category column leaks and a 3-category column hardly does. TargetEncoder avoids this by cross-fitting internally. That's not something a pipeline gives you automatically. It's something this class does on purpose, and it's worth knowing the difference between the two kinds of protection.

One caveat on the number: the leaky version uses a plain, unsmoothed category mean, while TargetEncoder also smooths small categories toward the overall rate. So the 0.020 measures the naive approach against the correct one, and a small part of it may come from smoothing rather than the leak alone.

# --- selection, twice: once where it is harmless, once where it is not ---
rng = np.random.default_rng(0)
noise_X = pd.DataFrame(rng.normal(size=(200, 2000)))
coin_y = pd.Series(rng.integers(0, 2, 200))

chosen = SelectKBest(f_classif, k=20).fit(noise_X, coin_y).get_support()
a = cross_val_score(new_lr(), noise_X.loc[:, chosen], coin_y, cv=CV, scoring='accuracy', n_jobs=1).mean()
b = cross_val_score(Pipeline([('select', SelectKBest(f_classif, k=20)), ('model', new_lr())]),
                    noise_X, coin_y, cv=CV, scoring='accuracy', n_jobs=1).mean()
leaks['selector (2000 noise columns)'] = (a, b)
print('selector on pure noise, accuracy: leaky %.4f  honest %.4f  difference %+.4f' % (a, b, a - b))

real = X[NUMERIC].assign(log_tenure=np.log1p(X['tenure']),
                         tenure_squared=X['tenure'] ** 2,
                         num_services=X[SERVICES].isin(['Yes', 'DSL', 'Fiber optic']).sum(axis=1))
chosen = SelectKBest(f_classif, k=3).fit(real, y).get_support()
a = cross_val_score(Pipeline([('scale', StandardScaler()), ('model', new_lr())]),
                    real.loc[:, chosen], y, cv=CV, scoring='roc_auc', n_jobs=1).mean()
b = cross_val_score(Pipeline([('select', SelectKBest(f_classif, k=3)),
                              ('scale', StandardScaler()), ('model', new_lr())]),
                    real, y, cv=CV, scoring='roc_auc', n_jobs=1).mean()
leaks['selector (real Telco columns)'] = (a, b)
print('selector on real columns,  AUC:   leaky %.4f  honest %.4f  difference %+.4f' % (a, b, a - b))

# why the two answers differ: do the folds agree on what to select?
noise_picks = [set(np.flatnonzero(SelectKBest(f_classif, k=20)
                                  .fit(noise_X.iloc[tr], coin_y.iloc[tr]).get_support()))
               for tr, _ in CV.split(noise_X, coin_y)]
real_picks = [tuple(real.columns[SelectKBest(f_classif, k=3)
                                 .fit(real.iloc[tr], y.iloc[tr]).get_support()])
              for tr, _ in CV.split(real, y)]
print()
print('noise: %d columns chosen by all 5 folds, %d chosen by at least one'
      % (len(set.intersection(*noise_picks)), len(set().union(*noise_picks))))
print('real:  %d distinct selections across the 5 folds, %s'
      % (len(set(real_picks)), list(real_picks[0])))

selector on pure noise, accuracy: leaky 0.7450  honest 0.4750  difference +0.2700
selector on real columns,  AUC:   leaky 0.7412  honest 0.7412  difference +0.0000

noise: 1 columns chosen by all 5 folds, 68 chosen by at least one
real:  1 distinct selections across the 5 folds, ['tenure', 'log_tenure', 'tenure_squared']

Same transformer, same mistake, two very different answers. On 2,000 pure-noise columns with a coin-flip target, selecting first reports 74.5% accuracy for a problem whose true answer is 50%. On six real Telco columns, the same mistake costs nothing. The last two lines of output show why, and it's the most useful idea in this section: the less stable the selection, the bigger the leak.

On the real columns, all five folds pick the same three, so choosing them in advance changes nothing. On the noise columns, one column survives every fold and 68 are picked by at least one. Each of those disagreements is the selector fitting itself to rows it's about to be tested on. Notebook 8 reached the same conclusion from the other direction, with the same two counts.

summary = (pd.DataFrame([(k, v[0], v[1], v[0] - v[1]) for k, v in leaks.items()],
                        columns=['transformer', 'leaky', 'honest', 'difference'])
             .sort_values('difference'))
print(summary.to_string(index=False, float_format=lambda v: '%.6f' % v))

fig, ax = plt.subplots(figsize=(9, 3.6))
labels = summary['transformer'].tolist()
vals = summary['difference'].clip(lower=1e-7).tolist()
colors = ['#4c72b0' if d < 0.01 else '#c44e52' for d in summary['difference']]
ax.barh(labels, vals, color=colors)
ax.set_xscale('log')
ax.set_xlim(1e-7, 1)
ax.set_xlabel('score the leak buys you (log scale)')
ax.set_title('The same mistake is worth five orders of magnitude more in some places than others')
for i, d in enumerate(summary['difference']):
    ax.text(max(d, 1e-7) * 1.4, i, ('0' if d <= 0 else '%.6f' % d).rstrip('0').rstrip('.'),
            va='center', fontsize=9)
plt.tight_layout()
plt.show()

That spread is the real argument, and it's different from the one usually given. If leaking a scaler cost two points of AUC, everyone would already be careful with scalers and nobody would need a pipeline to enforce it. Instead it costs a millionth of a point, so being careful with scalers feels like superstition. People then carry that carelessness over to encoders and selectors, which is exactly where it costs real money.

A pipeline is worth it because you can't tell which row of that table you're in until you've already done the work properly. And it costs nothing to be right about all six at once.

Notebook 6 made the same point: a leak's size depends on the setup, and the only way to measure it is to run the honest version too. The table above only exists because the honest version was calculated too.

Section 3ColumnTransformer, and seeing what came out of it#

A plain Pipeline applies every step to every column. That's fine when all columns are numbers, and useless as soon as they aren't

You can't scale PaymentMethod, and you can't one-hot encode MonthlyCharges. ColumnTransformer takes a list of (name, transformer, columns) entries, runs each one on its own group of columns, and joins the results side by side. Any entry can itself be a small Pipeline, which is how one group of columns gets more than one step.

NOMINAL   = ['gender', 'SeniorCitizen', 'Partner', 'Dependents',
             'PaperlessBilling', 'PaymentMethod', 'InternetService']
ORDINAL   = ['Contract']
CONTRACT_ORDER = [['Month-to-month', 'One year', 'Two year']]

for name, cols in [('numeric', NUMERIC), ('nominal', NOMINAL), ('ordinal', ORDINAL),
                   ('high cardinality', ['bundle']), ('constructed from', SERVICES)]:
    print('%-18s %d columns  %s' % (name, len(cols), cols if len(cols) < 8 else cols[:3] + ['...']))
print()
print('dropped entirely:  customerID (notebook 1: a perfect identifier, 7,032 distinct values)')

numeric            3 columns  ['tenure', 'MonthlyCharges', 'TotalCharges']
nominal            7 columns  ['gender', 'SeniorCitizen', 'Partner', 'Dependents', 'PaperlessBilling', 'PaymentMethod', 'InternetService']
ordinal            1 columns  ['Contract']
high cardinality   1 columns  ['bundle']
constructed from   9 columns  ['PhoneService', 'MultipleLines', 'InternetService', '...']

dropped entirely:  customerID (notebook 1: a perfect identifier, 7,032 distinct values)

def count_services(frame):
    """nb7's num_services: how many of the nine services this customer actually has."""
    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():
    return ColumnTransformer(
        transformers=[
            ('num', Pipeline([('impute', SimpleImputer(strategy='median')),
                              ('scale', StandardScaler())]), NUMERIC),
            ('svc', Pipeline([('count', service_counter),
                              ('scale', StandardScaler())]), SERVICES),
            ('ord', OrdinalEncoder(categories=CONTRACT_ORDER), ORDINAL),
            ('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_full_pipe():
    return Pipeline([('prep', make_preprocessor()), ('model', new_lr())])

make_full_pipe()

Pipeline(steps=[('prep',
                 ColumnTransformer(transformers=[('num',
                                                  Pipeline(steps=[('impute',
                                                                   SimpleImputer(strategy='median')),
                                                                  ('scale',
                                                                   StandardScaler())]),
                                                  ['tenure', 'MonthlyCharges',
                                                   'TotalCharges']),
                                                 ('svc',
                                                  Pipeline(steps=[('count',
                                                                   FunctionTransformer(feature_names_out=<function service_feature_names at 0x118f79120>,
                                                                                       func=<function count_services at 0x118f788...
                                                                              'year']]),
                                                  ['Contract']),
                                                 ('nom',
                                                  OneHotEncoder(drop='if_binary',
                                                                handle_unknown='ignore',
                                                                sparse_output=False),
                                                  ['gender', 'SeniorCitizen',
                                                   'Partner', 'Dependents',
                                                   'PaperlessBilling',
                                                   'PaymentMethod',
                                                   'InternetService']),
                                                 ('hi',
                                                  TargetEncoder(random_state=42),
                                                  ['bundle'])],
                                   verbose_feature_names_out=False)),
                ('model', LogisticRegression(max_iter=2000, random_state=42))])

That last line prints the pipeline itself, not a score. In a real notebook, scikit-learn shows it as a clickable diagram, which is the quickest way to confirm the nested structure is the one you meant to build. Next, the score, compared with the steps this series has built since notebook 1.

scores = {
    'numeric only (nb1)': X[NUMERIC],
    '+ log_tenure, tenure_squared (nb1)': X[NUMERIC].assign(
        log_tenure=np.log1p(X['tenure']), tenure_squared=X['tenure'] ** 2),
    '+ num_services (nb7)': X[NUMERIC].assign(
        log_tenure=np.log1p(X['tenure']), tenure_squared=X['tenure'] ** 2,
        num_services=X[SERVICES].isin(['Yes', 'DSL', 'Fiber optic']).sum(axis=1)),
}
for label, frame in scores.items():
    print('%-38s %.4f' % (label, cross_val_score(
        Pipeline([('scale', StandardScaler()), ('model', new_lr())]),
        frame, y, cv=CV, scoring='roc_auc', n_jobs=1).mean()))

full_cv = cross_val_score(make_full_pipe(), X, y, cv=CV, scoring='roc_auc', n_jobs=1).mean()
print('%-38s %.4f' % ('everything, one ColumnTransformer', full_cv))

numeric only (nb1)                     0.8087
+ log_tenure, tenure_squared (nb1)     0.8159
+ num_services (nb7)                   0.8241
everything, one ColumnTransformer      0.8416

0.8087, 0.8159 and 0.8241 match the series exactly, and the full pipeline reaches 0.8416. That's a real gain of +0.0329 over the numeric baseline, and all of it comes from columns earlier notebooks already worked out how to handle, not from anything clever added here. Contract in its proper order came from notebooks 1 and 4. num_services came from notebook 7. The one-hot settings came from notebook 4.

Section 7 checks this number on rows that played no part in producing it, because a five-fold average isn't yet something you can quote with confidence.

prep = make_preprocessor().set_output(transform='pandas')
transformed = prep.fit_transform(X, y)

print('%d raw columns in, %d features out' % (X.shape[1], transformed.shape[1]))
print()
print(transformed.head(3).round(3).to_string())

21 raw columns in, 18 features out

   tenure  MonthlyCharges  TotalCharges  num_services  Contract  gender_Male  SeniorCitizen_1  Partner_Yes  Dependents_Yes  PaperlessBilling_Yes  PaymentMethod_Bank transfer (automatic)  PaymentMethod_Credit card (automatic)  PaymentMethod_Electronic check  PaymentMethod_Mailed check  InternetService_DSL  InternetService_Fiber optic  InternetService_No  bundle
0  -1.280          -1.162        -0.994        -0.929       0.0          0.0              0.0          1.0             0.0                   1.0                                      0.0                                    0.0                             1.0                         0.0                  1.0                          0.0                 0.0   0.385
1   0.064          -0.261        -0.174        -0.064       1.0          1.0              0.0          0.0             0.0                   0.0                                      0.0                                    0.0                             0.0                         1.0                  1.0                          0.0                 0.0   0.138
2  -1.240          -0.364        -0.960        -0.064       0.0          1.0              0.0          0.0             0.0                   1.0                                      0.0                                    0.0                             0.0                         1.0                  1.0                          0.0                 0.0   0.227

set_output(transform='pandas') may be the most useful line in this notebook for everyday debugging. Without it, a ColumnTransformer returns an unlabeled numpy array, and finding the column that broke means counting on your fingers. With it, you get a real DataFrame with column names from get_feature_names_out.

Two details in that output are worth a closer look. verbose_feature_names_out=False keeps names as plain tenure instead of the default num__tenure. The longer default guarantees unique names, so only turn it off when you know your column groups can't clash.

And drop='if_binary' is why gender became one column instead of two. For a two-category column, the second column is just the first one flipped, so it adds no information and would only give a linear model a redundant coefficient.

One small inconsistency is worth admitting. Two columns reach the model unscaled: Contract as 0, 1, 2, and the target-encoded bundle, a churn rate between 0 and 1. Notebook 3 said a penalized model should see columns on a common scale. With the default penalty and several thousand rows the effect here is small, but a stricter pipeline would scale those two as well.

Here's one trap to meet on purpose, because the error message doesn't make the cause obvious.

sparse_prep = ColumnTransformer(
    [('nom', OneHotEncoder(handle_unknown='ignore'), NOMINAL)],   # sparse_output defaults to True
    remainder='drop').set_output(transform='pandas')
try:
    sparse_prep.fit_transform(X)
except Exception as err:
    print('%s: %s' % (type(err).__name__, err))

ValueError: Pandas output does not support sparse data. Set sparse_output=False to output pandas dataframes or disable Pandas output via` ohe.set_output(transform="default").

OneHotEncoder returns a sparse matrix by default, and scikit-learn's pandas output doesn't support sparse data, so the two don't work together. For a 322-category column, sparse is the right choice, and it's worth keeping even if you lose the DataFrame view. For seven columns with few categories it makes no real difference, and sparse_output=False is worth it there.

The other setting to choose on purpose is remainder. It defaults to 'drop', which throws away any column you didn't list. That sounds worrying, but it's actually the safest default.

passthrough = ColumnTransformer([('num', StandardScaler(), NUMERIC)],
                                remainder='passthrough', verbose_feature_names_out=False)
try:
    Pipeline([('prep', passthrough), ('model', new_lr())]).fit(X, y)
except Exception as err:
    print('remainder="passthrough": %s: %s' % (type(err).__name__, str(err).split(chr(10))[0]))
print()
print('remainder="drop" instead: customerID never reaches the model, and we said so on purpose.')

remainder="passthrough": ValueError: could not convert string to float: '7590-VHVEG'

remainder="drop" instead: customerID never reaches the model, and we said so on purpose.

With passthrough, customerID goes straight to the model, and logistic regression can't read a string like 7590-VHVEG. Here that's a loud, immediate failure, which is the good outcome. Notebook 1 showed the bad one: give the same column to a decision tree and it trains happily, scoring 100% on training rows and 61.3% on test, worse than the 73.5% from always predicting the majority class.

Listing your columns explicitly takes a few extra lines, and what you get is a schema: a new column that appears in next month's file won't silently sneak into your model.

Section 4Custom transformers, stateless and stateful#

Eventually you'll need a transformation scikit-learn doesn't have. The interface is small, and one question decides how to write it

The interface is small: a fit that returns self, and a transform that returns the transformed data. BaseEstimator and TransformerMixin add get_params, set_params and fit_transform for free. But before writing a class, ask one question: does this transformation actually learn anything?

num_services counts how many of nine columns say yes for a row. That count depends only on that one row: no training data, nothing stored, nothing that could leak. A transformation like that doesn't need a class, which is why it already appeared as FunctionTransformer(count_services, feature_names_out=service_feature_names).

A stateless transformer is the one case where doing the work outside the pipeline is actually safe, since there's no fit step to run on the wrong rows. It still belongs in the pipeline, though, so the saved model knows how to build that column when serving, without anyone having to remember the formula.

Now a transformer that does learn something. Notebook 4 found that one-hot encoding all 322 bundles scored 0.8007, worse than dropping the column, and that grouping rare categories together fixed it. But to group rare categories, you have to know which ones are rare, and that's a fact about the training rows, not about the column in general.

class RareCategoryGrouper(BaseEstimator, TransformerMixin):
    """Keep categories seen at least `min_count` times in training, send the rest to 'Other'."""

    def __init__(self, min_count=20, other_label='Other'):
        self.min_count = min_count
        self.other_label = other_label

    def fit(self, X, y=None):
        counts = pd.DataFrame(X).iloc[:, 0].value_counts()
        self.frequent_ = set(counts[counts >= self.min_count].index)
        self.n_seen_ = len(counts)
        return self

    def transform(self, X):
        column = pd.DataFrame(X).iloc[:, 0]
        return column.where(column.isin(self.frequent_), self.other_label).to_frame('bundle_grouped')

    def get_feature_names_out(self, input_features=None):
        return np.array(['bundle_grouped'])


on_everything = RareCategoryGrouper(min_count=20).fit(X[['bundle']])
X_dev, X_hold, y_dev, y_hold = train_test_split(X, y, test_size=0.2, random_state=SEED, stratify=y)
on_train_only = RareCategoryGrouper(min_count=20).fit(X_dev[['bundle']])

print('fitted on all 7,032 rows: %d categories kept, plus Other' % len(on_everything.frequent_))
print('fitted on the 5,625 training rows: %d categories kept, plus Other' % len(on_train_only.frequent_))
print()
moved = on_train_only.transform(X_hold[['bundle']])['bundle_grouped']
print('held-out rows landing in Other: %d of %d' % ((moved == 'Other').sum(), len(moved)))

fitted on all 7,032 rows: 76 categories kept, plus Other
fitted on the 5,625 training rows: 61 categories kept, plus Other

held-out rows landing in Other: 461 of 1407

76 categories against 61. Fitting the grouper on everything doesn't just nudge a few numbers. It changes the width of the feature matrix, using rows the model will be tested on. This case deserves a close look, because people miss it most often. RareCategoryGrouper never touches y, so it can't leak the target, and in that narrow sense it's safe.

But it can still make your list of features depend on the test rows, and that's what happened here: fifteen extra columns exist only because those rows were visible during fitting. Notebook 8 found the same issue with VarianceThreshold, which also ignores the target but still decides which columns exist.

Part of this gap is simply that min_count=20 is an absolute count: with 25% more rows, more categories reach 20. A threshold given as a fraction of rows would move less, but it would still be learned from whatever rows it sees.

Now that the class is written, the obvious next step is to check whether scikit-learn already does the same job.

def build_with(high_cardinality_step):
    return Pipeline([
        ('prep', ColumnTransformer([
            ('num', Pipeline([('impute', SimpleImputer(strategy='median')),
                              ('scale', StandardScaler())]), NUMERIC),
            ('svc', Pipeline([('count', service_counter), ('scale', StandardScaler())]), SERVICES),
            ('ord', OrdinalEncoder(categories=CONTRACT_ORDER), ORDINAL),
            ('nom', OneHotEncoder(handle_unknown='ignore', drop='if_binary',
                                  sparse_output=False), NOMINAL),
            ('hi',  high_cardinality_step, ['bundle']),
        ], remainder='drop', verbose_feature_names_out=False)),
        ('model', new_lr()),
    ])

options = {
    'our RareCategoryGrouper, then one-hot': Pipeline([
        ('group', RareCategoryGrouper(20)),
        ('encode', OneHotEncoder(handle_unknown='ignore', sparse_output=False))]),
    'OneHotEncoder(min_frequency=20)': OneHotEncoder(handle_unknown='ignore', min_frequency=20,
                                                     sparse_output=False),
    'TargetEncoder': TargetEncoder(random_state=SEED),
    'drop bundle entirely': 'drop',
}
for label, step in options.items():
    print('%-38s %.4f' % (label, cross_val_score(build_with(step), X, y, cv=CV,
                                                 scoring='roc_auc', n_jobs=1).mean()))

our RareCategoryGrouper, then one-hot  0.8416
OneHotEncoder(min_frequency=20)        0.8415
TargetEncoder                          0.8416
drop bundle entirely                   0.8399

OneHotEncoder(min_frequency=20) does the same job with one argument, scoring 0.8415 against the hand-written version's 0.8416. Writing the class was still worth doing once: now you understand what the built-in setting does, and you have a template for the day you need a grouping rule scikit-learn doesn't offer. But shipping the hand-written version when a built-in setting matches it just gives someone code to maintain for no gain.

The more interesting line is at the bottom of that table. Dropping the 322-category column entirely costs only 0.0017. Notebook 4 found bundle was worth a lot on the bare numeric baseline: target-encoded, it scored 0.8250 against 0.8087. Here, with Contract, PaymentMethod, InternetService and num_services already in the model, most of that value is already covered by other columns.

So the choice of high-cardinality encoding, which notebook 4 spent a long section on, has basically stopped mattering. The hand-written grouper, OneHotEncoder(min_frequency=20) and a single target-encoded column all finish within 0.0001 of each other. The only difference left is the 0.0017 lost by dropping the column altogether. A feature isn't worth a fixed amount. It's worth whatever it adds to the columns already there.

Section 5Selection is a pipeline step like any other#

Notebook 8 was built around one demonstration: choosing features on the whole dataset, then cross-validating, reports 74.5% accuracy on a coin flip. Section 2 reproduced it here. The fix is one line in the list

SelectKBest has fit and transform like everything else in this notebook, so it fits between the preprocessor and the model and gets refit automatically on every fold.

for k in [3, 5, 8, 12, 15, 'all']:
    picked = Pipeline([('prep', make_preprocessor()),
                       ('select', SelectKBest(f_classif, k=k)),
                       ('model', new_lr())])
    print('  k = %-4s %.4f' % (k, cross_val_score(picked, X, y, cv=CV,
                                                  scoring='roc_auc', n_jobs=1).mean()))

  k = 3    0.8205
  k = 5    0.8348
  k = 8    0.8395
  k = 12   0.8418
  k = 15   0.8415
  k = all  0.8416

Keeping twelve of the eighteen features scores 0.8418, against 0.8416 for keeping all of them. Selection didn't improve the model in any meaningful way, and it was never going to. Notebook 8 showed that filters help most when a feature set is mostly junk, and this set was built carefully, one column group at a time, with almost no junk left by the time selection ran.

What selection does give you here is six fewer columns to compute, check and explain, for the same score. That's the usual, unexciting and perfectly valid reason to select features, and it's worth saying plainly instead of hunting for a gain that isn't there.

The drop-off below twelve deserves exact numbers too. Going down to eight costs 0.0023, already more than the 0.0017 for dropping bundle entirely. Five costs 0.0070, and three costs 0.0213. These columns are cheap to keep, and they only stop being free below twelve.

Section 6Searching over preprocessing, not just over the model#

Notebook 1 said a later notebook would go deeper into pipelines, and this is one of the biggest payoffs. Once preprocessing is inside the model object, its settings become hyperparameters you can tune like any other

A double-underscore path reaches through the nesting: prep__num__scale points to the scale step of the num group inside the prep transformer. So questions the previous eight notebooks argued about for whole sections become things you can search directly: which scaler, whether the high-cardinality column is worth keeping, and how many features to keep.

From here on, everything uses only the 5,625 development rows. The 1,407 held-out rows are set aside for section 7, and nothing in this section can see them.

One honest caveat: sections 2 to 5 cross-validated on all 7,032 rows, held-out rows included, and some options in this grid, like keeping 8 or 12 features, were chosen after seeing those results. So the hold-out isn't perfectly untouched. The nested cross-validation at the end of section 7 is the cleaner check.

searchable = Pipeline([('prep', make_preprocessor()),
                       ('select', SelectKBest(f_classif, k='all')),
                       ('model', new_lr())])

grid = {
    'prep__num__scale': [StandardScaler(), MinMaxScaler(), RobustScaler()],
    'prep__hi': [TargetEncoder(random_state=SEED),
                 OneHotEncoder(handle_unknown='ignore', min_frequency=20, sparse_output=False),
                 'drop'],
    'select__k': [8, 12, 'all'],
    'model__C': [0.03, 0.3, 3.0],
}

search = GridSearchCV(searchable, grid, cv=CV, scoring='roc_auc', n_jobs=1, refit=True)
search.fit(X_dev, y_dev)

print('%d combinations, best cross-validated AUC %.4f' % (len(search.cv_results_['params']),
                                                          search.best_score_))
for key, value in search.best_params_.items():
    print('   %-18s %s' % (key.split('__', 1)[1], type(value).__name__ if hasattr(value, 'fit') else value))

81 combinations, best cross-validated AUC 0.8425
   C                  3.0
   hi                 TargetEncoder
   num__scale         StandardScaler
   k                  12

The best combination reaches 0.8425: StandardScaler, target encoding for bundle, twelve features and a light regularization penalty (C = 3.0). Note that 3.0 is the largest C in the grid, and the average score rises steadily with C (0.8377, 0.8398, 0.8405). When the best value sits at the edge of a grid, the true best may lie beyond it, so a real search should extend the range. Given how flat these results are, it's unlikely to matter much here.

But look at the overall picture before treating one winning combination as the story.

results = pd.DataFrame(search.cv_results_)
print('all %d combinations span %.4f to %.4f, a range of %.4f'
      % (len(results), results['mean_test_score'].min(), results['mean_test_score'].max(),
         results['mean_test_score'].max() - results['mean_test_score'].min()))
print()
print('average AUC for each setting, over every combination of the others:')
for column, label in [('param_prep__num__scale', 'scaler'), ('param_prep__hi', 'high cardinality'),
                      ('param_select__k', 'k'), ('param_model__C', 'C')]:
    tidy = results[column].astype(str).str.replace(r'\(.*\)', '', regex=True)
    means = results.groupby(tidy)['mean_test_score'].mean()
    print('  %-17s %s' % (label, {name: round(v, 4) for name, v in means.items()}))

all 81 combinations span 0.8354 to 0.8425, a range of 0.0072

average AUC for each setting, over every combination of the others:
  scaler            {'MinMaxScaler': 0.8388, 'RobustScaler': 0.8393, 'StandardScaler': 0.8398}
  high cardinality  {'OneHotEncoder': 0.8394, 'TargetEncoder': 0.8402, 'drop': 0.8383}
  k                 {'12': 0.8398, '8': 0.8382, 'all': 0.84}
  C                 {'0.03': 0.8377, '0.3': 0.8398, '3.0': 0.8405}

fig, ax = plt.subplots(figsize=(9, 3.4))
ranked = results['mean_test_score'].sort_values().to_numpy()
ax.plot(range(len(ranked)), ranked, 'o', ms=4, color='#4c72b0')
ax.axhline(ranked.max(), color='#c44e52', ls='--', lw=1)
ax.annotate('best of 81: %.4f' % ranked.max(), xy=(2, ranked.max()),
            xytext=(2, ranked.max() - 0.0012), fontsize=9, color='#c44e52')
ax.set_ylim(0.830, 0.846)
ax.set_xlabel('the 81 preprocessing and regularisation combinations, worst to best')
ax.set_ylabel('5-fold AUC on the development rows')
ax.set_title('Every combination we searched, on one axis, spanning %.4f' % np.ptp(ranked))
plt.tight_layout()
plt.show()

All 81 combinations score between 0.8354 and 0.8425. The whole search, covering three scalers, three ways to handle a 322-category column, three feature counts and three regularization strengths, moves the score by 0.0072 end to end.

The averages per setting show where that small amount comes from. The biggest single factor is C, worth 0.0028 between its worst and best values, and C is a model setting, not preprocessing. The number of features kept is worth 0.0018. Of the real preprocessing choices, the scaler is worth 0.0010 and the high-cardinality treatment 0.0019.

This is useful to have measured, because on real projects people often spend their effort the other way around. It's easy to fuss over preprocessing choices and then reach for a bigger model to make up the difference. On this dataset, once the columns are handled sensibly, the remaining preprocessing choices barely matter. What really mattered was deciding that Contract is ordinal and that nine service columns can be counted, and no grid search would ever have suggested either of those.

Section 7What you are allowed to claim#

We now have a best cross-validated score of 0.8425, and it's tempting to write it down as the final answer. It isn't, for the same reason notebook 8 ended on

A score used to choose between options is no longer an unbiased estimate. From here on there are three different numbers to keep apart.

def holdout_auc(fitted, columns=None):
    frame = X_hold if columns is None else X_hold[columns]
    return roc_auc_score(y_hold, fitted.predict_proba(frame)[:, 1])

numeric_pipe = Pipeline([('impute', SimpleImputer(strategy='median')),
                         ('scale', StandardScaler()), ('model', new_lr())])
default_pipe = make_full_pipe()

rows = [
    ('numeric only',
     cross_val_score(numeric_pipe, X_dev[NUMERIC], y_dev, cv=CV, scoring='roc_auc', n_jobs=1).mean(),
     holdout_auc(numeric_pipe.fit(X_dev[NUMERIC], y_dev), NUMERIC)),
    ('full pipeline, nothing tuned',
     cross_val_score(default_pipe, X_dev, y_dev, cv=CV, scoring='roc_auc', n_jobs=1).mean(),
     holdout_auc(default_pipe.fit(X_dev, y_dev))),
    ('full pipeline, best of 81',
     search.best_score_,
     holdout_auc(search)),
]
table = pd.DataFrame(rows, columns=['pipeline', 'development CV', 'held-out'])
table['gap'] = table['development CV'] - table['held-out']
print(table.to_string(index=False, float_format=lambda v: '%.4f' % v))

                    pipeline  development CV  held-out     gap
                numeric only          0.8072    0.8111 -0.0039
full pipeline, nothing tuned          0.8419    0.8360  0.0059
   full pipeline, best of 81          0.8425    0.8365  0.0060

Compare the two score columns and the obvious story is that cross-validation was optimistic by about 0.006, and the hold-out corrected it. That story is wrong, and the middle row gives it away. The untuned pipeline shows an almost identical 0.0059 gap, and none of its settings were chosen by looking at that score.

A gap that appears whether or not anything was tuned isn't evidence that tuning caused it. So the next step is to find out how much a single hold-out of 1,407 rows moves on its own, with nothing about the pipeline changing.

gaps = []
for seed in range(8):
    Xd, Xh, yd, yh = train_test_split(X, y, test_size=0.2, random_state=seed, stratify=y)
    pipe = make_full_pipe()
    cv_score = cross_val_score(pipe, Xd, yd, cv=CV, scoring='roc_auc', n_jobs=1).mean()
    held = roc_auc_score(yh, pipe.fit(Xd, yd).predict_proba(Xh)[:, 1])
    gaps.append((seed, cv_score, held, cv_score - held))
    print('  split %d   dev CV %.4f   held-out %.4f   gap %+.4f' % gaps[-1])

gap_values = np.array([g[3] for g in gaps])
print()
print('mean gap %+.4f, standard deviation %.4f, from %+.4f to %+.4f'
      % (gap_values.mean(), gap_values.std(), gap_values.min(), gap_values.max()))
print('spread of the 5-fold means:   %.4f' % np.ptp([g[1] for g in gaps]))
print('spread of the single hold-outs: %.4f' % np.ptp([g[2] for g in gaps]))

  split 0   dev CV 0.8403   held-out 0.8471   gap -0.0068
  split 1   dev CV 0.8450   held-out 0.8265   gap +0.0185
  split 2   dev CV 0.8385   held-out 0.8529   gap -0.0144
  split 3   dev CV 0.8434   held-out 0.8301   gap +0.0134
  split 4   dev CV 0.8429   held-out 0.8375   gap +0.0054
  split 5   dev CV 0.8407   held-out 0.8419   gap -0.0012
  split 6   dev CV 0.8414   held-out 0.8403   gap +0.0011
  split 7   dev CV 0.8384   held-out 0.8530   gap -0.0146

mean gap +0.0002, standard deviation 0.0113, from -0.0146 to +0.0185
spread of the 5-fold means:   0.0065
spread of the single hold-outs: 0.0266

fig, ax = plt.subplots(figsize=(9, 4))
seeds = [g[0] for g in gaps]
ax.plot(seeds, [g[1] for g in gaps], 'o-', label='5-fold CV on the development rows', color='#4c72b0')
ax.plot(seeds, [g[2] for g in gaps], 's-', label='single held-out split', color='#c44e52')
ax.axhline(np.mean([g[1] for g in gaps]), color='#4c72b0', ls=':', lw=1)
ax.set_xlabel('which 20% we happened to hold out (random_state)')
ax.set_ylabel('ROC AUC')
ax.set_title('The same untuned pipeline, eight different hold-outs')
ax.legend(loc='lower right', fontsize=9)
plt.tight_layout()
plt.show()

The gap has an average of +0.0002 and a standard deviation of 0.0113, ranging from −0.0146 to +0.0185, depending only on which 20% of customers landed in the hold-out. It's not a bias. It's noise, and that noise is bigger than every effect the whole grid search in section 6 found.

Two uncomfortable things follow. First, the 0.8360 hold-out figure from earlier was never precise. Quoting four decimal places from one 1,407-row test suggests a precision that doesn't exist. The cross-validated line is about four times steadier, spanning 0.0065 across the eight splits against 0.0266 for the single hold-out, because it averages over five test folds instead of one. That's the real argument for cross-validating instead of trusting one split.

Second, the original question is still open. If that gap is just noise, how much did searching 81 combinations inflate the best score? To answer that properly, the search itself has to run inside each fold. It's the same trick as everything else in this notebook, just applied one level higher.

inner = StratifiedKFold(n_splits=5, shuffle=True, random_state=SEED)
outer = StratifiedKFold(n_splits=5, shuffle=True, random_state=SEED)

nested = cross_val_score(
    GridSearchCV(searchable, grid, cv=inner, scoring='roc_auc', n_jobs=1),
    X, y, cv=outer, scoring='roc_auc', n_jobs=1)

flat = GridSearchCV(searchable, grid, cv=inner, scoring='roc_auc', n_jobs=1).fit(X, y)
untuned = cross_val_score(make_full_pipe(), X, y, cv=outer, scoring='roc_auc', n_jobs=1).mean()

print('nested CV, the search refitted inside every outer fold  %.4f' % nested.mean())
print('  outer folds:', nested.round(4))
print('best-of-81 quoted from the same data it chose on        %.4f' % flat.best_score_)
print('optimism from quoting the winner                        %+.4f' % (flat.best_score_ - nested.mean()))
print()
print('the untuned pipeline, no search at all                  %.4f' % untuned)
print('what the entire search bought                           %+.4f' % (nested.mean() - untuned))

nested CV, the search refitted inside every outer fold  0.8418
  outer folds: [0.8518 0.8439 0.8403 0.8372 0.8358]
best-of-81 quoted from the same data it chose on        0.8419
optimism from quoting the winner                        +0.0001

the untuned pipeline, no search at all                  0.8416
what the entire search bought                           +0.0002

Quoting the winner inflated the score by just +0.0001, and the whole search was worth +0.0002. Both numbers tell us more than a big number would have. Notebook 8's rule explains the first one: the inflation from picking a winner grows with how much the candidates disagree. All 81 candidates here sit inside a narrow 0.0072 band, so there's almost nothing to overfit, and picking the best one is nearly free.

In notebook 8's noise problem, where 2,000 candidate columns disagreed wildly, picking the "best" ones inflated accuracy by 0.27.

The second number is the one to take to work. Eighty-one fits bought two ten-thousandths of AUC, while deciding what the columns meant, back in notebooks 1, 4 and 7, was worth 0.0329. And even that +0.0002 shouldn't be taken seriously as a gain: the five outer folds range from 0.8358 to 0.8518, so a difference in the fourth decimal place is well inside the noise. The honest summary: the search bought nothing measurable, and the column decisions bought over three points of AUC.

Section 8Saving the thing, and what breaks#

A fitted pipeline is one Python object holding the medians, means, category lists, encodings and coefficients. joblib writes it to a file and reads it back, and that's deployment in two lines

It's worth being clear about what those two lines guarantee, and what they don't.

production = make_full_pipe().fit(X_dev, y_dev)
reference = production.predict_proba(X_hold)[:, 1]

path = 'churn_pipeline.joblib'
joblib.dump(production, path)
print('saved %.1f KB' % (os.path.getsize(path) / 1024))

reloaded = joblib.load(path)
print('reloaded predictions bit-for-bit identical:',
      np.array_equal(reference, reloaded.predict_proba(X_hold)[:, 1]))
print('saved from sklearn', sklearn.__version__, '- the file records nothing that enforces this')
os.remove(path)

saved 24.0 KB
reloaded predictions bit-for-bit identical: True
saved from sklearn 1.8.0 - the file records nothing that enforces this

The predictions are identical after saving and loading, which is exactly what you want. The version line printed next to it matters too.

# our service counter uses a named function; the same thing with a lambda cannot be saved
with_lambda = FunctionTransformer(count_services, feature_names_out=lambda t, f: ['num_services'])
try:
    joblib.dump(with_lambda, 'doomed.joblib')
except Exception as err:
    print('%s: %s' % (type(err).__name__, err))
finally:
    # joblib opens the file before it discovers it cannot pickle the lambda
    if os.path.exists('doomed.joblib'):
        os.remove('doomed.joblib')

PicklingError: Can't pickle <function <lambda> at 0x119025760>: it's not found as __main__.<lambda>

This is a common failure, and it tends to show up at the very end, when you save a model after a long training run, instead of during development, when it would be cheap to fix. Pickle, which joblib uses, stores functions and classes by reference, not by value. A lambda has no importable name, so it can't be referenced, and it can't be saved.

Libraries like cloudpickle can save lambdas, but loading still depends on having the same code available. That's why service_feature_names in section 3 is a named function instead of the lambda it wants to be.

The same rule applies to RareCategoryGrouper. The saved file holds a pointer to the class, not the class itself, so whatever loads it needs the same class importable from the same place. A custom transformer defined in a notebook cell won't load anywhere else. Put custom transformers in a proper module and ship it with the model file.

The scikit-learn version matters for a similar reason. The saved file does record which version wrote it, and loading it under a different version raises an InconsistentVersionWarning. But that's only a warning, nothing stops the load. The result can work perfectly, fail loudly or quietly behave differently. So pin your versions next to the model file, every time.

Next, the failures that happen while serving real traffic, which a well-built pipeline handles well.

unseen = X_hold.iloc[[0]].copy()
unseen.loc[:, 'PaymentMethod'] = 'Crypto wallet'          # a category invented after training
print('handle_unknown="ignore": %.4f   (unchanged row: %.4f)'
      % (production.predict_proba(unseen)[0, 1], reference[0]))

strict = Pipeline([
    ('prep', ColumnTransformer([('nom', OneHotEncoder(handle_unknown='error'), NOMINAL)],
                               remainder='drop')),
    ('model', new_lr())]).fit(X_dev, y_dev)
try:
    strict.predict_proba(unseen)
except Exception as err:
    print('handle_unknown="error":  %s: %s' % (type(err).__name__, str(err).split(chr(10))[0]))

handle_unknown="ignore": 0.0211   (unchanged row: 0.0154)
handle_unknown="error":  ValueError: Found unknown categories ['Crypto wallet'] in column 5 during transform
/Library/Frameworks/Python.framework/Versions/3.13/lib/python3.13/site-packages/sklearn/preprocessing/_encoders.py:261: UserWarning: Found unknown categories in columns [5] during transform. These unknown categories will be encoded as all zeros
  warnings.warn(msg, UserWarning)

new_bundle = X_hold.iloc[[0]].copy()
new_bundle.loc[:, 'bundle'] = 'a combination nobody has bought before'
encoder = production.named_steps['prep'].named_transformers_['hi']
print('TargetEncoder value for an unseen bundle: %.4f' % encoder.transform(new_bundle[['bundle']])[0, 0])
print('training churn rate:                      %.4f' % y_dev.mean())
print('prediction for that row:                  %.4f' % production.predict_proba(new_bundle)[0, 1])

TargetEncoder value for an unseen bundle: 0.2658
training churn rate:                      0.2658
prediction for that row:                  0.0210

missing_value = X_hold.iloc[[0]].copy()
missing_value.loc[:, 'TotalCharges'] = np.nan
print('a NaN at serving time, with the imputer:  %.4f'
      % production.predict_proba(missing_value)[0, 1])

no_imputer = Pipeline([
    ('prep', ColumnTransformer([('num', StandardScaler(), NUMERIC)], remainder='drop')),
    ('model', new_lr())]).fit(X_dev, y_dev)
try:
    no_imputer.predict_proba(missing_value)
except Exception as err:
    print('the same row, without the imputer step: %s: %s' % (type(err).__name__, err))

print()
print('column order shuffled: ',
      'identical' if np.array_equal(reference, production.predict_proba(
          X_hold.sample(frac=1, axis=1, random_state=1))[:, 1]) else 'DIFFERENT')
extra = X_hold.copy(); extra['some_new_column'] = 1
print('an extra column:       ',
      'identical' if np.array_equal(reference, production.predict_proba(extra)[:, 1]) else 'DIFFERENT')
try:
    production.predict_proba(X_hold.drop(columns=['gender']))
except Exception as err:
    print('a missing column:       %s: %s' % (type(err).__name__, err))

a NaN at serving time, with the imputer:  0.0064
the same row, without the imputer step: ValueError: Input X contains NaN.
LogisticRegression does not accept missing values encoded as NaN natively. For supervised learning, you might want to consider sklearn.ensemble.HistGradientBoostingClassifier and Regressor which accept missing values encoded as NaNs natively. Alternatively, it is possible to preprocess the data, for instance by using an imputer transformer in a pipeline or drop samples with missing values. See https://scikit-learn.org/stable/modules/impute.html You can find a list of all estimators that handle NaN values at the following page: https://scikit-learn.org/stable/modules/impute.html#estimators-that-handle-nan-values

column order shuffled:  identical
an extra column:        identical
a missing column:       ValueError: columns are missing: {'gender'}

This is the payoff for being explicit earlier. An unseen PaymentMethod becomes a row of zeros instead of a crash, because handle_unknown='ignore' was set on purpose. The prediction moves from 0.0154 to 0.0211, which is just the model relying on everything else it knows about that customer. Scikit-learn raises a warning when this happens, which is the right amount of noise: keep the warning, log it, and alert if it starts happening more often.

An unseen bundle gets the overall training churn rate, 0.2658, the only sensible guess for a category nobody has seen. A NaN in TotalCharges gets the training median. The imputer never filled anything during training, since those rows were filtered out in the first cell. It's there for this moment: without it, the same row fails with Input X contains NaN instead of getting a prediction.

But keeping the model running isn't the same as giving a good answer. Filling the median moved this customer's prediction from 0.0154 to 0.0064. And as notebook 1 warned, a blank TotalCharges on a brand-new customer should be 0, not the median. The imputer protects against crashes. Choosing the right fill value is still your job.

And columns can arrive in any order, because ColumnTransformer selects them by name, not position. An extra column is quietly ignored because of remainder='drop', and a missing column fails right away and names which one is missing. That's a free schema check, just from listing your columns explicitly.