11 - Metrics and ML Model Evaluation

A model that always answers «no fraud» has an accuracy of 97.4% on a set of transactions in which 2.6% are fraud, and never uncovers a single fraud. The number you looked at was right and told you nothing. Evaluating a model means choosing which question you ask it: how many mistakes it makes overall, how often it is right when it says «positive», how many positives it manages to find, how reliable the probabilities it states are, how far off it is in euros or in degrees.

Every question has its own metric and each can mislead in its own way. This article, the eleventh in the ML Fundamentals series, computes the main metrics by hand and with scikit-learn, shows with an experiment where each one lies, and tells you how far to trust a single number.

Metrics are used on a test set that the model has not seen, as you did with cross-validation and with the models of the series, from logistic regression onward. The data are synthetic with a fixed seed, declared case by case, and the cost figures are hypothetical and said to be so: no real datasets are downloaded. The principles on the process (split, leakage) are in the Machine Learning Engineer series article on split, data leakage and metrics that lie.

The confusion matrix and why accuracy deceives

The confusion matrix counts the four outcomes of a binary classifier: true positives (TP, you said yes and it was yes), false positives (FP, you said yes and it was no), false negatives (FN, you said no and it was yes) and true negatives (TN). All the classification metrics are ratios between these four numbers.

Accuracy is $(TP + TN)/N$; precision is $TP/(TP + FP)$, that is the fraction of true alarms; recall is $TP/(TP + FN)$, the fraction of positives found; F1 is their harmonic mean $2PR/(P + R)$, which is low if either of the two is low. Balanced accuracy is the mean of the recall on each class.

The listing generates 20,000 examples with 10 features (about 2.6% positives), keeps half for the test and compares two models: a classifier that always answers «no» (DummyClassifier with the most_frequent strategy) and a logistic regression. Then it recomputes the four metrics from the confusion matrix of the logistic regression, with the formulas just written, to check that the numbers add up.

# Imbalanced classes (about 2.6% positives): whoever always answers "no" has 97.4% accuracy
import numpy as np
from sklearn.datasets import make_classification
from sklearn.dummy import DummyClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import (accuracy_score, balanced_accuracy_score, confusion_matrix,
                             f1_score, precision_score, recall_score)
from sklearn.model_selection import train_test_split

X, y = make_classification(n_samples=20000, n_features=10, n_informative=5, weights=[0.98, 0.02],
                           flip_y=0.01, class_sep=0.8, random_state=0)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.5, random_state=0, stratify=y)
print("positives in the test set:", int(yte.sum()), "out of", len(yte))

models = {"always 'no'": DummyClassifier(strategy="most_frequent"),
          "logistic regression": LogisticRegression(max_iter=1000)}
for name, m in models.items():
    p = m.fit(Xtr, ytr).predict(Xte)
    tn, fp, fn, tp = confusion_matrix(yte, p).ravel()
    print(f"\n{name}: TN={tn} FP={fp} FN={fn} TP={tp}")
    print(f"  accuracy {accuracy_score(yte, p):.4f}  balanced accuracy {balanced_accuracy_score(yte, p):.4f}")
    print(f"  precision {precision_score(yte, p, zero_division=0):.4f}  recall {recall_score(yte, p):.4f}  F1 {f1_score(yte, p):.4f}")

# The same metrics recomputed by hand from the confusion matrix of the logistic regression
precision = tp / (tp + fp); recall = tp / (tp + fn)
f1 = 2 * precision * recall / (precision + recall)
acc = (tp + tn) / (tp + tn + fp + fn)
print("\nby hand: accuracy %.4f  precision %.4f  recall %.4f  F1 %.4f" % (acc, precision, recall, f1))
print("maximum gap from the sklearn values:", max(abs(acc - accuracy_score(yte, p)), abs(precision - precision_score(yte, p)),
      abs(recall - recall_score(yte, p)), abs(f1 - f1_score(yte, p))))

# Output:
# positives in the test set: 256 out of 10000
#
# always 'no': TN=9744 FP=0 FN=256 TP=0
#   accuracy 0.9744  balanced accuracy 0.5000
#   precision 0.0000  recall 0.0000  F1 0.0000
#
# logistic regression: TN=9742 FP=2 FN=179 TP=77
#   accuracy 0.9819  balanced accuracy 0.6503
#   precision 0.9747  recall 0.3008  F1 0.4597
#
# by hand: accuracy 0.9819  precision 0.9747  recall 0.3008  F1 0.4597
# maximum gap from the sklearn values: 5.551115123125783e-17

