Hyperparameter Optimisation and Interpretability

Week 14 of 16 · Production · 7 days

Week 14 · Production

Hyperparameter Optimisation and Interpretability

Week 14 · Day 1 of 7

What to Tune, and How to Search

Which hyperparameters matter, and why random search beats a grid

Every model so far has been fitted with whatever settings scikit-learn ships. Sometimes that is close to optimal and sometimes it is not, and the difference between a productive week of tuning and a wasted one is almost entirely about which settings you touch.

Not all hyperparameters matter

import numpy as np
import pandas as pd

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

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

d = df.copy()
d['is_monthly'] = (d['contract'] == 'Month-To-Month').astype(int)
d['is_fibre'] = (d['internet_service'] == 'Fibre optic').astype(int)
FEATS = ['tenure_months', 'support_calls', 'monthly_charges',
         'total_charges', 'is_monthly', 'is_fibre']

X_tr, X_te, y_tr, y_te = train_test_split(
    d[FEATS], d['churned'], test_size=0.25, stratify=d['churned'],
    random_state=42)

# Impute after the split, with medians learned from the training rows
# only -- week 4's rule still applies. A Pipeline is the right way to do
# this; these pages need a plain frame to inspect column by column, so
# the two lines are written out instead.
medians = X_tr.median()
X_tr = X_tr.fillna(medians)
X_te = X_te.fillna(medians)
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.model_selection import cross_val_score

def score(**kw):
    m = HistGradientBoostingClassifier(random_state=42, **kw)
    return cross_val_score(m, X_tr, y_tr, cv=5, scoring='roc_auc').mean()

base = score()
print('defaults %.4f\n' % base)
print('%-22s %-18s %9s %9s' % ('parameter', 'values', 'best', 'spread'))
sweeps = {
    'learning_rate': [0.01, 0.05, 0.1, 0.3],
    'max_leaf_nodes': [7, 15, 31, 63],
    'min_samples_leaf': [5, 20, 50, 100],
    'l2_regularization': [0.0, 0.1, 1.0, 10.0],
    'max_features': [0.4, 0.7, 1.0],
}
for name, values in sweeps.items():
    got = [score(**{name: v}) for v in values]
    print('%-22s %-18s %9.4f %9.4f'
          % (name, str(values)[:18], max(got), max(got) - min(got)))

defaults 0.7945

