31 - The ML project process: splits, data leakage and metrics that lie

Imagine a student preparing for an exam with the papers of past years. If the real exam contains exactly those papers, the student will get top marks without having understood anything. A machine learning model does the same thing: if you evaluate it on the same data you trained it on, it will always tell you it is brilliant. The whole method of an ML project, from the split of the data to the choice of metrics, serves a single purpose: to stop you from telling yourself a lie about the quality of the model.

This lesson opens the module on classical machine learning. There is still no sophisticated algorithm: there is the process you will use with all the algorithms of the following lessons, and the three traps that ruin most real projects: data leakage, metrics that lie on imbalanced classes and overfitting hidden by a good accuracy.

What machine learning really is

A traditional program is a function written by you: rules in, answers out. In machine learning you reverse the relationship: you supply examples (data and expected answers) and an algorithm derives the function. The result of training is an object that takes a new input and returns a prediction. Two families are enough for now:

• Supervised learning: every example has a label (the correct answer). If the label is a category ("spam" / "not spam") we speak of classification; if it is a number ("price of a house") of regression.
• Unsupervised learning: there are no labels, you look for structure in the data (groups, directions of maximum variation). You will see it with clustering and dimensionality reduction.

Formally, in the supervised case you have pairs (x, y) drawn from an unknown distribution and you look for a function f such that f(x) ≈ y even on pairs never seen. The last part of the sentence is the heart of everything: the goal is not to remember the training data, it is to generalise.

The basic vocabulary, defined once and for all

• Example (or sample, or row): a single observation. In a dataset of customers, one customer.
• Feature: an input value that describes the example (age, amount of the last purchase). The set of the features of all the examples is a table with one row per example and one column per feature, the matrix X.
• Label (or target): the answer we want to predict, the vector y.
• Model: the function, with its parameters, that turns features into a prediction. Training (fit) means choosing the parameters by looking at the data; inference (predict) means applying the trained model to new data.
• Hyperparameter: a choice you make before training and that the algorithm does not learn (for example the maximum depth of a tree). Parameters are learned, hyperparameters are chosen with the validation set.

The five entries split into two groups: the data (example, feature, label) and what transforms them (model, hyperparameter). The distinction that really matters is the last one, between parameters that the model learns and hyperparameters that you choose by looking at the validation set: if you choose a hyperparameter by looking at the test set, the test set has already stopped being an exam.

The process, in six steps

1. Formulate the problem: what do I predict, with what cost of error? A false alarm and a missed diagnosis do not cost the same.
2. Collect and inspect the data: distribution of the classes, missing values, duplicates.
3. Split the data into training, validation and test, before any transformation.
4. Prepare the features (scaling, encoding of categories) learning the parameters only from the training set.
5. Train and compare models using the validation set.
6. Evaluate once on the test set and only at the end, with the metric consistent with the cost of the error.

The order of the six steps is not a preference: each depends on the previous one, and swapping two has a precise consequence. Splitting before preparing avoids leakage, evaluating on the test set only once avoids optimism, the initial formulation decides which metric will make sense. The step never to skip is the first, because it fixes the cost of the error.

The drawing shows the typical split. The percentages are not sacred: with a million examples a few percentage points are enough for validation and test, with a thousand examples cross-validation is preferable (further on).

Step 3 comes before step 4 and this order is not negotiable: it is the source of data leakage, which we will see shortly.

Train, validation and test: three different roles

The analogy of the student holds up well. The training set is the textbook with the worked exercises: you study on it. The validation set is the mock exam: you use it to understand what to fix (which algorithm, which hyperparameters). The test set is the real exam: you open it only once. If after seeing the mock exam mark you go back to studying "aiming" at those questions, the mock exam is no longer independent. The same goes for the test set: every time you look at it and then change the model, you use up some of its honesty.

With imbalanced classes the split must be stratified, that is, keep the same proportion of classes in every part. Let us generate a problem where the positives are about 5%:

Note on the code: the Python blocks of this article must be run in order, in a single script or notebook. Each one reuses the imports and variables defined in the previous blocks.

import numpy as np
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split

X, y = make_classification(
    n_samples=4000, n_features=10, n_informative=4,
    weights=[0.95, 0.05], flip_y=0.02, class_sep=1.0, random_state=0)

X_tr, X_te, y_tr, y_te = train_test_split(
    X, y, test_size=0.25, stratify=y, random_state=0)