Out of 10,000 test examples the positives are 256. Whoever always answers «no» has accuracy 0.9744, but zero recall, zero F1 and balanced accuracy 0.5, the value of chance. The logistic regression has accuracy 0.9819 and seems only slightly better: just 0.75 points higher.

The precision, however, is 0.9747 (out of 79 alarms, 77 true) and the recall 0.3008: it finds 77 positives out of 256 and misses 179. This is the real question, and accuracy did not see it. The hand formulas coincide with scikit-learn (gap of 6e-17), so you can compute everything from a confusion matrix printed on a sheet of paper.

ROC, AUC and precision-recall: what changes with prevalence

A classifier with probabilities produces a score, and the decision depends on a threshold. The ROC curve shows, as the threshold varies, the recall (fraction of positives found) against the false positive rate (fraction of negatives mistaken for positives); the AUC is the area under the curve, and equals the probability that a randomly chosen positive has a higher score than a randomly chosen negative.

The precision-recall curve shows precision against recall, and its summary is the average precision (AP). Saito and Rehmsmeier (PLOS ONE, 2015) argue that the ROC can mislead on imbalanced data because of a wrong interpretation of specificity, while precision-recall evaluates the fraction of true positives among the positive predictions.

The listing verifies it without models: it generates scores for the positives from a normal with mean 1.5 and for the negatives from a normal with mean 0 (variance 1), and changes only the proportion of positives, from 50% to 0.1%, with 400,000 examples in all. The «classifier» is the same on every row: the separation between the two populations does not change.

# Same classifier, different prevalence: the AUC does not move, the precision does
import numpy as np
from sklearn.metrics import average_precision_score, precision_recall_curve, roc_auc_score

rng = np.random.default_rng(0)
N = 400_000
print(f"{'prevalence':>11s} {'positives':>9s} {'ROC-AUC':>8s} {'AP':>7s} {'random AP':>11s} {'precision at recall>=0.8':>26s}")
for prev in (0.5, 0.1, 0.01, 0.001):
    n_pos = int(N * prev)
    s = np.r_[rng.normal(1.5, 1.0, n_pos), rng.normal(0.0, 1.0, N - n_pos)]     # scores: positives and negatives
    y = np.r_[np.ones(n_pos), np.zeros(N - n_pos)]
    P, R, _ = precision_recall_curve(y, s)
    prec_80 = P[R >= 0.8].max()
    print(f"{prev:11.3f} {n_pos:9d} {roc_auc_score(y, s):8.3f} {average_precision_score(y, s):7.3f} {prev:11.3f} {prec_80:26.3f}")

# Output:
#  prevalence positives  ROC-AUC      AP   random AP   precision at recall>=0.8
#       0.500    200000    0.855   0.854       0.500                      0.757
#       0.100     40000    0.856   0.480       0.100                      0.259
#       0.010      4000    0.853   0.111       0.010                      0.031
#       0.001       400    0.863   0.017       0.001                      0.003

The AUC stays at 0.85-0.86 on all rows: it depends on the scores of the positives against those of the negatives, not on how many there are. The average precision collapses from 0.854 to 0.017, and the precision you get if you want to find 80% of the positives goes from 0.757 to 0.003: with one fraud in a thousand transactions, almost all the alarms of a system that found 80% of the frauds would be false.

The baseline of a random classification is the prevalence (the «random AP» column): 0.017 is 17 times chance, and must be read that way. So the AUC tells you how well the model ranks; the AP, or the precision at a fixed recall, tells you how useful it is when positives are rare. Report both, and report the prevalence.

The threshold, the cost and the F-beta

The 0.5 threshold is a convention, and is fine only if the two errors cost the same and the classes are balanced. A lower threshold finds more positives and produces more false alarms. If missing a positive costs much more than a false alarm, the right threshold is lower, and it is chosen by the cost, not by the model.