parameter              values                  best    spread
learning_rate          [0.01, 0.05, 0.1,     0.8012    0.0207
max_leaf_nodes         [7, 15, 31, 63]       0.8120    0.0238
min_samples_leaf       [5, 20, 50, 100]      0.8069    0.0144
l2_regularization      [0.0, 0.1, 1.0, 10    0.8032    0.0117
max_features           [0.4, 0.7, 1.0]       0.7992    0.0063

Read the spread column. One or two parameters move the score by something worth having and the rest barely register, and that pattern, not the specific winner, is what generalises. It is why the next two pages are about search strategy rather than about parameter values.

Model	Tune first	Tune if you have budget	Leave alone
Gradient boosting	`learning_rate`, `max_leaf_nodes`, `n_estimators`	`min_samples_leaf`, `l2_regularization`, subsampling	Almost everything else
Random forest	`max_features`, `min_samples_leaf`	`max_depth`	`n_estimators`: set it as high as you can afford
Linear / logistic	`C` or `alpha`	penalty type, `class_weight`	The solver, usually
SVM (RBF)	`C`, `gamma`	kernel choice	-
Neural network	learning rate	width, depth, dropout, batch size	The optimiser, use Adam

`n_estimators` is not a hyperparameter in a forest

More trees in a random forest never makes it worse, it only makes it slower. The variance of the average falls and then flattens. So there is nothing to search: pick the largest number you can afford and move on. In boosting it is entirely different, because each tree corrects the last and too many will overfit. Same argument name, opposite behaviour.

What a grid actually costs

params = {'learning_rate': [0.01, 0.05, 0.1, 0.3],
          'max_leaf_nodes': [7, 15, 31, 63],
          'min_samples_leaf': [5, 20, 50],
          'l2_regularization': [0.0, 0.1, 1.0]}

total = 1
for name, values in params.items():
    total *= len(values)
    print('%-22s %d values, running total %d' % (name, len(values), total))

for folds in [5, 10]:
    fits = total * folds
    print('\n%d-fold: %d fits' % (folds, fits))
    for secs in [0.2, 2.0]:
        print('   at %.1fs per fit: %.1f minutes' % (secs, fits * secs / 60))

learning_rate          4 values, running total 4
max_leaf_nodes         4 values, running total 16
min_samples_leaf       3 values, running total 48
l2_regularization      3 values, running total 144

5-fold: 720 fits
   at 0.2s per fit: 2.4 minutes
   at 2.0s per fit: 24.0 minutes

10-fold: 1440 fits
   at 0.2s per fit: 4.8 minutes
   at 2.0s per fit: 48.0 minutes

Adding one more value to one more parameter multiplies the cost. A grid is the only strategy whose price grows exponentially in the number of things you are curious about, which is a strange property for the default choice to have.

Random search, and why it wins

import numpy as np
import pandas as pd

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

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

d = df.copy()
d['is_monthly'] = (d['contract'] == 'Month-To-Month').astype(int)
d['is_fibre'] = (d['internet_service'] == 'Fibre optic').astype(int)
FEATS = ['tenure_months', 'support_calls', 'monthly_charges',
         'total_charges', 'is_monthly', 'is_fibre']

X_tr, X_te, y_tr, y_te = train_test_split(
    d[FEATS], d['churned'], test_size=0.25, stratify=d['churned'],
    random_state=42)

# Impute after the split, with medians learned from the training rows
# only -- week 4's rule still applies. A Pipeline is the right way to do
# this; these pages need a plain frame to inspect column by column, so
# the two lines are written out instead.
medians = X_tr.median()
X_tr = X_tr.fillna(medians)
X_te = X_te.fillna(medians)
import time
from scipy.stats import loguniform, randint
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.model_selection import GridSearchCV, RandomizedSearchCV

grid = {'learning_rate': [0.01, 0.05, 0.1, 0.3],
        'max_leaf_nodes': [7, 15, 31, 63],
        'min_samples_leaf': [5, 20, 50, 100]}
dists = {'learning_rate': loguniform(0.005, 0.4),
         'max_leaf_nodes': randint(4, 80),
         'min_samples_leaf': randint(3, 120)}

est = HistGradientBoostingClassifier(random_state=42)
for label, search in [
        ('grid, 64 points', GridSearchCV(est, grid, cv=5,
                                         scoring='roc_auc')),
        ('random, 25 draws', RandomizedSearchCV(est, dists, n_iter=25, cv=5,
                                                scoring='roc_auc',
                                                random_state=0))]:
    t = time.time()
    search.fit(X_tr, y_tr)
    print('%-18s best %.4f in %5.1fs (%d fits)'
          % (label, search.best_score_, time.time() - t,
             len(search.cv_results_['params']) * 5))

grid, 64 points    best 0.8155 in  71.8s (320 fits)
random, 25 draws   best 0.8160 in  28.4s (125 fits)

import numpy as np

# Two parameters, only one of which matters. A 5 x 5 grid tries five
# distinct values of each. Twenty-five random draws try twenty-five.
rng = np.random.default_rng(0)
print('%-16s %22s %22s' % ('', 'distinct values tried', 'of the one that matters'))
for n in [9, 16, 25, 36]:
    side = int(round(n ** 0.5))
    print('%-16s %22d %22d'
          % ('grid %dx%d' % (side, side), n, side))
    print('%-16s %22d %22d' % ('random %d' % n, n, n))

                  distinct values tried of the one that matters
grid 3x3                              9                      3
random 9                              9                      9
grid 4x4                             16                      4
random 16                            16                     16
grid 5x5                             25                      5
random 25                            25                     25
grid 6x6                             36                      6
random 36                            36                     36

The argument in one sentence

A grid spends its budget trying the same few values of the parameter that matters, over and over, paired with values of parameters that do not. Random search spends the same budget on distinct values of everything. Since you rarely know in advance which parameters matter, random search is the better default, and it lets you stop at any point rather than at a multiple of the grid size.

Day 1 takeaway

Find out which hyperparameters actually move your score before you search, and search those. Grid cost grows exponentially in the number of parameters and wastes most of its budget on repeats; prefer RandomizedSearchCV with distributions, and use loguniform for anything that spans orders of magnitude.

Week 14 · Day 2 of 7

Spending the Budget Well

Successive halving, early stopping, and Bayesian optimisation from scratch

By Pawan Rama Mali1426 words

Random search still gives every candidate the full budget, including the hopeless ones. Two better ideas: stop bad candidates early, and let the model itself decide when it has had enough.

Successive halving

Successive halving: Evaluate many candidates on a small slice of the data, keep the best fraction, double their resources, and repeat. Bad candidates are eliminated cheaply and the survivors get evaluated properly, so the same total cost covers far more candidates.

import numpy as np
import pandas as pd

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

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

d = df.copy()
d['is_monthly'] = (d['contract'] == 'Month-To-Month').astype(int)
d['is_fibre'] = (d['internet_service'] == 'Fibre optic').astype(int)
FEATS = ['tenure_months', 'support_calls', 'monthly_charges',
         'total_charges', 'is_monthly', 'is_fibre']

X_tr, X_te, y_tr, y_te = train_test_split(
    d[FEATS], d['churned'], test_size=0.25, stratify=d['churned'],
    random_state=42)

# Impute after the split, with medians learned from the training rows
# only -- week 4's rule still applies. A Pipeline is the right way to do
# this; these pages need a plain frame to inspect column by column, so
# the two lines are written out instead.
medians = X_tr.median()
X_tr = X_tr.fillna(medians)
X_te = X_te.fillna(medians)
import time
from scipy.stats import loguniform, randint
from sklearn.experimental import enable_halving_search_cv   # still required
from sklearn.model_selection import HalvingRandomSearchCV
from sklearn.ensemble import HistGradientBoostingClassifier

dists = {'learning_rate': loguniform(0.005, 0.4),
         'max_leaf_nodes': randint(4, 80),
         'min_samples_leaf': randint(3, 120)}

search = HalvingRandomSearchCV(
    HistGradientBoostingClassifier(random_state=42), dists,
    n_candidates=60, factor=3, resource='n_samples',
    # Without this the first round scores on 20 rows, where a 5-fold
    # split can hand a fold a single class and roc_auc returns nan.
    min_resources=250, cv=5,
    scoring='roc_auc', random_state=0)
t = time.time()
search.fit(X_tr, y_tr)

res = pd.DataFrame(search.cv_results_)
print('%-8s %12s %12s %12s' % ('round', 'candidates', 'rows each', 'best'))
for i, g in res.groupby('iter'):
    print('%-8d %12d %12d %12.4f'
          % (i, len(g), g['n_resources'].iloc[0], g['mean_test_score'].max()))
print('\n60 candidates explored in %.1fs, best %.4f'
      % (time.time() - t, search.best_score_))

round      candidates    rows each         best
                60          250       0.7730
                20          750       0.8011
                 7         2250       0.8160

candidates explored in 43.5s, best 0.8160

Sixty candidates screened where day 1's random search managed twenty-five, and the score climbs each round as the survivors get more data, 250 rows, then 750, then all 2,250. It arrives at 0.8160, the same figure random search reached on day 1.

The same answer, for more wall-clock. Measure before you believe

Halving took longer here than the plain random search it matched. That is not a criticism of the algorithm, it is a statement about the dataset: fitting this model on 2,250 rows takes a fraction of a second, so there is nothing to save by fitting it on 250 instead, and the extra rounds of bookkeeping are pure overhead.

Successive halving is built for the opposite situation, hundreds of thousands of rows, or a model that takes minutes per fit, where eliminating forty hopeless candidates cheaply is worth a great deal. Note also min_resources: left at its default this search began on 20 rows, where a five-fold split can hand a fold a single class and the AUC comes back as nan. If your first round is scoring on almost no data, its ranking is noise and the candidates it discards were discarded at random.

Let the model stop itself

import numpy as np
import pandas as pd

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

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

d = df.copy()
d['is_monthly'] = (d['contract'] == 'Month-To-Month').astype(int)
d['is_fibre'] = (d['internet_service'] == 'Fibre optic').astype(int)
FEATS = ['tenure_months', 'support_calls', 'monthly_charges',
         'total_charges', 'is_monthly', 'is_fibre']

X_tr, X_te, y_tr, y_te = train_test_split(
    d[FEATS], d['churned'], test_size=0.25, stratify=d['churned'],
    random_state=42)

# Impute after the split, with medians learned from the training rows
# only -- week 4's rule still applies. A Pipeline is the right way to do
# this; these pages need a plain frame to inspect column by column, so
# the two lines are written out instead.
medians = X_tr.median()
X_tr = X_tr.fillna(medians)
X_te = X_te.fillna(medians)
import time
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.metrics import roc_auc_score

print('%-28s %8s %10s %8s' % ('', 'trees', 'test auc', 'seconds'))
for label, kw in [('max_iter=1000, no stopping',
                   dict(max_iter=1000, early_stopping=False)),
                  ('max_iter=1000, early stopping',
                   dict(max_iter=1000, early_stopping=True,
                        n_iter_no_change=20, validation_fraction=0.15))]:
    t = time.time()
    m = HistGradientBoostingClassifier(random_state=42, **kw).fit(X_tr, y_tr)
    print('%-28s %8d %10.4f %8.1f'
          % (label, m.n_iter_,
             roc_auc_score(y_te, m.predict_proba(X_te)[:, 1]),
             time.time() - t))

                                trees   test auc  seconds
max_iter=1000, no stopping       1000     0.7287      3.9
max_iter=1000, early stopping       47     0.7866      0.2

n_estimators is the one hyperparameter you should almost never search for a boosted model, because the model can find it itself for the price of a single fit. Set it high, turn on early stopping, and spend your search budget on the parameters that cannot be discovered that way.

Bayesian optimisation, from scratch

Bayesian optimisation: Fit a cheap model of the relationship between hyperparameters and score, use it to decide where to look next, evaluate there, update. Instead of sampling blindly it spends the next evaluation where the expected improvement is largest.

import numpy as np
import pandas as pd

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

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

d = df.copy()
d['is_monthly'] = (d['contract'] == 'Month-To-Month').astype(int)
d['is_fibre'] = (d['internet_service'] == 'Fibre optic').astype(int)
FEATS = ['tenure_months', 'support_calls', 'monthly_charges',
         'total_charges', 'is_monthly', 'is_fibre']

X_tr, X_te, y_tr, y_te = train_test_split(
    d[FEATS], d['churned'], test_size=0.25, stratify=d['churned'],
    random_state=42)

# Impute after the split, with medians learned from the training rows
# only -- week 4's rule still applies. A Pipeline is the right way to do
# this; these pages need a plain frame to inspect column by column, so
# the two lines are written out instead.
medians = X_tr.median()
X_tr = X_tr.fillna(medians)
X_te = X_te.fillna(medians)
import numpy as np
from scipy.stats import norm
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern, ConstantKernel
from sklearn.model_selection import cross_val_score

def objective(log_lr):
    m = HistGradientBoostingClassifier(learning_rate=float(10 ** log_lr),
                                       random_state=42)
    return cross_val_score(m, X_tr, y_tr, cv=5, scoring='roc_auc').mean()

rng = np.random.default_rng(0)
seen_x = list(rng.uniform(-2.3, -0.4, size=5))       # 0.005 to 0.4
seen_y = [objective(x) for x in seen_x]
print('%d random starting points, best %.4f\n' % (len(seen_x), max(seen_y)))

XI = 0.002      # how much improvement counts as worth having
grid = np.linspace(-2.3, -0.4, 60).reshape(-1, 1)
for step in range(6):
    gp = GaussianProcessRegressor(ConstantKernel(1.0) * Matern(nu=2.5),
                                  alpha=1e-6, normalize_y=True,
                                  random_state=0)
    gp.fit(np.array(seen_x).reshape(-1, 1), seen_y)
    mu, sd = gp.predict(grid, return_std=True)
    best = max(seen_y) + XI            # XI keeps it exploring
    z = (mu - best) / np.maximum(sd, 1e-9)
    ei = (mu - best) * norm.cdf(z) + sd * norm.pdf(z)   # expected improvement
    nxt = float(grid[int(ei.argmax()), 0])
    seen_x.append(nxt)
    seen_y.append(objective(nxt))
    print('step %d: tried lr=%.4f -> %.4f   (best so far %.4f)'
          % (step + 1, 10 ** nxt, seen_y[-1], max(seen_y)))

i = int(np.argmax(seen_y))
print('\nbest learning rate %.4f scoring %.4f after %d evaluations'
      % (10 ** seen_x[i], seen_y[i], len(seen_y)))

5 random starting points, best 0.8045

step 1: tried lr=0.0297 -> 0.8054   (best so far 0.8054)
step 2: tried lr=0.3981 -> 0.7788   (best so far 0.8054)
step 3: tried lr=0.0430 -> 0.8031   (best so far 0.8054)
step 4: tried lr=0.0221 -> 0.8053   (best so far 0.8054)
step 5: tried lr=0.0105 -> 0.8007   (best so far 0.8054)
step 6: tried lr=0.0579 -> 0.7969   (best so far 0.8054)

best learning rate 0.0297 scoring 0.8054 after 11 evaluations

When this is worth the complexity

Bayesian optimisation earns its keep when a single evaluation is expensive, minutes or hours, because then it is worth thinking hard about where to spend the next one. When a fit takes a second, the thinking costs more than the fitting and random search with more draws wins on wall-clock. Libraries such as Optuna implement this properly, with pruning and parallelism; the version above is here so the idea is not a black box.

Day 2 takeaway

Successive halving eliminates weak candidates on cheap evaluations and explores far more of the space per unit of time. Early stopping removes n_estimators from the search entirely. Bayesian optimisation is worth it only when each evaluation is genuinely expensive.

Week 14 · Day 3 of 7

Tuning the Whole Pipeline

Preprocessing as a hyperparameter, caching, and the bias in best_score_

By Pawan Rama Mali1109 words

The choices made before the model sees the data, which imputer, which encoder, whether to scale, are hyperparameters too. They are usually decided once, by hand, and never revisited, which is odd given how much they can matter.

Tuning the preprocessing

import numpy as np
import pandas as pd

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

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

prep = ColumnTransformer([
    ('num', Pipeline([('i', SimpleImputer(strategy='median')),
                      ('s', StandardScaler())]), NUM),
    ('cat', Pipeline([('i', SimpleImputer(strategy='most_frequent')),
                      ('o', OneHotEncoder(handle_unknown='ignore'))]), CAT),
])

X, y = df[NUM + CAT], df['churned']
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.25,
                                          stratify=y, random_state=42)
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import RandomizedSearchCV
from scipy.stats import loguniform

pipe = Pipeline([('prep', prep),
                 ('clf', LogisticRegression(max_iter=2000))])

space = {
    'prep__num__i__strategy': ['median', 'mean'],
    'prep__num__s': [StandardScaler(), 'passthrough'],
    'prep__cat__o__min_frequency': [1, 20, 60],
    'clf__C': loguniform(0.01, 100),
}
search = RandomizedSearchCV(pipe, space, n_iter=20, cv=5,
                            scoring='roc_auc', random_state=0)
search.fit(X_tr, y_tr)
print('best cross-validated AUC %.4f' % search.best_score_)
for k, v in sorted(search.best_params_.items()):
    print('  %-28s %s' % (k, v))

best cross-validated AUC 0.8219
  clf__C                       0.19795388574374428
  prep__cat__o__min_frequency  20
  prep__num__i__strategy       median
  prep__num__s                 StandardScaler()

The double underscore is a path

prep__num__i__strategy reads as: the step named prep, inside it the transformer named num, inside that the step named i, and its strategy argument. Any nesting depth works, and passing 'passthrough' in place of a step is how you make the existence of that step itself a thing to be tuned.

Caching the repeated work, and when not to

import numpy as np
import pandas as pd

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

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

prep = ColumnTransformer([
    ('num', Pipeline([('i', SimpleImputer(strategy='median')),
                      ('s', StandardScaler())]), NUM),
    ('cat', Pipeline([('i', SimpleImputer(strategy='most_frequent')),
                      ('o', OneHotEncoder(handle_unknown='ignore'))]), CAT),
])

X, y = df[NUM + CAT], df['churned']
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.25,
                                          stratify=y, random_state=42)