print("train shape:", X_tr.shape, " test shape:", X_te.shape)
print("share of positives total/train/test: %.3f %.3f %.3f" % (y.mean(), y_tr.mean(), y_te.mean()))

train shape: (3000, 10)  test shape: (1000, 10)
share of positives total/train/test: 0.059 0.059 0.059

What the code does. make_classification generates a synthetic dataset: 4,000 examples, 10 features of which 4 are informative, with about 95% negative examples (weights) and 2% of the labels swapped at random (flip_y, to simulate annotation errors). train_test_split sets aside 25% of the data for the test; stratify=y guarantees the same share of positives in the two parts, and random_state makes the split reproducible. The three prints confirm the shapes (3,000 training rows, 1,000 test rows) and that the share of positives is 0.059 everywhere: about one example in 17 is positive.

Without stratify=y, with few positives you could end up with a test set almost devoid of the class you care about, and metrics that are unstable from one run to the next.

Metrics that lie: accuracy on imbalanced classes

Accuracy is the fraction of correct predictions: (right predictions) / (all predictions). It seems innocent. Now build the "laziest possible model": it always says "negative". With the positives at 5.9% of the data as here, that model has 94.1% accuracy and finds no positive. Let us check.

from sklearn.dummy import DummyClassifier
from sklearn.metrics import accuracy_score, recall_score

pigro = DummyClassifier(strategy="most_frequent").fit(X_tr, y_tr)
pred = pigro.predict(X_te)
print("accuracy of the lazy model: %.3f" % accuracy_score(y_te, pred))
print("recall on the positives:    %.3f" % recall_score(y_te, pred))

accuracy of the lazy model: 0.941
recall on the positives:    0.000

The code uses a DummyClassifier with the "most_frequent" strategy: it completely ignores the features and always predicts the most common class in the training set. It is the cheapest baseline, and it should be the first comparison of every project: a model that does not beat the lazy baseline has learned nothing. Here the accuracy is 0.941 and the recall is 0.000: no positive found.

To see what happens you need the confusion matrix, which crosses reality and prediction. Let us call "positive" the rare class we want to find:

                    predicted positive    predicted negative
actual positive           TP                    FN
actual negative           FP                    TN

TP are the true positives, FN the false negatives (positives that slipped past you), FP the false positives (false alarms), TN the true negatives. From here come the metrics that matter:

accuracy  = (TP + TN) / (TP + TN + FP + FN)
precision = TP / (TP + FP)     "when I say positive, how often am I right?"
recall    = TP / (TP + FN)     "how many of the true positives did I find?"
F1        = 2 · precision · recall / (precision + recall)

Precision is the seriousness of your alarms, recall is the coverage. F1 is the harmonic mean: the choice of the harmonic mean is not a whim, it serves to punish imbalances. With precision 1.0 and recall 0.01 the arithmetic mean would give 0.505, the harmonic one about 0.02, which describes a useless model better. Let us compute everything by hand from the matrix and compare with scikit-learn, using a real model (logistic regression, which you will see in detail in the next module):

from sklearn.linear_model import LogisticRegression
from sklearn.metrics import confusion_matrix, precision_score, f1_score

modello = LogisticRegression(max_iter=1000).fit(X_tr, y_tr)
pred = modello.predict(X_te)

vn, fp, fn, vp = confusion_matrix(y_te, pred).ravel()
print("TP=%d FP=%d FN=%d TN=%d" % (vp, fp, fn, vn))

acc  = (vp + vn) / (vp + vn + fp + fn)
prec = vp / (vp + fp) if (vp + fp) else 0.0
rec  = vp / (vp + fn)
f1   = 2 * prec * rec / (prec + rec)
print("by hand   -> acc %.3f prec %.3f rec %.3f f1 %.3f" % (acc, prec, rec, f1))
print("sklearn   -> acc %.3f prec %.3f rec %.3f f1 %.3f" % (
    accuracy_score(y_te, pred), precision_score(y_te, pred),
    recall_score(y_te, pred), f1_score(y_te, pred)))

TP=8 FP=1 FN=51 TN=940
by hand   -> acc 0.948 prec 0.889 rec 0.136 f1 0.235
sklearn   -> acc 0.948 prec 0.889 rec 0.136 f1 0.235