The F-beta generalizes F1 with a weight: $F_\beta = (1+\beta^2)PR/(\beta^2 P + R)$; with $\beta = 2$ the recall counts four times more than the precision in the denominator, with $\beta = 0.5$ the precision.

The listing uses the same data as before and splits the training set into two parts: one to train the logistic regression and a validation one to choose the threshold, which is then evaluated on the test set never touched. The cost is hypothetical, declared in the code: missing a positive costs 50 units, a false alarm 1. It compares four thresholds: 0.5, the one that maximizes F1 on the validation set, the one that maximizes F2, and the one that minimizes the cost.

# The 0.5 threshold is not sacred: choosing the threshold with a cost function, and F-beta computed by hand
import numpy as np
from sklearn.datasets import make_classification
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import confusion_matrix, f1_score, fbeta_score
from sklearn.model_selection import train_test_split

X, y = make_classification(n_samples=20000, n_features=10, n_informative=5, weights=[0.98, 0.02],
                           flip_y=0.01, class_sep=0.8, random_state=0)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.5, random_state=0, stratify=y)
Xfit, Xval, yfit, yval = train_test_split(Xtr, ytr, test_size=0.5, random_state=1, stratify=ytr)   # validation set for the threshold

m = LogisticRegression(max_iter=1000).fit(Xfit, yfit)
p_val, p_te = m.predict_proba(Xval)[:, 1], m.predict_proba(Xte)[:, 1]
C_FN, C_FP = 50.0, 1.0       # HYPOTHETICAL costs: missing a positive costs 50 times a false alarm

def cost(y, p, t):
    tn, fp, fn, tp = confusion_matrix(y, p >= t).ravel()
    return C_FN * fn + C_FP * fp

thresholds = np.linspace(0.005, 0.95, 190)
t_cost = thresholds[np.argmin([cost(yval, p_val, t) for t in thresholds])]
t_f1 = thresholds[np.argmax([f1_score(yval, p_val >= t) for t in thresholds])]
t_f2 = thresholds[np.argmax([fbeta_score(yval, p_val >= t, beta=2) for t in thresholds])]
print(f"thresholds chosen on validation: minimum cost -> {t_cost:.3f}, maximum F1 -> {t_f1:.3f}, maximum F2 -> {t_f2:.3f}")

print(f"\n{'threshold':>18s} {'TP':>4s} {'FP':>5s} {'FN':>4s} {'test cost':>11s} {'F1':>6s} {'F2':>6s}")
for name, t in [("0.5 (default)", 0.5), ("maximum F1", t_f1), ("maximum F2", t_f2), ("minimum cost", t_cost)]:
    tn, fp, fn, tp = confusion_matrix(yte, p_te >= t).ravel()
    print(f"{name:>18s} {tp:4d} {fp:5d} {fn:4d} {cost(yte, p_te, t):11.0f} {f1_score(yte, p_te >= t):6.3f} {fbeta_score(yte, p_te >= t, beta=2):6.3f}")

# F-beta by hand: (1 + b^2) P R / (b^2 P + R)
tn, fp, fn, tp = confusion_matrix(yte, p_te >= t_cost).ravel()
P, R = tp / (tp + fp), tp / (tp + fn)
for b in (0.5, 1, 2):
    print(f"F{b}: by hand {(1 + b * b) * P * R / (b * b * P + R):.4f}  sklearn {fbeta_score(yte, p_te >= t_cost, beta=b):.4f}")

# Output:
# thresholds chosen on validation: minimum cost -> 0.055, maximum F1 -> 0.225, maximum F2 -> 0.135
#
#          threshold   TP    FP   FN   test cost     F1     F2
#      0.5 (default)   83     2  173        8652  0.487  0.374
#         maximum F1  123    47  133        6697  0.577  0.515
#         maximum F2  136   142  120        6142  0.509  0.522
#       minimum cost  167   631   89        5081  0.317  0.458
# F0.5: by hand 0.2422  sklearn 0.2422
# F1: by hand 0.3169  sklearn 0.3169
# F2: by hand 0.4583  sklearn 0.4583

The threshold that minimizes the cost on validation is 0.055, that of F1 is 0.225 and that of F2 is 0.135. On the test set, with the 0.5 threshold the cost is 8,652 (2 false alarms and 173 positives lost); with the minimum-cost threshold it falls to 5,081, 41% less, with 167 positives found but 631 false alarms.