import shutil
import tempfile
import time
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import GridSearchCV

space = {'clf__C': [0.01, 0.1, 1.0, 10.0, 100.0]}
cache = tempfile.mkdtemp()

for label, memory in [('no cache', None), ('memory=cache', cache)]:
    pipe = Pipeline([('prep', prep),
                     ('clf', LogisticRegression(max_iter=2000))],
                    memory=memory)
    t = time.time()
    GridSearchCV(pipe, space, cv=5, scoring='roc_auc').fit(X_tr, y_tr)
    print('%-14s %.2fs' % (label, time.time() - t))

shutil.rmtree(cache, ignore_errors=True)
print('\nthe preprocessing is identical for every value of C, so the')
print('cache turns 25 fits of it into 5 -- plus 25 trips to disk.')

no cache       0.44s
memory=cache   0.88s

the preprocessing is identical for every value of C, so the
cache turns 25 fits of it into 5 -- plus 25 trips to disk.

Measure it; here the cache made things slower

Caching trades recomputation for serialisation and disk I/O. This ColumnTransformer is an imputer, a scaler and a one-hot encoder over 2,250 rows. It costs a few milliseconds, and writing the result to disk and reading it back costs more than that. So the cached version loses.

It wins, and wins enormously, when the preprocessing is genuinely expensive: fitting a TfidfVectorizer over a large corpus, a KNNImputer, a target encoder, anything reading from disk. The rule is the same one as for successive halving: the optimisation is only an optimisation if you timed it.