What the code does: it trains a logistic regression (a linear classifier that produces a probability; you will study it in full later in the module), then derives the four cells of the matrix with confusion_matrix(...).ravel() and computes the metrics with the formulas, line by line. The condition if (vp + fp) else 0.0 avoids a division by zero when the model never predicts "positive". The two prints coincide: the formulas by hand and those of scikit-learn are the same thing.

Look at the result: the "real" model has 94.8% accuracy against the 94.1% of the lazy one. It seems a minimal improvement, but the recall tells the truth: it finds 8 positives out of 59 (13.6%). The precision is high (8 alarms out of 9 are right) because the model exposes itself only in the obvious cases. An accuracy of 94.8% hides the fact that 51 positives out of 59 went unnoticed.

Note the order of confusion_matrix(...).ravel(): scikit-learn returns [[TN, FP], [FN, TP]], that is with class 0 first. Swapping that order is a very common mistake, and produces "inverted" precision and recall without any error message.

Decision threshold: the slider between precision and recall

A probabilistic classifier produces a probability; "positive" means probability ≥ threshold, and the default threshold 0.5 is a convention, not a law. Lowering it finds more positives (recall up) but generates more false alarms (precision down). Which compromise to choose depends on the cost of the errors, not on statistics: in a medical screening you accept some false alarms rather than miss a sick patient; in a filter that blocks accounts you do the opposite.

In the experiment below, going down from threshold 0.5 to 0.1 brings the recall from 0.136 to 0.390 but makes the precision collapse from 0.889 to 0.261. There is no "right" value in absolute terms: there is the one that is right for your cost of errors.

p = modello.predict_proba(X_te)[:, 1]
for soglia in (0.5, 0.3, 0.2, 0.1):
    pr = (p >= soglia).astype(int)
    print("threshold %.1f  precision %.3f  recall %.3f" % (
        soglia, precision_score(y_te, pr, zero_division=0), recall_score(y_te, pr)))

threshold 0.5  precision 0.889  recall 0.136
threshold 0.3  precision 0.800  recall 0.271
threshold 0.2  precision 0.667  recall 0.305
threshold 0.1  precision 0.261  recall 0.390

Data leakage: the trap that makes a useless model look brilliant

Data leakage happens when, directly or indirectly, information enters the training process that in the real world you would not have at the moment of the prediction. The typical forms:

• Preprocessing on the whole dataset: you compute the mean and standard deviation (or select the features) on all the data, then split. The test set has already influenced the training.
• Features that contain the answer: a column "days since cancellation" to predict the cancellation.
• Duplicates or near-duplicates that ended up both in train and in test.
• Mixed-up temporal data: predicting the future using data that come later.

First, two terms. Feature selection: choosing, among many columns, the few that seem most useful (here, those with the highest F statistic with respect to the label). Cross-validation: repeating training and evaluation on different portions of the data; we explain it better below, for now it is enough to know that it returns an estimate of accuracy on "unseen" data.

To make the point evident we build an extreme but instructive case: 100 samples, 5,000 features, all pure noise, and a random label. No honest model can do better than chance (50%). Now we commit the classic mistake: choosing the 20 features "most correlated with the label" by looking at all the data, and only afterwards doing the cross-validation.

from sklearn.feature_selection import SelectKBest, f_classif
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import make_pipeline

rng = np.random.default_rng(0)
Xr = rng.normal(size=(100, 5000))
yr = rng.integers(0, 2, size=100)

# WRONG: the selection has already seen the labels of every fold
Xsel = SelectKBest(f_classif, k=20).fit_transform(Xr, yr)
sbagliato = cross_val_score(LogisticRegression(max_iter=1000), Xsel, yr, cv=5).mean()

# RIGHT: the selection sits INSIDE the pipeline, redone in every fold
pipe = make_pipeline(SelectKBest(f_classif, k=20), LogisticRegression(max_iter=1000))
giusto = cross_val_score(pipe, Xr, yr, cv=5).mean()

print("CV accuracy with leakage:    %.2f" % sbagliato)
print("CV accuracy without leakage: %.2f" % giusto)

CV accuracy with leakage:    0.85
CV accuracy without leakage: 0.48

What the code does. We generate 100 rows and 5,000 columns of random numbers and a random label: by construction no column has any relationship with the label. In the WRONG case SelectKBest is fitted once on all 100 samples and returns the 20 "best" columns; among 5,000 columns of noise, by pure chance, twenty correlate with the label, and the subsequent cross-validation finds them "predictive", because the labels of the test samples have already entered the choice. In the RIGHT case make_pipeline joins selection and classifier in a single object: at each fold the selection is redone looking only at the training samples of that fold.