On that row F1 is the lowest in the table (0.317) and F2 is not the highest (0.458 against 0.522): each threshold is optimal for its own metric and worse for the others. The three F-betas computed by hand coincide with scikit-learn (0.2422, 0.3169, 0.4583). The question that decides is not «which metric is better» but «how much does each error cost». If you do not know, the metrics do not know it in your place.

Calibrated probabilities

A model is calibrated if, among the times it states 0.7, the positive happens about 70% of the time. It matters when you use the probability and not just the class (risk, price, priority). The Brier score is the mean squared error between the probability and the 0/1 outcome; the scikit-learn documentation warns that it mixes calibration and ability to discriminate, so a lower Brier does not necessarily mean a better calibrated model.

The ECE (Expected Calibration Error) divides the predictions into bins of equal width (Guo and colleagues use 15) and computes the mean, weighted by the number of predictions, of the absolute difference between the fraction of positives and the mean probability of each bin.

The listing compares three models on the same test set: logistic regression, a random forest of 100 trees, and the same forest with isotonic calibration (CalibratedClassifierCV, method="isotonic", 3 folds). Then it prints a five-bin reliability table for the uncalibrated forest, to see where it goes wrong.

# Calibration: if the model says 0.8, is it true 80% of the time? Brier score, ECE (15 bins) and reliability table
import numpy as np
from sklearn.calibration import CalibratedClassifierCV, calibration_curve
from sklearn.datasets import make_classification
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import brier_score_loss, log_loss
from sklearn.model_selection import train_test_split

X, y = make_classification(n_samples=12000, n_features=10, n_informative=5, flip_y=0.1, class_sep=1.0, random_state=0)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.4, random_state=0, stratify=y)

def ece(y, p, bins=15):
    """Expected Calibration Error: weighted mean of |fraction of positives - mean probability| per bin."""
    idx = np.minimum((p * bins).astype(int), bins - 1)
    tot = 0.0
    for b in range(bins):
        m = idx == b
        if m.any():
            tot += m.mean() * abs(y[m].mean() - p[m].mean())
    return tot

models = {
    "logistic regression": LogisticRegression(max_iter=1000),
    "random forest (100 trees)": RandomForestClassifier(n_estimators=100, random_state=0, n_jobs=1),
    "forest + isotonic calibration": CalibratedClassifierCV(RandomForestClassifier(n_estimators=100, random_state=0, n_jobs=1),
                                                             method="isotonic", cv=3),
}
print(f"{'model':34s} {'Brier':>7s} {'log-loss':>9s} {'ECE':>7s} {'accuracy':>12s}")
prob = {}
for name, m in models.items():
    p = m.fit(Xtr, ytr).predict_proba(Xte)[:, 1]
    prob[name] = p
    print(f"{name:34s} {brier_score_loss(yte, p):7.4f} {log_loss(yte, p):9.4f} {ece(yte, p):7.4f} {((p > 0.5) == yte).mean():12.4f}")

print("\nreliability table of the forest (uncalibrated): mean declared probability against actual fraction of positives")
frac, mean = calibration_curve(yte, prob["random forest (100 trees)"], n_bins=5)
for f, m_ in zip(frac, mean):
    print(f"  declared {m_:.3f}  ->  actual {f:.3f}")

# Output:
# model                                Brier  log-loss     ECE     accuracy
# logistic regression                 0.1349    0.4361  0.0389       0.8187
# random forest (100 trees)           0.0967    0.3577  0.0384       0.8806
# forest + isotonic calibration       0.0944    0.3265  0.0182       0.8792
#
# reliability table of the forest (uncalibrated): mean declared probability against actual fraction of positives
#   declared 0.086  ->  actual 0.078
#   declared 0.289  ->  actual 0.192
#   declared 0.501  ->  actual 0.530
#   declared 0.716  ->  actual 0.779
#   declared 0.930  ->  actual 0.933

The forest has a better Brier than the logistic regression (0.0967 against 0.1349) and a higher accuracy (0.881 against 0.819), but the reliability table shows where its probabilities do not add up: when it states 0.289 the positive happens 19.2% of the time. The 15-bin ECE, however, is almost identical for the two models (0.0384 and 0.0389): the summary number does not say where each one goes wrong or in which direction, so you also need the table.