How much of the gain is real?

A search reports the score of the best candidate it found. That figure is the maximum of many noisy numbers, and a maximum is biased upwards: some of what the winner won by is genuine and some of it is the fold split happening to suit it.

import numpy as np
import pandas as pd

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

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

d = df.copy()
d['is_monthly'] = (d['contract'] == 'Month-To-Month').astype(int)
d['is_fibre'] = (d['internet_service'] == 'Fibre optic').astype(int)
FEATS = ['tenure_months', 'support_calls', 'monthly_charges',
         'total_charges', 'is_monthly', 'is_fibre']

X_tr, X_te, y_tr, y_te = train_test_split(
    d[FEATS], d['churned'], test_size=0.25, stratify=d['churned'],
    random_state=42)

# Impute after the split, with medians learned from the training rows
# only -- week 4's rule still applies. A Pipeline is the right way to do
# this; these pages need a plain frame to inspect column by column, so
# the two lines are written out instead.
medians = X_tr.median()
X_tr = X_tr.fillna(medians)
X_te = X_te.fillna(medians)
import numpy as np
from scipy.stats import loguniform, randint
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.model_selection import RandomizedSearchCV
from sklearn.metrics import roc_auc_score

dists = {'learning_rate': loguniform(0.005, 0.4),
         'max_leaf_nodes': randint(4, 80),
         'min_samples_leaf': randint(3, 120)}

print('%-12s %14s %12s %10s' % ('candidates', 'best cv score', 'test score',
                                'gap'))
for n in [5, 20, 60]:
    s = RandomizedSearchCV(HistGradientBoostingClassifier(random_state=42),
                           dists, n_iter=n, cv=5, scoring='roc_auc',
                           random_state=0).fit(X_tr, y_tr)
    test = roc_auc_score(y_te, s.predict_proba(X_te)[:, 1])
    print('%-12d %14.4f %12.4f %10.4f'
          % (n, s.best_score_, test, s.best_score_ - test))

candidates    best cv score   test score        gap
                  0.8118       0.7976     0.0142
                 0.8160       0.8038     0.0121
                 0.8160       0.8038     0.0121

The gap is positive at every budget: the cross-validated figure is optimistic by more than a point of AUC each time. Note also that trying sixty candidates instead of twenty bought nothing at all here, the same winner, the same score. A search that has stopped improving is telling you the ceiling is set by the features, not the settings, and that is a signal to go back to week 8 rather than to widen the search.

Never report `best_score_` as your result

It is the winner of a competition you ran on your own validation folds, so it carries the same selection bias as any other "best of many" statistic. Report the held-out score, or the nested cross-validation score from week 7 day 1 if you have no held-out set to spare. best_score_ is for comparing candidates against each other, and for nothing else.

Day 3 takeaway

Preprocessing choices belong inside the search, addressed through the double-underscore path, with 'passthrough' to make a step optional. Add memory= so shared preprocessing is not recomputed. And always report a score from data the search never touched.

Week 14 · Day 4 of 7

Partial Dependence and ICE

How the prediction responds to a feature, and when the average lies

By Pawan Rama Mali1575 words

A tuned model that nobody understands is difficult to defend and impossible to debug. Week 6 established which features a model relies on. This is the next question: how does the prediction change as a feature changes?

Partial dependence, by hand first

Partial dependence: Take every row in the data, set one feature to a fixed value, predict, and average. Repeat across a range of values. The result is the model's average response to that feature with everything else held as it is in the data.

import numpy as np
import pandas as pd

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

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

d = df.copy()
d['is_monthly'] = (d['contract'] == 'Month-To-Month').astype(int)
d['is_fibre'] = (d['internet_service'] == 'Fibre optic').astype(int)
FEATS = ['tenure_months', 'support_calls', 'monthly_charges',
         'total_charges', 'is_monthly', 'is_fibre']

X_tr, X_te, y_tr, y_te = train_test_split(
    d[FEATS], d['churned'], test_size=0.25, stratify=d['churned'],
    random_state=42)

# Impute after the split, with medians learned from the training rows
# only -- week 4's rule still applies. A Pipeline is the right way to do
# this; these pages need a plain frame to inspect column by column, so
# the two lines are written out instead.
medians = X_tr.median()
X_tr = X_tr.fillna(medians)
X_te = X_te.fillna(medians)
import numpy as np
from sklearn.ensemble import HistGradientBoostingClassifier

m = HistGradientBoostingClassifier(random_state=42).fit(X_tr, y_tr)

grid = [1, 6, 12, 24, 36, 48, 60, 72]
print('%14s %14s' % ('tenure_months', 'mean P(churn)'))
for value in grid:
    copy = X_tr.copy()
    copy['tenure_months'] = value           # everybody, at this tenure
    print('%14d %14.4f' % (value, m.predict_proba(copy)[:, 1].mean()))

 tenure_months  mean P(churn)
             1         0.3687
             6         0.2390
            12         0.3096
            24         0.1227
            36         0.1526
            48         0.0202
            60         0.0079
            72         0.0082

import numpy as np
import pandas as pd

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

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

d = df.copy()
d['is_monthly'] = (d['contract'] == 'Month-To-Month').astype(int)
d['is_fibre'] = (d['internet_service'] == 'Fibre optic').astype(int)
FEATS = ['tenure_months', 'support_calls', 'monthly_charges',
         'total_charges', 'is_monthly', 'is_fibre']

X_tr, X_te, y_tr, y_te = train_test_split(
    d[FEATS], d['churned'], test_size=0.25, stratify=d['churned'],
    random_state=42)

# Impute after the split, with medians learned from the training rows
# only -- week 4's rule still applies. A Pipeline is the right way to do
# this; these pages need a plain frame to inspect column by column, so
# the two lines are written out instead.
medians = X_tr.median()
X_tr = X_tr.fillna(medians)
X_te = X_te.fillna(medians)
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.inspection import partial_dependence

m = HistGradientBoostingClassifier(random_state=42).fit(X_tr, y_tr)

for method in ['decision_function', 'predict_proba']:
    res = partial_dependence(m, X_tr, ['tenure_months'],
                             grid_resolution=5, kind='average',
                             # 'recursion' is the fast tree-specific path
                             # and only speaks log-odds; asking for
                             # probabilities means doing it the slow way.
                             method='brute', response_method=method)
    print('response_method=%s' % method)
    print('  ' + '  '.join('%8.4f' % a for a in res['average'][0]))
print('\ngrid: ' + '  '.join('%8.1f' % v
                             for v in res['grid_values'][0]))

response_method=decision_function
   -1.2841   -1.9846   -3.0956   -5.3464   -5.9621
response_method=predict_proba
    0.3101    0.2245    0.1138    0.0141    0.0079

grid:      4.0      19.2      34.5      49.8      65.0

Check the scale before you read the numbers

For a classifier, response_method='auto' uses decision_function when the estimator has one, so the default output is in log-odds, not probability, and will not match the hand-rolled version above. Log-odds is the better scale for seeing shape, because it is the scale the model is additive on. Probability is the better scale for showing anybody else. Pick deliberately, and label the axis.