The first number is a mirage (0.85 on data that are noise): on data that are literally noise, the model "discovers" a nonexistent structure because the feature selection used the labels of the samples that would later have been test samples. The second number (0.48, that is chance) is the truth. The practical rule that follows is just one: any step that "learns" something from the data (scaler, imputation, selection, PCA, target encoding) goes inside a Pipeline, so that it is redone on every fold using only the training part of that fold.

Cross-validation: evaluating without wasting data

With little data, setting aside 25% for a single evaluation is costly and noisy. k-fold cross-validation divides the training set into k parts; k times it uses k−1 to train and one to validate, then averages the results. You get a more stable estimate and also a measure of its variability (the standard deviation across the folds). The price is computing time: k trainings instead of one. For imbalanced classes use StratifiedKFold, which cross_val_score already applies by default in classification problems.

Overfitting masked by good accuracy

Overfitting is when the model memorises the noise of the training set instead of the signal. The symptom is a gap between performance on the training set and on unseen data. A decision tree with no depth limit shows it very well:

from sklearn.tree import DecisionTreeClassifier

albero = DecisionTreeClassifier(random_state=0).fit(X_tr, y_tr)
print("train: acc %.3f  recall %.3f" % (albero.score(X_tr, y_tr), recall_score(y_tr, albero.predict(X_tr))))
print("test:  acc %.3f  recall %.3f" % (albero.score(X_te, y_te), recall_score(y_te, albero.predict(X_te))))

train: acc 1.000  recall 1.000
test:  acc 0.945  recall 0.492