Isotonic calibration halves the ECE (0.0182) and improves the log-loss from 0.3577 to 0.3265, and the Brier from 0.0967 to 0.0944 (2.4% less), without worsening the accuracy appreciably (0.879 against 0.881). The improvement on the Brier is small: calibration does not change the ordering much, it corrects the scale.

More classes: micro, macro and weighted

With more than two classes F1 is computed for each class (the class against all the others) and then averaged in three ways. The macro mean is the simple mean of the classes, so every class weighs the same; the weighted one weighs each class by its number of examples; the micro one sums TP, FP and FN of all the classes and computes F1 just once.

The listing generates 3 classes with proportions 80%, 15% and 5%, trains a logistic regression and prints precision, recall and F1 per class and the three means, recomputing macro and weighted by hand.

# Three imbalanced classes (80% / 15% / 5%): micro, macro and weighted tell different stories
import numpy as np
from sklearn.datasets import make_classification
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score, f1_score, precision_recall_fscore_support
from sklearn.model_selection import train_test_split

X, y = make_classification(n_samples=20000, n_features=10, n_informative=4, n_redundant=0, n_classes=3,
                           n_clusters_per_class=1, weights=[0.8, 0.15, 0.05], class_sep=0.9, random_state=0)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.5, random_state=0, stratify=y)
p = LogisticRegression(max_iter=2000).fit(Xtr, ytr).predict(Xte)

P, R, F, S = precision_recall_fscore_support(yte, p)
print("class  precision  recall     F1  support")
for k in range(3):
    print(f"{k:5d}  {P[k]:9.3f} {R[k]:7.3f} {F[k]:6.3f} {S[k]:8d}")

print(f"\nF1 macro    (simple mean of the 3 classes): {f1_score(yte, p, average='macro'):.4f}   by hand {F.mean():.4f}")
print(f"F1 weighted (mean weighted by support):      {f1_score(yte, p, average='weighted'):.4f}   by hand {(F * S).sum() / S.sum():.4f}")
print(f"F1 micro    (counts TP, FP, FN of all classes): {f1_score(yte, p, average='micro'):.4f}")
print(f"accuracy: {accuracy_score(yte, p):.4f}  -> micro-F1 = accuracy in single-label classification: {np.isclose(f1_score(yte, p, average='micro'), accuracy_score(yte, p))}")

# Output:
# class  precision  recall     F1  support
#     0      0.976   0.997  0.986     7953
#     1      0.977   0.974  0.976     1517
#     2      0.939   0.636  0.758      530
#
# F1 macro    (simple mean of the 3 classes): 0.9067   by hand 0.9067
# F1 weighted (mean weighted by support):      0.9725   by hand 0.9725
# F1 micro    (counts TP, FP, FN of all classes): 0.9745
# accuracy: 0.9745  -> micro-F1 = accuracy in single-label classification: True

The rare class, class 2, has precision 0.939 but recall 0.636 and F1 0.758: the model misses it one time in three. The three means tell different stories: the weighted one gives 0.9725, dominated by the two large classes; the micro one 0.9745, which coincides with accuracy (it is always so when every example has a single label); the macro one 0.9067, the lowest, because the small class weighs a third.

If you care about the rare classes, the macro and the per-class recall are the metrics that see them; if you care about the fraction of correct cases, accuracy.

Regression: RMSE, MAE and R²

For a numeric value the error of each prediction is the residual $e_i = y_i - \hat{y}_i$. The MAE is the mean of $|e_i|$: in the unit of the target, every error weighs in proportion. The RMSE is $\sqrt{\text{mean}(e_i^2)}$, also in the unit of the target but with the squares under the root: large errors weigh more than small ones.

R² is $1 - \sum e_i^2 / \sum (y_i - \overline{y})^2$, that is the fraction of variance explained relative to a model that predicts the mean. The listing uses the diabetes dataset (442 patients, target: disease progression index), compares a linear regression with the model that always predicts the training mean, recomputes R² by hand and then changes a single value of the test target by adding 400.

# Regression: RMSE, MAE and R^2 on the diabetes dataset, and what a single huge error does
import numpy as np
from sklearn.datasets import load_diabetes
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
from sklearn.model_selection import train_test_split