The average can hide the story

ICE curve: Individual Conditional Expectation: the same calculation done per row instead of averaged. One line per customer. Partial dependence is the mean of these lines, and a mean is a poor summary when the lines disagree about direction.

import numpy as np
import pandas as pd

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

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

d = df.copy()
d['is_monthly'] = (d['contract'] == 'Month-To-Month').astype(int)
d['is_fibre'] = (d['internet_service'] == 'Fibre optic').astype(int)
FEATS = ['tenure_months', 'support_calls', 'monthly_charges',
         'total_charges', 'is_monthly', 'is_fibre']

X_tr, X_te, y_tr, y_te = train_test_split(
    d[FEATS], d['churned'], test_size=0.25, stratify=d['churned'],
    random_state=42)

# Impute after the split, with medians learned from the training rows
# only -- week 4's rule still applies. A Pipeline is the right way to do
# this; these pages need a plain frame to inspect column by column, so
# the two lines are written out instead.
medians = X_tr.median()
X_tr = X_tr.fillna(medians)
X_te = X_te.fillna(medians)
import numpy as np
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.inspection import partial_dependence

m = HistGradientBoostingClassifier(random_state=42).fit(X_tr, y_tr)
res = partial_dependence(m, X_tr, ['monthly_charges'],
                         grid_resolution=6, kind='both')
grid = res['grid_values'][0]
ice = res['individual'][0]

slope = ice[:, -1] - ice[:, 0]
print('%d customers rise with price, %d fall'
      % (int((slope > 0).sum()), int((slope < 0).sum())))
print('\n%14s %10s %10s %10s' % ('charges', 'average', 'p10', 'p90'))
for j, g in enumerate(grid):
    col = ice[:, j]
    print('%14.1f %10.4f %10.4f %10.4f'
          % (g, col.mean(), np.percentile(col, 10),
             np.percentile(col, 90)))

2123 customers rise with price, 127 fall

       charges    average        p10        p90
          15.0     0.1230     0.0017     0.3105
          30.2     0.1276     0.0022     0.3310
          45.4     0.1177     0.0007     0.3135
          60.6     0.2636     0.0014     0.6268
          75.8     0.2301     0.0034     0.5721
          91.0     0.3635     0.0063     0.8330

Two features at once

import numpy as np
import pandas as pd

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

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

d = df.copy()
d['is_monthly'] = (d['contract'] == 'Month-To-Month').astype(int)
d['is_fibre'] = (d['internet_service'] == 'Fibre optic').astype(int)
FEATS = ['tenure_months', 'support_calls', 'monthly_charges',
         'total_charges', 'is_monthly', 'is_fibre']

X_tr, X_te, y_tr, y_te = train_test_split(
    d[FEATS], d['churned'], test_size=0.25, stratify=d['churned'],
    random_state=42)

# Impute after the split, with medians learned from the training rows
# only -- week 4's rule still applies. A Pipeline is the right way to do
# this; these pages need a plain frame to inspect column by column, so
# the two lines are written out instead.
medians = X_tr.median()
X_tr = X_tr.fillna(medians)
X_te = X_te.fillna(medians)
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.inspection import partial_dependence

m = HistGradientBoostingClassifier(random_state=42).fit(X_tr, y_tr)
# Two names in one list asks for the joint surface, not two curves.
res = partial_dependence(m, X_tr, ['tenure_months', 'is_monthly'],
                         grid_resolution=6, kind='average')
tenure, monthly = res['grid_values']
surface = res['average'][0]

print('rows = tenure, columns = contract type')
print('%10s %12s %12s %12s' % ('tenure', 'longer term', 'month-to-month',
                               'difference'))
for i, t in enumerate(tenure):
    a, b = surface[i, 0], surface[i, 1]
    print('%10.1f %12.4f %12.4f %12.4f' % (t, a, b, b - a))

rows = tenure, columns = contract type
    tenure  longer term month-to-month   difference
       4.0      -0.9477       0.6599       1.6076
      16.2      -1.7094       0.2208       1.9301
      28.4      -3.9701      -0.0364       3.9337
      40.6      -3.9225      -1.7032       2.2193
      52.8      -4.6525      -2.6154       2.0370
      65.0      -5.9758      -4.0526       1.9232

If the difference column were constant, the two features would act independently and the model would be additive in them. Where it changes, there is an interaction: the effect of tenure depends on the contract. That is precisely the structure a linear model cannot represent without being told, and the reason week 8 spent a day on interaction features.

Partial dependence averages over combinations that do not exist

Setting tenure_months=1 for every customer includes customers whose total_charges is £3,000, which is a person who has spent three thousand pounds in their first month. The model is being asked about a region of the input space with no data in it, and its answer there is extrapolation dressed as insight. The stronger the correlation between your features, the less trustworthy the curve, check the correlation before you present one.

import numpy as np
import pandas as pd

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

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

d = df.copy()
d['is_monthly'] = (d['contract'] == 'Month-To-Month').astype(int)
d['is_fibre'] = (d['internet_service'] == 'Fibre optic').astype(int)
FEATS = ['tenure_months', 'support_calls', 'monthly_charges',
         'total_charges', 'is_monthly', 'is_fibre']

X_tr, X_te, y_tr, y_te = train_test_split(
    d[FEATS], d['churned'], test_size=0.25, stratify=d['churned'],
    random_state=42)

# Impute after the split, with medians learned from the training rows
# only -- week 4's rule still applies. A Pipeline is the right way to do
# this; these pages need a plain frame to inspect column by column, so
# the two lines are written out instead.
medians = X_tr.median()
X_tr = X_tr.fillna(medians)
X_te = X_te.fillna(medians)
print('correlation with tenure_months:')
print(X_tr.corr()['tenure_months'].drop('tenure_months').round(3)
      .sort_values(key=abs, ascending=False).to_string())
print('\ntotal_charges is nearly determined by tenure, so a partial')
print('dependence plot for one of them holds the other at impossible values.')

correlation with tenure_months:
total_charges      0.829
is_monthly        -0.463
support_calls      0.033
is_fibre           0.014
monthly_charges    0.007

total_charges is nearly determined by tenure, so a partial
dependence plot for one of them holds the other at impossible values.

Day 4 takeaway

Partial dependence shows the model's average response to a feature; ICE curves show whether that average represents anybody. Use a two-way plot to see interactions. And treat both with suspicion when the feature is strongly correlated with another, because the calculation walks off the edge of your data.

Week 14 · Day 5 of 7

SHAP

Per-prediction contributions that add up, and what they cost

By Pawan Rama Mali1153 words

Partial dependence describes the model. It does not tell a particular customer why they were flagged, and that is usually the question that gets asked.

Shapley value: From cooperative game theory: the fair division of a payout among players who contributed unequally. Treat the prediction as the payout and each feature as a player, and you get a per-feature contribution for a single prediction, unique, given a few reasonable axioms.

One prediction, explained

import numpy as np
import pandas as pd

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

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

d = df.copy()
d['is_monthly'] = (d['contract'] == 'Month-To-Month').astype(int)
d['is_fibre'] = (d['internet_service'] == 'Fibre optic').astype(int)
FEATS = ['tenure_months', 'support_calls', 'monthly_charges',
         'total_charges', 'is_monthly', 'is_fibre']