On the training set the tree is perfect (accuracy and recall at 1.000); on the test set the recall on the positives collapses to 0.492, that is it finds only about half of the positives. But the overall accuracy on the test set stays at 0.945, a high and reassuring number because it is dominated by the negatives: it is overfitting in disguise. Two defences: always compare the right metric between train and validation, and prefer simpler models until validation proves otherwise (Occam's razor in practical form). The systematic remedy, regularisation, is the topic of the lesson on overfitting and regularisation.

Trade-offs to keep in mind

• More test data = a more reliable estimate, but less data to train on.
• Precision versus recall: they are not maximised together, you choose according to the cost of the errors.
• Cross-validation versus hold-out: the former is more stable and slower; with huge datasets hold-out is enough.
• Simple model versus complex: the complex one almost always wins on the training set, and often loses outside it.

The four compromises have the same structure: every gain is paid for with a loss on another axis (training data, recall, computing time, ability to generalise). None is correct in absolute terms; they must be chosen according to the cost of the error and the amount of data, and declared in the result.

Checklist of the traps

• Did you split the data before scaling, imputing or selecting features?
• Is the split stratified (imbalanced classes) or temporal (time series)?
• Did you look at precision and recall, not only at accuracy? Do you have a comparison with the lazy model?
• Was the test set used only once, at the end?
• Are there duplicates between train and test? Does a feature already "know" the answer?
• Is the order of confusion_matrix what you think it is?

The checklist is ordered from the point where mistakes are most frequent to the least frequent: split before transformation, type of split, metric, use of the test set, duplicates and features that know the answer, order of the cells. If a check seems obvious to you, it is the one you are about to violate; it is worth running through it before every evaluation, not after.

Beyond the threshold: ROC curve, precision-recall curve and AUC

So far you have evaluated a model at a fixed threshold. But a probabilistic classifier does not have "one" threshold: it has infinitely many, and each produces a different pair of precision and recall. Two curves summarise the whole compromise. The ROC curve (Receiver Operating Characteristic) plots, as the threshold varies, the true positive rate TPR = TP/(TP + FN) (which is the recall) against the false positive rate FPR = FP/(FP + TN) (the fraction of negatives you mistake for positives). A useless model sits on the diagonal; a perfect one passes through the top-left corner. The AUC (area under the ROC curve) is a single number between 0.5 (chance) and 1 (perfect) and has an elegant interpretation: it is the probability that a randomly chosen positive receives a higher score than a randomly chosen negative. We verify it by computing the AUC in two independent ways, by hand.

from sklearn.metrics import roc_auc_score, average_precision_score, precision_recall_curve

scores = LogisticRegression(max_iter=1000).fit(X_tr, y_tr).predict_proba(X_te)[:, 1]

# 1) ROC curve by hand: one threshold per distinct score, from the highest to the lowest
order = np.argsort(-scores)
y_sorted = y_te[order]
tpr = np.concatenate([[0], np.cumsum(y_sorted) / y_sorted.sum()])
fpr = np.concatenate([[0], np.cumsum(1 - y_sorted) / (1 - y_sorted).sum()])
auc_trapezoids = np.trapezoid(tpr, fpr)                              # area with the trapezoidal rule

# 2) AUC as a probability: compare every positive with every negative (ties = 1/2)
pos, neg = scores[y_te == 1], scores[y_te == 0]
auc_pairs = ((pos[:, None] > neg[None, :]).sum() + 0.5 * (pos[:, None] == neg[None, :]).sum()) / (len(pos) * len(neg))

print("AUC with trapezoids: %.4f   AUC as a probability over pairs: %.4f   scikit-learn: %.4f" % (auc_trapezoids, auc_pairs, roc_auc_score(y_te, scores)))
print("positive-negative pairs compared: %d × %d = %d" % (len(pos), len(neg), len(pos) * len(neg)))

# comparison with the precision-recall curve
ap = average_precision_score(y_te, scores)
print("ROC AUC = %.3f  but  average precision (area under the precision-recall curve) = %.3f   (random model: %.3f)" % (roc_auc_score(y_te, scores), ap, y_te.mean()))
prec, rec, thr = precision_recall_curve(y_te, scores)
for target in (0.2, 0.4, 0.6):
    print("with recall ≥ %.1f the maximum precision reachable is %.3f" % (target, prec[:len(rec)][rec >= target].max()))

AUC with trapezoids: 0.5886   AUC as a probability over pairs: 0.5886   scikit-learn: 0.5886
positive-negative pairs compared: 59 × 941 = 55519
ROC AUC = 0.589  but  average precision (area under the precision-recall curve) = 0.363   (random model: 0.059)
with recall ≥ 0.2 the maximum precision reachable is 0.933
with recall ≥ 0.4 the maximum precision reachable is 0.180
with recall ≥ 0.6 the maximum precision reachable is 0.060

What the code does. First the scores of a logistic regression on the test set. Method 1 sorts the samples from the highest score: lowering the threshold progressively amounts to walking down this list and declaring the first k positive; at each step TPR and FPR are the cumulative fractions of positives and negatives met. The area under the curve is computed with the trapezoidal rule (np.trapezoid). Method 2 draws nothing: it compares all the 59 × 941 = 55,519 (positive, negative) pairs and counts in how many the positive has a higher score, with ties counted as half. The three numbers coincide (0.5886): the AUC is that probability.

Look at the rest of the output, because it tells two different stories about the same model. The ROC AUC is 0.589, just above chance: if you take a positive and a negative at random, the model orders them correctly in only 59% of the pairs. The average precision (area under the precision-recall curve) is 0.363, more than six times the reference of a random model (0.059, the share of positives). And the table says why: with recall of at least 0.2 the maximum precision reachable is 0.93, with recall 0.4 it drops to 0.18 and with recall 0.6 it is 0.06, that is the one of chance. The model is very sure about roughly 20% of the positives (very high scores, almost all true positives) and for the rest it does not tell them apart from the negatives. The AUC weighs all the pairs equally, including the very many in which the positive is one of the "difficult" ones; the average precision looks mostly at the head of the ranking.

Which metric is right depends on the use: if an analyst can review only the first 100 alarms, the head of the ranking counts and the average precision is far more informative; if you need to order the whole population well, the AUC counts. A second reason for caution: with very imbalanced classes the FPR has all the negatives (941) in the denominator, so even dozens of false alarms produce a tiny FPR and the ROC curve looks more optimistic than an operator, who sees only the alarms, experiences.

k-fold cross-validation, written by hand

You have used cross_val_score as a black box. It is worth writing it once, because validation bugs are among the most insidious. The idea: divide the examples into k groups (folds) of similar size; for each group, train on the other k − 1 and evaluate on the one left out; then report the mean and standard deviation of the k scores. In the stratified case, each group keeps the proportion of the classes. The crucial point, repeated: everything that learns from the data (the scaler, the selection, the model) must be trained inside the loop, on the training part of the fold only.

from sklearn.model_selection import StratifiedKFold, cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

def stratified_folds(y, k, seed=0):
    rng = np.random.default_rng(seed)
    fold_of = np.empty(len(y), dtype=int)
    for c in np.unique(y):
        idx = rng.permutation(np.where(y == c)[0])                    # shuffle the examples of each class
        fold_of[idx] = np.arange(len(idx)) % k                        # and deal them out in turn to the k groups
    return fold_of

folds = stratified_folds(y_tr, k=5)
print("examples per group: %s   positives per group: %s" % (np.bincount(folds), [int(y_tr[folds == f].sum()) for f in range(5)]))

results = []
for f in range(5):
    train_idx, val_idx = folds != f, folds == f
    pipe = make_pipeline(StandardScaler(), LogisticRegression(max_iter=1000)).fit(X_tr[train_idx], y_tr[train_idx])    # the scaler sees only the training part of the fold
    p_val = pipe.predict_proba(X_tr[val_idx])[:, 1]
    results.append(average_precision_score(y_tr[val_idx], p_val))
print("by hand:        average precision per fold %s   mean %.3f ± %.3f" % (np.round(results, 3), np.mean(results), np.std(results)))

reference = cross_val_score(make_pipeline(StandardScaler(), LogisticRegression(max_iter=1000)), X_tr, y_tr,
                            cv=StratifiedKFold(5, shuffle=True, random_state=0), scoring="average_precision")
print("scikit-learn:   average precision per fold %s   mean %.3f ± %.3f" % (np.round(reference, 3), reference.mean(), reference.std()))

examples per group: [601 600 600 600 599]   positives per group: [36, 35, 35, 35, 35]
by hand:        average precision per fold [0.497 0.572 0.389 0.438 0.453]   mean 0.470 ± 0.062
scikit-learn:   average precision per fold [0.502 0.51  0.414 0.557 0.363]   mean 0.469 ± 0.070

What the code does. stratified_folds assigns each example a group number from 0 to 4: for each class it shuffles the indices and deals them out "in turn", so that each group receives 35 or 36 positives (177 in all) and 599-601 examples. The loop trains a Pipeline (scaler and model) on the four training groups and evaluates the average precision on the group left out. The result by hand (0.470 ± 0.062) and that of scikit-learn (0.469 ± 0.070) are not identical fold by fold, because the two methods shuffle differently and so the groups do not coincide, but they are statistically equivalent: same mean within the uncertainty. The deviation between groups is the valuable datum: the five values range from 0.36 to 0.56, a wide interval, because the positives per group are only 35 and it takes just a few "difficult" ones to move the score. If the improvement between two models is smaller than that deviation, you have no evidence that it is real: cross-validation serves to measure not only the mean but also how far to trust it.

From metric to cost: choosing the threshold with money, not with 0.5

In real life errors have different costs, and the right threshold is the one that minimises the expected cost, not the one that maximises an abstract metric. Imagine a system that flags suspicious transactions: a false negative (an undetected fraud) costs on average 50 units, a false positive (a useless block of an honest customer, with the support call) costs 1. The numbers are invented for the exercise: what matters is the method. The correct procedure is to choose the threshold on the training set, using out-of-sample predictions (out-of-fold, that is obtained with cross-validation), and only then to measure the cost on the test set once.

from sklearn.model_selection import cross_val_predict

COST_FN, COST_FP = 50.0, 1.0                                          # illustrative costs, in arbitrary units
cv_scores = cross_val_predict(make_pipeline(StandardScaler(), LogisticRegression(max_iter=1000)), X_tr, y_tr,
                              cv=StratifiedKFold(5, shuffle=True, random_state=0), method="predict_proba")[:, 1]

def total_cost(y, p, threshold):
    pred = p >= threshold
    return COST_FN * ((~pred) & (y == 1)).sum() + COST_FP * (pred & (y == 0)).sum()

thresholds = np.logspace(-3, np.log10(0.99), 150)                  # from 0.001 to 0.99, on a logarithmic scale
costs_cv = np.array([total_cost(y_tr, cv_scores, t) for t in thresholds])
best_t = thresholds[costs_cv.argmin()]

final = make_pipeline(StandardScaler(), LogisticRegression(max_iter=1000)).fit(X_tr, y_tr)
p_test = final.predict_proba(X_te)[:, 1]
print("threshold chosen on the training set with cross-validation: %.3f   (mean cost per example: %.3f)" % (best_t, costs_cv.min() / len(y_tr)))
print("%-30s %-12s %-12s %-12s" % ("", "threshold", "recall", "cost/example"))
theory_t = COST_FP / (COST_FP + COST_FN)                          # theoretical optimal threshold for calibrated probabilities
for name, t in (("threshold 0.5 (default)", 0.5), ("threshold chosen by cost (CV)", best_t), ("theoretical C_FP/(C_FP+C_FN)", theory_t), ("always flag (threshold 0)", 0.0)):
    pred = p_test >= t
    print("%-30s %-12.3f %-12.3f %-12.3f" % (name, t, recall_score(y_te, pred), total_cost(y_te, p_test, t) / len(y_te)))
print("cost of a model that never flags anything: %.3f per example" % (COST_FN * y_te.sum() / len(y_te)))

threshold chosen on the training set with cross-validation: 0.003   (mean cost per example: 0.938)
                               threshold    recall       cost/example
threshold 0.5 (default)        0.500        0.119        2.601       
threshold chosen by cost (CV)  0.003        0.966        1.038       
theoretical C_FP/(C_FP+C_FN)   0.020        0.746        1.504       
always flag (threshold 0)      0.000        1.000        0.941       
cost of a model that never flags anything: 2.950 per example

What the code does. cross_val_predict returns, for each example in the training set, the probability predicted by a model that has not seen that example: honest predictions on which one can choose. The function total_cost sums the costs of the two types of error. We scan 150 thresholds on a logarithmic scale (the interesting probabilities are near zero) and pick the one with minimum cost on the training set: 0.003. Only then is the test set opened, and the result is instructive on three levels. First, the default threshold 0.5 finds only 12% of the positives and costs 2.60 per example, almost as much as doing nothing (2.95). Second, the threshold chosen by cross-validation (0.003) finds 97% of the positives and costs 1.04: a great improvement, but the model is in effect flagging almost everything.

Third, and most important: the trivial strategy "always flag" costs 0.941 per example, less than any threshold of the model. The model, with AUC 0.59, adds no value over that baseline, because with a cost ratio of 50 to 1 and 6% positives, calling everything positive is already worthwhile. The result is not a flaw of the exercise: it is exactly what a cost-based evaluation must discover.

There is also a formula that predicts the threshold. For calibrated probabilities, flagging is worthwhile if (1 − p)·C_FP < p·C_FN, that is if p > C_FP/(C_FP + C_FN) = 1/51 ≈ 0.020. The theoretical threshold 0.020 gives a cost of 1.50: worse than 0.003 because the model's probabilities are not well calibrated (lesson on logistic regression) and the validation corrected them empirically. Moral: before evaluating a model, compare it with the trivial strategies (do nothing, always act), and choose the threshold with the costs of your problem.

Three exercises

Exercise 1. Change the costs (FN = 10, then FN = 200) and observe how the optimal threshold moves. What happens to the recall? Exercise 2. Repeat the ROC versus precision-recall comparison with a non-imbalanced dataset (weights=[0.5, 0.5] in make_classification): does the difference between ROC AUC and average precision shrink? Exercise 3. Deliberately insert leakage into the k-fold loop (standardise before splitting) and measure how much the average precision changes: it is a good estimate of the "optimism" of the leakage for this problem.

The code exercises of this lesson are in a single file, esercizi.py (about 16 KB). Three code exercises on a new domain, the inspection of welds: the metrics derived from the four cells, the minimum-cost threshold found by scanning, and a split by group that exposes the leakage from duplicate rows. It requires Python 3.12+ and NumPy 2.3 - no scikit-learn. Run it with python3 esercizi.py: it also runs with the parts to complete left empty and tells you which ones are missing.

Previous lesson in the course: Markov Chains: The Property That Forgets the Past, and the Distribution That Stops Changing. Other lessons in the series that link to this one: Linear Algebra 2: Eigenvalues, SVD and PCA, Support Vector Machine: the maximum margin and the kernel trick, KNN: classifying by looking at neighbours (and the curse of dimensionality), Random Forest: why a hundred trees err less than one, Reinforcement Learning: The Agent That Learns by Trial and Error (and the Exploration-Exploitation Dilemma).

In the next lessons you will apply this same process to every algorithm, starting from the simplest to understand: Naive Bayes.