X, y = load_diabetes(return_X_y=True)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=0)
p = LinearRegression().fit(Xtr, ytr).predict(Xte)
base = np.full_like(yte, ytr.mean())                     # model that always predicts the training mean

def row(name, y_true, y_pred):
    rmse = np.sqrt(mean_squared_error(y_true, y_pred)); mae = mean_absolute_error(y_true, y_pred)
    print(f"{name:36s} RMSE {rmse:7.2f}  MAE {mae:7.2f}  R2 {r2_score(y_true, y_pred):7.3f}")

print("standard deviation of the target on the test set: %.2f" % yte.std())
row("always predicts the mean", yte, base)
row("linear regression", yte, p)

# R^2 by hand: 1 - sum of squared residuals / sum of squares around the test mean
r2 = 1 - ((yte - p) ** 2).sum() / ((yte - yte.mean()) ** 2).sum()
print("R2 by hand: %.3f" % r2)

# A single large error: I add 400 to the target of a single test sample
y_dirty = yte.copy(); y_dirty[0] += 400
row("linear, one wrong target +400", y_dirty, p)

# Output:
# standard deviation of the target on the test set: 71.42
# always predicts the mean             RMSE   71.42  MAE   59.53  R2  -0.000
# linear regression                    RMSE   55.65  MAE   44.62  R2   0.393
# R2 by hand: 0.393
# linear, one wrong target +400        RMSE   69.20  MAE   47.63  R2   0.345

The model that predicts the mean has RMSE 71.42 (the standard deviation of the target) and practically null R²; the linear regression drops to RMSE 55.65, MAE 44.62 and R² 0.393. The hand R² coincides. The dirty sample changes the two metrics differently: the RMSE goes from 55.65 to 69.20 (+24%), the MAE from 44.62 to 47.63 (+7%).

One huge error out of 133 samples moves the RMSE by 24% and the MAE by 7%. If large errors really are more serious, the RMSE is the right metric; if they are noise in the data, the MAE is more stable. R² has no «good» value in absolute terms: 0.393 here can be excellent or poor, depending on the problem.

How far to trust a number: the bootstrap interval

Every metric computed on a test set is an estimate, and it would change if the test set were another. The bootstrap reproduces this on the computer: it resamples the test set with replacement (same size), recomputes the metric, repeats 500 times, and takes the 2.5th and 97.5th percentiles as the 95% interval.

The listing draws test sets of 200, 1,000, 5,000 and 20,000 rows from a pool, trains a logistic regression and a random forest of 60 trees, and for each test set prints the AUC of the two models with its interval and the bootstrap interval of the difference between the two, resampling the same rows for both.

# How reliable is an AUC? 95% bootstrap interval as the test set grows, and comparison between two models
import numpy as np
from sklearn.datasets import make_classification
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score
from sklearn.model_selection import train_test_split

X, y = make_classification(n_samples=30000, n_features=10, n_informative=5, flip_y=0.1, class_sep=1.0, random_state=0)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=20000, random_state=0, stratify=y)
lr = LogisticRegression(max_iter=1000).fit(Xtr, ytr)
rf = RandomForestClassifier(n_estimators=60, random_state=0, n_jobs=1).fit(Xtr, ytr)
p_lr, p_rf = lr.predict_proba(Xte)[:, 1], rf.predict_proba(Xte)[:, 1]
rng = np.random.default_rng(0)

def bootstrap(n, B=500):
    """Draws a test set of n rows and resamples it with replacement B times."""
    idx0 = rng.choice(len(yte), n, replace=False)
    a, b, t = p_lr[idx0], p_rf[idx0], yte[idx0]
    auc_lr, auc_rf = [], []
    for _ in range(B):
        i = rng.integers(0, n, n)
        if t[i].min() == t[i].max():
            continue
        auc_lr.append(roc_auc_score(t[i], a[i])); auc_rf.append(roc_auc_score(t[i], b[i]))
    auc_lr, auc_rf = np.array(auc_lr), np.array(auc_rf)
    return roc_auc_score(t, a), np.percentile(auc_lr, [2.5, 97.5]), roc_auc_score(t, b), np.percentile(auc_rf, [2.5, 97.5]), np.percentile(auc_rf - auc_lr, [2.5, 97.5])