X_tr, X_te, y_tr, y_te = train_test_split(
    d[FEATS], d['churned'], test_size=0.25, stratify=d['churned'],
    random_state=42)

# Impute after the split, with medians learned from the training rows
# only -- week 4's rule still applies. A Pipeline is the right way to do
# this; these pages need a plain frame to inspect column by column, so
# the two lines are written out instead.
medians = X_tr.median()
X_tr = X_tr.fillna(medians)
X_te = X_te.fillna(medians)
import shap
from sklearn.ensemble import HistGradientBoostingClassifier

m = HistGradientBoostingClassifier(random_state=42).fit(X_tr, y_tr)
explainer = shap.TreeExplainer(m)
sv = explainer(X_te.iloc[:50])

i = int(m.predict_proba(X_te.iloc[:50])[:, 1].argmax())   # most at risk
row = X_te.iloc[i]
print('customer predicted at %.3f probability of churn'
      % m.predict_proba(X_te.iloc[[i]])[0, 1])
print('\n%-18s %10s %12s' % ('feature', 'value', 'contribution'))
order = np.argsort(-np.abs(sv.values[i]))
for j in order:
    print('%-18s %10.2f %+12.4f'
          % (X_te.columns[j], row.iloc[j], sv.values[i][j]))
print('%-18s %10s %+12.4f' % ('base value', '', sv.base_values[i]))
print('%-18s %10s %+12.4f'
      % ('sum', '', sv.values[i].sum() + sv.base_values[i]))
print('model log-odds output %+.4f' % m.decision_function(X_te.iloc[[i]])[0])

customer predicted at 0.814 probability of churn

feature                 value contribution
is_monthly               1.00      +2.0362
support_calls            4.00      +1.4048
monthly_charges         88.05      +1.0216
total_charges         2590.66      -0.5712
tenure_months           28.00      -0.1467
is_fibre                 1.00      +0.0414
base value                         -2.3128
sum                                +1.4733
model log-odds output +1.4733

Additivity is the property that makes this usable

The contributions plus the base value equal the model's output exactly, not approximately. That is what lets you put the numbers in front of somebody and say "these five things, adding to this" without hand-waving. Note the scale: for a classifier, contributions are in log-odds, not probability, so they add up on a scale most audiences do not think in. Convert carefully, or explain the scale.

Global importance, assembled from local explanations

import numpy as np
import pandas as pd

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

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

d = df.copy()
d['is_monthly'] = (d['contract'] == 'Month-To-Month').astype(int)
d['is_fibre'] = (d['internet_service'] == 'Fibre optic').astype(int)
FEATS = ['tenure_months', 'support_calls', 'monthly_charges',
         'total_charges', 'is_monthly', 'is_fibre']

X_tr, X_te, y_tr, y_te = train_test_split(
    d[FEATS], d['churned'], test_size=0.25, stratify=d['churned'],
    random_state=42)

# Impute after the split, with medians learned from the training rows
# only -- week 4's rule still applies. A Pipeline is the right way to do
# this; these pages need a plain frame to inspect column by column, so
# the two lines are written out instead.
medians = X_tr.median()
X_tr = X_tr.fillna(medians)
X_te = X_te.fillna(medians)
import shap
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.inspection import permutation_importance

m = HistGradientBoostingClassifier(random_state=42).fit(X_tr, y_tr)
sv = shap.TreeExplainer(m)(X_te)
shap_rank = np.abs(sv.values).mean(axis=0)

perm = permutation_importance(m, X_te, y_te, n_repeats=10,
                              random_state=0, scoring='roc_auc')

print('%-18s %14s %16s' % ('feature', 'mean |shap|', 'permutation'))
for j in np.argsort(-shap_rank):
    print('%-18s %14.4f %16.4f'
          % (X_te.columns[j], shap_rank[j], perm.importances_mean[j]))

feature               mean |shap|      permutation
is_monthly                 1.3276           0.0984
tenure_months              0.9738           0.0804
monthly_charges            0.5814           0.0245
total_charges              0.5090          -0.0027
support_calls              0.4216           0.0191
is_fibre                   0.0885          -0.0015

The two rankings usually agree, and when they disagree it is informative rather than alarming: permutation importance measures how much the score falls when a feature is scrambled, so it is about predictive usefulness. Mean absolute SHAP measures how much the prediction moves, so a feature can swing predictions strongly in both directions and still not improve accuracy.

What it costs

import numpy as np
import pandas as pd

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

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

d = df.copy()
d['is_monthly'] = (d['contract'] == 'Month-To-Month').astype(int)
d['is_fibre'] = (d['internet_service'] == 'Fibre optic').astype(int)
FEATS = ['tenure_months', 'support_calls', 'monthly_charges',
         'total_charges', 'is_monthly', 'is_fibre']

X_tr, X_te, y_tr, y_te = train_test_split(
    d[FEATS], d['churned'], test_size=0.25, stratify=d['churned'],
    random_state=42)

# Impute after the split, with medians learned from the training rows
# only -- week 4's rule still applies. A Pipeline is the right way to do
# this; these pages need a plain frame to inspect column by column, so
# the two lines are written out instead.
medians = X_tr.median()
X_tr = X_tr.fillna(medians)
X_te = X_te.fillna(medians)
import time
import shap
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.linear_model import LogisticRegression

m = HistGradientBoostingClassifier(random_state=42).fit(X_tr, y_tr)
lr = LogisticRegression(max_iter=2000).fit(X_tr, y_tr)

t = time.time()
shap.TreeExplainer(m)(X_te)
tree = time.time() - t

t = time.time()
shap.KernelExplainer(lr.predict_proba,
                     shap.sample(X_tr, 50, random_state=0)).shap_values(
    X_te.iloc[:25], nsamples=100, silent=True)
kernel = time.time() - t

print('%-16s %6s %8s %14s' % ('', 'rows', 'seconds', 'ms per row'))
print('%-16s %6d %8.2f %14.2f'
      % ('TreeExplainer', len(X_te), tree, 1000 * tree / len(X_te)))
print('%-16s %6d %8.2f %14.2f'
      % ('KernelExplainer', 25, kernel, 1000 * kernel / 25))
print('\ncompare the last column, not the middle one.')

                   rows  seconds     ms per row
TreeExplainer       750     0.22           0.29
KernelExplainer      25     0.06           2.59

compare the last column, not the middle one.

Per row, the sampling explainer is an order of magnitude more expensive, and that is with only six features, a hundred samples and a fifty-row background set. All three of those have to grow for a real problem, and the cost grows with them. TreeExplainer avoids all of it by exploiting the tree structure to compute the exact answer directly, which is why tree models are the ones people actually explain at scale.

Correlated features share credit arbitrarily

If two features carry the same information, any split of the credit between them produces the same prediction, so the axioms do not pin down a unique answer, the implementation picks one. A stakeholder reading "total_charges contributed most" may act on it, when tenure_months would have served identically. Check your correlations before you present attributions, and consider reporting correlated features as a group.

Day 5 takeaway

SHAP gives per-prediction contributions that sum exactly to the output, which is what makes individual decisions explainable. Aggregate them for a global view, but remember they are in log-odds for a classifier, that TreeExplainer is the only fast exact option, and that correlated features make the split of credit ambiguous.

Week 14 · Day 6 of 7

Counterfactuals, Surrogates and Causality

What would have to change, readable approximations, and the limit

By Pawan Rama Mali1456 words