print(f"{'n test':>7s}  {'logistic AUC [95% bootstrap]':>32s}  {'forest AUC [95% bootstrap]':>30s}  {'difference forest - logistic [95%]':>38s}")
for n in (200, 1000, 5000, 20000):
    a, ca, b, cb, cd = bootstrap(n)
    print(f"{n:7d}  {a:.3f} [{ca[0]:.3f}, {ca[1]:.3f}] (width {ca[1]-ca[0]:.3f})  {b:.3f} [{cb[0]:.3f}, {cb[1]:.3f}]  [{cd[0]:+.3f}, {cd[1]:+.3f}]")

# Output:
#  n test      logistic AUC [95% bootstrap]      forest AUC [95% bootstrap]      difference forest - logistic [95%]
#     200  0.893 [0.841, 0.942] (width 0.101)  0.922 [0.877, 0.962]  [+0.004, +0.058]
#    1000  0.933 [0.914, 0.947] (width 0.033)  0.955 [0.940, 0.967]  [+0.011, +0.034]
#    5000  0.915 [0.906, 0.923] (width 0.017)  0.938 [0.930, 0.945]  [+0.018, +0.029]
#   20000  0.915 [0.910, 0.919] (width 0.009)  0.939 [0.935, 0.942]  [+0.021, +0.027]

The interval of the logistic regression's AUC narrows from a width of 0.101 with 200 rows to 0.009 with 20,000, roughly with the square root of the number of samples (4 times the samples, about half the width). With 200 rows the intervals of the two AUCs (0.893 and 0.922) overlap, but the interval of the difference computed on the same rows, from +0.004 to +0.058, excludes zero by a small margin: looking at the two intervals separately would say «indistinguishable», the paired comparison says the forest is better, by a little.

With 5,000 rows the interval of the difference narrows to [+0.018; +0.029]. The test sets of different size are drawn separately, so the point AUCs are not monotonic (0.893, 0.933, 0.915): that is why a number alone cannot be read. A sentence like «model B is better than A» needs at least an interval on the difference, computed on the same examples.

Worked exercise

A classifier, on 1,000 examples, produces 40 true positives, 10 false positives, 20 false negatives and 930 true negatives. Compute by hand accuracy, precision, recall, specificity, F1 and F2, and the accuracy of whoever always answers «no». What does the comparison tell you?

# Exercise: from the confusion matrix to the metrics
tp, fp, fn, tn = 40, 10, 20, 930
n = tp + fp + fn + tn
print("samples:", n, " true positives:", tp + fn)
print("accuracy:", (tp + tn) / n)
P, R = tp / (tp + fp), tp / (tp + fn)
print("precision:", round(P, 4), " recall:", round(R, 4), " specificity:", round(tn / (tn + fp), 4))
print("F1:", round(2 * P * R / (P + R), 4), " F2:", round(5 * P * R / (4 * P + R), 4))
print("accuracy of whoever always answers 'no':", (tn + fp) / n)

# Output:
# samples: 1000  true positives: 60
# accuracy: 0.97
# precision: 0.8  recall: 0.6667  specificity: 0.9894
# F1: 0.7273  F2: 0.6897
# accuracy of whoever always answers 'no': 0.94

What you can do now

After this article you can do five things, and you can check them with a computation:

• read a confusion matrix and derive accuracy, precision, recall, F1 and F-beta by hand;
• explain why the AUC does not change with prevalence and the average precision does, and which one to report;
• choose a threshold with a declared cost function, on a validation set;
• measure calibration with Brier, ECE and a reliability table, and correct it if needed;
• attach a bootstrap interval to every metric you use to compare two models.

What the five items have in common is that no metric is worth anything alone: each answers a question and the choice of the question comes before the model. The next step is to take the model out of the notebook and keep measuring it after release: deployment in production is the last article of the series. If you want to revisit the model we evaluated most here, logistic regression is article 3; for tree models, trees and forests. The same metrics serve elsewhere in the series: SVMs need calibrated probabilities, clustering is evaluated with ARI and silhouette, PCA is judged by the reconstruction error and neural networks by the learning curve.

The same series also includes: 01 - Introduction to Machine Learning: Taxonomy and Workflow.