Two more tools, both aimed at the gap between what a model computes and what a person can act on.

Counterfactuals: what would have to change

Counterfactual explanation: The smallest change to a customer's features that would flip the model's decision. It answers the question people actually ask, not "why" but "what would I have to do", and it is the form of explanation regulators tend to have in mind.

import numpy as np
import pandas as pd

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

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

d = df.copy()
d['is_monthly'] = (d['contract'] == 'Month-To-Month').astype(int)
d['is_fibre'] = (d['internet_service'] == 'Fibre optic').astype(int)
FEATS = ['tenure_months', 'support_calls', 'monthly_charges',
         'total_charges', 'is_monthly', 'is_fibre']

X_tr, X_te, y_tr, y_te = train_test_split(
    d[FEATS], d['churned'], test_size=0.25, stratify=d['churned'],
    random_state=42)

# Impute after the split, with medians learned from the training rows
# only -- week 4's rule still applies. A Pipeline is the right way to do
# this; these pages need a plain frame to inspect column by column, so
# the two lines are written out instead.
medians = X_tr.median()
X_tr = X_tr.fillna(medians)
X_te = X_te.fillna(medians)
import numpy as np
from sklearn.ensemble import HistGradientBoostingClassifier

m = HistGradientBoostingClassifier(random_state=42).fit(X_tr, y_tr)
proba = m.predict_proba(X_te)[:, 1]
i = int(proba.argmax())
row = X_te.iloc[[i]]
print('customer starts at P(churn) = %.3f\n' % proba[i])

# Only the levers the business can actually pull.
moves = {'monthly_charges': [-20, -10, -5],
         'support_calls': [-3, -2, -1],
         'is_monthly': [0]}

print('%-20s %10s %12s %10s' % ('change', 'new value', 'P(churn)', 'drop'))
for feature, deltas in moves.items():
    for delta in deltas:
        candidate = row.copy()
        new = (delta if feature == 'is_monthly'
               else max(0, float(row[feature].iloc[0]) + delta))
        candidate[feature] = new
        p = m.predict_proba(candidate)[0, 1]
        print('%-20s %10.1f %12.3f %10.3f'
              % ('%s %+g' % (feature, delta) if feature != 'is_monthly'
                 else 'move off monthly', new, p, proba[i] - p))

customer starts at P(churn) = 0.923

change                new value     P(churn)       drop
monthly_charges -20        69.6        0.693      0.230
monthly_charges -10        79.6        0.619      0.304
monthly_charges -5         84.6        0.781      0.142
support_calls -3            0.0        0.545      0.378
support_calls -2            0.0        0.545      0.378
support_calls -1            1.0        0.871      0.052
move off monthly            0.0        0.199      0.724

Moving this customer off a month-to-month contract is worth more than any price change on offer, which is an actionable finding rather than a statistical one. Notice too that the price column is not monotone: a £10 reduction helps more than a £20 one. That is not a bug. A boosted tree is a step function, so moving a customer across one threshold can help while moving them further lands them in a different leaf. Search the whole range rather than assuming more of a good thing is better.

Only offer changes somebody can make

A counterfactual saying "if this customer had been with us for forty more months" is useless, because nobody can act on it. Restrict the search to actionable features, price, contract type, the number of unresolved support calls, and the output becomes a retention offer rather than a curiosity. It also keeps you away from features that would be unlawful to act on.

Surrogate models

import numpy as np
import pandas as pd

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

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

d = df.copy()
d['is_monthly'] = (d['contract'] == 'Month-To-Month').astype(int)
d['is_fibre'] = (d['internet_service'] == 'Fibre optic').astype(int)
FEATS = ['tenure_months', 'support_calls', 'monthly_charges',
         'total_charges', 'is_monthly', 'is_fibre']

X_tr, X_te, y_tr, y_te = train_test_split(
    d[FEATS], d['churned'], test_size=0.25, stratify=d['churned'],
    random_state=42)

# Impute after the split, with medians learned from the training rows
# only -- week 4's rule still applies. A Pipeline is the right way to do
# this; these pages need a plain frame to inspect column by column, so
# the two lines are written out instead.
medians = X_tr.median()
X_tr = X_tr.fillna(medians)
X_te = X_te.fillna(medians)
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.tree import DecisionTreeRegressor, export_text

m = HistGradientBoostingClassifier(random_state=42).fit(X_tr, y_tr)
target = m.predict_proba(X_tr)[:, 1]        # the black box's own output

print('%8s %16s' % ('depth', 'fidelity (R2)'))
for depth in [2, 3, 4, 6]:
    s = DecisionTreeRegressor(max_depth=depth,
                              random_state=0).fit(X_tr, target)
    print('%8d %16.4f' % (depth, s.score(X_tr, target)))

s = DecisionTreeRegressor(max_depth=3, random_state=0).fit(X_tr, target)
print('\n' + export_text(s, feature_names=list(X_tr.columns),
                         decimals=2)[:1100])

   depth    fidelity (R2)
       2           0.4484
       3           0.5437
       4           0.6177
       6           0.7213

|--- is_monthly <= 0.50
|   |--- tenure_months <= 17.50
|   |   |--- support_calls <= 2.50
|   |   |   |--- value: [0.18]
|   |   |--- support_calls >  2.50
|   |   |   |--- value: [0.43]
|   |--- tenure_months >  17.50
|   |   |--- monthly_charges <= 75.74
|   |   |   |--- value: [0.02]
|   |   |--- monthly_charges >  75.74
|   |   |   |--- value: [0.07]
|--- is_monthly >  0.50
|   |--- monthly_charges <= 50.98
|   |   |--- total_charges <= 237.06
|   |   |   |--- value: [0.34]
|   |   |--- total_charges >  237.06
|   |   |   |--- value: [0.20]
|   |--- monthly_charges >  50.98
|   |   |--- tenure_months <= 32.50
|   |   |   |--- value: [0.52]
|   |   |--- tenure_months >  32.50
|   |   |   |--- value: [0.09]

Fidelity is the number that matters, and it is usually omitted

A surrogate is an explanation of the model only to the extent that it agrees with the model. Quote the R² whenever you show one: at 0.9 it is a fair summary, at 0.5 it is a different model that happens to be readable, and presenting the second as an explanation is worse than presenting nothing. Fidelity is also local. A surrogate can agree well on average and disagree exactly on the cases you care about.

By that standard, the readable depth-3 tree printed above does not qualify: it accounts for about half the variation in the model's output. It is a reasonable description of the broad shape and it should not be presented as what the model does. Depth 6 reaches roughly 0.72 and is already too large to read aloud, which is the trade in a single sentence, readable and faithful are in tension, and a surrogate that is comfortable to present is usually the unfaithful one.

Explanations are not causal

This is the most consequential misunderstanding in the whole subject, so it is worth demonstrating rather than asserting. Add a column that is pure consequence. A retention call that the company makes because a customer looks like they are leaving, and watch every interpretability tool point at it.

import numpy as np
import pandas as pd

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

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

d = df.copy()
d['is_monthly'] = (d['contract'] == 'Month-To-Month').astype(int)
d['is_fibre'] = (d['internet_service'] == 'Fibre optic').astype(int)
FEATS = ['tenure_months', 'support_calls', 'monthly_charges',
         'total_charges', 'is_monthly', 'is_fibre']

X_tr, X_te, y_tr, y_te = train_test_split(
    d[FEATS], d['churned'], test_size=0.25, stratify=d['churned'],
    random_state=42)

# Impute after the split, with medians learned from the training rows
# only -- week 4's rule still applies. A Pipeline is the right way to do
# this; these pages need a plain frame to inspect column by column, so
# the two lines are written out instead.
medians = X_tr.median()
X_tr = X_tr.fillna(medians)
X_te = X_te.fillna(medians)
import numpy as np
import shap
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.metrics import roc_auc_score

rng = np.random.default_rng(0)
# The retention team rings customers who are about to leave. The call is
# caused by the churn risk; it does not cause it.
call_tr = (y_tr.to_numpy() * (rng.random(len(y_tr)) < 0.8)).astype(int)
call_te = (y_te.to_numpy() * (rng.random(len(y_te)) < 0.8)).astype(int)

A = X_tr.assign(retention_call=call_tr)
B = X_te.assign(retention_call=call_te)
m = HistGradientBoostingClassifier(random_state=42).fit(A, y_tr)

print('test AUC with the new column %.4f'
      % roc_auc_score(y_te, m.predict_proba(B)[:, 1]))
sv = shap.TreeExplainer(m)(B)
rank = np.abs(sv.values).mean(axis=0)
print('\nmean |shap|:')
for j in np.argsort(-rank):
    print('  %-18s %.4f' % (A.columns[j], rank[j]))
print('\nSHAP is right: the model does rely on it.')
print('Acting on it -- stop making retention calls -- changes nothing.')

test AUC with the new column 0.9489

mean |shap|:
  retention_call     4.6405
  tenure_months      0.9912
  is_monthly         0.6624
  monthly_charges    0.5483
  total_charges      0.5097
  support_calls      0.3393
  is_fibre           0.0401

SHAP is right: the model does rely on it.
Acting on it -- stop making retention calls -- changes nothing.

What every interpretability tool is actually telling you

SHAP, permutation importance and partial dependence all answer the same kind of question: how does this model use this feature? None of them answers what happens in the world if we change this feature? The second question needs an experiment, or causal assumptions you are prepared to defend in writing. Every time an importance ranking is read as a list of things to go and change, this distinction is the one being skipped.

Day 6 takeaway

Counterfactuals over actionable features turn a score into a decision somebody can act on. Surrogates are readable approximations and are worthless without a quoted fidelity. And no explanation method is causal. They describe the model, not the world.

Week 14 · Day 7 of 7

An Interpretability Report

One function, and the sentences to use when presenting it

By Pawan Rama Mali731 words

One function that produces everything a stakeholder needs, and an honest account of what to say about it.

The report

import numpy as np
import pandas as pd

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

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

d = df.copy()
d['is_monthly'] = (d['contract'] == 'Month-To-Month').astype(int)
d['is_fibre'] = (d['internet_service'] == 'Fibre optic').astype(int)
FEATS = ['tenure_months', 'support_calls', 'monthly_charges',
         'total_charges', 'is_monthly', 'is_fibre']

X_tr, X_te, y_tr, y_te = train_test_split(
    d[FEATS], d['churned'], test_size=0.25, stratify=d['churned'],
    random_state=42)

# Impute after the split, with medians learned from the training rows
# only -- week 4's rule still applies. A Pipeline is the right way to do
# this; these pages need a plain frame to inspect column by column, so
# the two lines are written out instead.
medians = X_tr.median()
X_tr = X_tr.fillna(medians)
X_te = X_te.fillna(medians)
import numpy as np
import shap
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.inspection import partial_dependence
from sklearn.metrics import roc_auc_score

def explain(model, X, y, top=4):
    print('=' * 58)
    print('MODEL EXPLANATION REPORT')
    print('=' * 58)
    p = model.predict_proba(X)[:, 1]
    print('rows %d   AUC %.4f   mean predicted risk %.3f'
          % (len(X), roc_auc_score(y, p), p.mean()))

    sv = shap.TreeExplainer(model)(X)
    rank = np.abs(sv.values).mean(axis=0)
    order = np.argsort(-rank)[:top]

    print('\nWHAT THE MODEL USES')
    for j in order:
        direction = np.corrcoef(X.iloc[:, j], sv.values[:, j])[0, 1]
        arrow = 'higher raises risk' if direction > 0 else 'higher lowers risk'
        print('  %-18s %.4f   %s' % (X.columns[j], rank[j], arrow))

    print('\nSHAPE OF THE TOP EFFECT')
    j = int(order[0])
    res = partial_dependence(model, X, [X.columns[j]], grid_resolution=6,
                             kind='average')
    for v, a in zip(res['grid_values'][0], res['average'][0]):
        print('  %-14s %8.2f -> %.4f' % (X.columns[j], v, a))

    print('\nWHO IS MOST AT RISK')
    for i in np.argsort(-p)[:3]:
        top_j = int(np.argmax(np.abs(sv.values[i])))
        print('  row %-5d p=%.3f  driven by %s=%.1f'
              % (i, p[i], X.columns[top_j], X.iloc[i, top_j]))
    print('=' * 58)

m = HistGradientBoostingClassifier(random_state=42).fit(X_tr, y_tr)
explain(m, X_te, y_te)

==========================================================
MODEL EXPLANATION REPORT
==========================================================
rows 750   AUC 0.7688   mean predicted risk 0.249

WHAT THE MODEL USES
  is_monthly         1.3276   higher raises risk
  tenure_months      0.9738   higher lowers risk
  monthly_charges    0.5814   higher raises risk
  total_charges      0.5090   higher lowers risk

SHAPE OF THE TOP EFFECT
  is_monthly         0.00 -> -3.0695
  is_monthly         1.00 -> 0.0764

WHO IS MOST AT RISK
  row 60    p=0.923  driven by monthly_charges=89.6
  row 545   p=0.923  driven by tenure_months=2.0
  row 119   p=0.915  driven by monthly_charges=92.5
==========================================================

Saying it out loud

Do not say	Say instead	Why
"Tenure is the most important factor"	"Tenure moves this model's predictions most"	Importance is a property of the model, not the world
"Reducing price would cut churn by 8%"	"Customers at lower prices churn less; testing a price change would tell us whether lowering it helps"	The first claims a causal effect the model cannot support
"The model is 87% accurate"	"It ranks well, AUC 0.87, and at our threshold it catches 6 in 10 leavers with 1 in 3 flagged"	Accuracy alone hides the trade-off that actually gets chosen
"The model explains why they left"	"These features drove this score"	It explains a prediction, not a departure

Before you present a model

Report a score from data the tuning never saw, not best_score_.
Show which features the model uses, and by how much, with the caveat that correlated features share credit ambiguously.
Show the shape of the top two or three effects, not just their rank, and check the correlations before you trust the curve.
Explain two or three individual predictions end to end, including one the model got wrong.
State what the model would need in order to be causal, and that it is not.
Say what the model does when it has never seen this kind of row before.
Agree the decision threshold with whoever bears the cost of the errors, using week 7's framing.
Write down the date, the data window and the version. Week 15 explains why that last one is not administrative pedantry.

Day 7 takeaway

Tuning and interpretability are the same discipline seen from two sides: one decides what the model is, the other decides whether anybody should believe it. Report held-out numbers, show the shape of the effects and not only their ranking, explain individual cases, and be precise about the line between what the model uses and what causes anything.