Scikit-Learn Model Evaluation

A model is only as good as the way it is measured. Every example here ran with scikit-learn 1.9.1, mostly on the built-in breast cancer dataset (569 tumours, 212 malignant, 357 benign).

Cross-Validation

cross_val_score trains and scores the model once per fold and returns every score.

from sklearn.datasets import load_breast_cancer
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import StratifiedKFold, cross_val_score

X, y = load_breast_cancer(return_X_y=True)
model = RandomForestClassifier(n_estimators=100, random_state=42)

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
scores = cross_val_score(model, X, y, cv=cv)
print("accuracy per fold:", scores.round(4))
print(f"mean {scores.mean():.4f}  std {scores.std():.4f}")

$ python cv.py
accuracy per fold: [0.9649 0.9386 0.9561 0.9474 0.9735]
mean 0.9561  std 0.0123

The five folds range from 0.939 to 0.974, a spread of more than three points from shuffling alone. A single train/test split would have reported one of those numbers as “the” accuracy. The standard deviation tells you how much weight to put on small differences between models.

cross_validate scores several metrics in the same pass and also records timing.

from sklearn.datasets import load_breast_cancer
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import StratifiedKFold, cross_validate

X, y = load_breast_cancer(return_X_y=True)
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
res = cross_validate(RandomForestClassifier(n_estimators=100, random_state=42), X, y,
                     cv=cv, scoring=["accuracy", "f1", "roc_auc"])
for metric in ["accuracy", "f1", "roc_auc"]:
    vals = res[f"test_{metric}"]
    print(f"{metric:9s} {vals.mean():.4f} +/- {vals.std():.4f}")
print("keys:", sorted(res))

$ python cv_multi.py
accuracy  0.9561 +/- 0.0123
f1        0.9651 +/- 0.0097
roc_auc   0.9885 +/- 0.0083
keys: ['fit_time', 'score_time', 'test_accuracy', 'test_f1', 'test_roc_auc']

The accuracy row matches cv.py exactly because the folds and seed are the same. f1 and roc_auc come from the same five fits, so trying another metric costs nothing extra. fit_time and score_time are included too, which is useful when comparing slow models.

Precision and Recall, and Which Class They Measure

In this dataset, label 0 is malignant and label 1 is benign. scikit-learn’s binary metrics report label 1 unless told otherwise.

from sklearn.datasets import load_breast_cancer
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.metrics import accuracy_score, precision_score, recall_score
from sklearn.model_selection import train_test_split

X, y = load_breast_cancer(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42, stratify=y)
pred = GradientBoostingClassifier(random_state=42).fit(X_train, y_train).predict(X_test)

print(f"accuracy                        {accuracy_score(y_test, pred):.4f}")
print(f"precision (default: benign)     {precision_score(y_test, pred):.4f}")
print(f"recall    (default: benign)     {recall_score(y_test, pred):.4f}")
print(f"precision (pos_label=0: malignant) {precision_score(y_test, pred, pos_label=0):.4f}")
print(f"recall    (pos_label=0: malignant) {recall_score(y_test, pred, pos_label=0):.4f}")

$ python metrics.py
accuracy                        0.9561
precision (default: benign)     0.9467
recall    (default: benign)     0.9861
precision (pos_label=0: malignant) 0.9744
recall    (pos_label=0: malignant) 0.9048

The default metrics describe benign tumours, the class nobody is worried about missing. Asked about malignant ones, recall drops to 0.905: about 1 malignant tumour in 10 was called benign, a figure that recall_score’s 0.986 default hid. Always check which label is the positive class before reading precision or recall.

The Confusion Matrix

Every binary metric is built from four counts, which confusion_matrix returns directly.

from sklearn.datasets import load_breast_cancer
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.metrics import confusion_matrix
from sklearn.model_selection import train_test_split

X, y = load_breast_cancer(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42, stratify=y)
pred = GradientBoostingClassifier(random_state=42).fit(X_train, y_train).predict(X_test)

cm = confusion_matrix(y_test, pred)
print(cm)
tn, fp, fn, tp = cm.ravel()
print(f"tn={tn} fp={fp} fn={fn} tp={tp}   (positive = label 1 = benign)")

$ python confusion.py
[[38  4]
 [ 1 71]]
tn=38 fp=4 fn=1 tp=71   (positive = label 1 = benign)

Rows are the true labels and columns the predictions, in label order (0 = malignant, 1 = benign). The 4 in the top-right are malignant tumours predicted benign, the costly mistake here. Because benign is label 1, scikit-learn names them “false positives”, a reminder that “positive” just means label 1.

ROC Curves and Choosing an Operating Point

roc_curve lists, for every threshold, the true-positive and false-positive rates.

import numpy as np
from sklearn.datasets import load_breast_cancer
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score, roc_curve
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

X, y = load_breast_cancer(return_X_y=True)
y = 1 - y                                   # make malignant the positive class
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=0, stratify=y)
model = make_pipeline(StandardScaler(), LogisticRegression(max_iter=1000))
prob = model.fit(X_train, y_train).predict_proba(X_test)[:, 1]

fpr, tpr, thr = roc_curve(y_test, prob)
print(f"ROC AUC {roc_auc_score(y_test, prob):.4f},  {len(thr)} thresholds")
i = np.argmax(tpr >= 0.99)                 # first threshold catching 99% of malignant
print(f"to catch 99% of malignant: threshold {thr[i]:.3f}, false-positive rate {fpr[i]:.3f}")

$ python roc.py
ROC AUC 0.9917,  14 thresholds
to catch 99% of malignant: threshold 0.032, false-positive rate 0.224

Flipping the labels makes malignant the positive class, so the rates now describe cancers. Catching 99% of them requires dropping the threshold to about 0.03, which also flags about 22% of benign tumours. ROC AUC summarises all thresholds in one number, but a deployment has to pick one point on the curve, and this is how to read that trade-off.

Regression Metrics

Different regression metrics weigh errors differently, so report more than one.

from sklearn.datasets import fetch_california_housing
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.metrics import (mean_absolute_error, mean_absolute_percentage_error,
                             r2_score, root_mean_squared_error)
from sklearn.model_selection import train_test_split

X, y = fetch_california_housing(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
pred = HistGradientBoostingRegressor(random_state=0).fit(X_train, y_train).predict(X_test)

print(f"MAE   {mean_absolute_error(y_test, pred):.4f}")
print(f"RMSE  {root_mean_squared_error(y_test, pred):.4f}")
print(f"MAPE  {mean_absolute_percentage_error(y_test, pred):.2%}")
print(f"R2    {r2_score(y_test, pred):.4f}")
cheap = y_test < 0.6
print(f"MAPE on the cheapest districts only: {mean_absolute_percentage_error(y_test[cheap], pred[cheap]):.2%}")

$ python regression_metrics.py
MAE   0.3095
RMSE  0.4621
MAPE  17.96%
R2    0.8370
MAPE on the cheapest districts only: 62.97%

MAE is the average miss in the target’s units ($31,000 here), and RMSE is larger because it squares errors and so weighs big misses more. MAPE says 18% overall, but on the cheapest districts it jumps to 63%: dividing by a small true value inflates percentage errors. Prefer MAE or RMSE unless every target is comfortably away from zero.

Imbalanced Classes

With 1% positives, the default settings catch few of them. Class weights and a lower threshold both shift the balance toward recall.

from sklearn.datasets import make_classification
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import precision_score, recall_score
from sklearn.model_selection import train_test_split

X, y = make_classification(n_samples=10_000, weights=[0.99, 0.01], random_state=42)
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, stratify=y, random_state=42)
print("positives in test set:", y_test.sum())

for cw in [None, "balanced"]:
    prob = RandomForestClassifier(class_weight=cw, random_state=42).fit(
        X_train, y_train).predict_proba(X_test)[:, 1]
    for t in [0.5, 0.2]:
        pred = prob >= t
        print(f"class_weight={str(cw):8s} threshold {t}: recall {recall_score(y_test, pred):.2f}"
              f"  precision {precision_score(y_test, pred, zero_division=0):.2f}")

$ python imbalanced.py
positives in test set: 29
class_weight=None     threshold 0.5: recall 0.14  precision 1.00
class_weight=None     threshold 0.2: recall 0.31  precision 0.53
class_weight=balanced threshold 0.5: recall 0.14  precision 0.80
class_weight=balanced threshold 0.2: recall 0.38  precision 0.22

With only 29 positives, the default forest at 0.5 catches 4 of them (recall 0.14), though everything it flags is real. Lowering the threshold to 0.2 roughly doubles recall at the cost of precision. class_weight="balanced" barely moves recall here (0.14, unchanged): a random forest still grows each tree to pure leaves. The threshold is the stronger lever here, and it can be chosen after training.

Calibration

A model can rank examples well while its probability values are off. GaussianNB tends to push probabilities toward 0 and 1 more than the data justifies.

from sklearn.calibration import CalibratedClassifierCV
from sklearn.datasets import make_classification
from sklearn.metrics import brier_score_loss
from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB

X, y = make_classification(n_samples=3000, n_features=20, n_informative=15,
                           n_redundant=5, random_state=42)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

uncal_prob = GaussianNB().fit(X_train, y_train).predict_proba(X_test)[:, 1]
cal = CalibratedClassifierCV(GaussianNB(), method="isotonic", cv=5).fit(X_train, y_train)
cal_prob = cal.predict_proba(X_test)[:, 1]

print(f"Uncalibrated Brier score: {brier_score_loss(y_test, uncal_prob):.4f}")
print(f"Calibrated Brier score:   {brier_score_loss(y_test, cal_prob):.4f}")

$ python calibration_brier.py
Uncalibrated Brier score: 0.1456
Calibrated Brier score:   0.1410

The Brier score is the mean squared error between the predicted probability and the actual 0 or 1, so lower is better. Calibration lowers it from 0.1456 to 0.1410. CalibratedClassifierCV fits GaussianNB on cross-validation folds and learns the correction on held-out predictions, so the calibration step never sees a prediction on the base model’s own training rows.

The Brier score doesn’t show where the miscalibration is. A reliability table bins the predictions and compares each bin’s average prediction with its actual positive rate.

from sklearn.calibration import CalibratedClassifierCV, calibration_curve
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB

X, y = make_classification(n_samples=3000, n_features=20, n_informative=15,
                           n_redundant=5, random_state=42)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

models = {"Uncalibrated": GaussianNB(),
          "Calibrated": CalibratedClassifierCV(GaussianNB(), method="isotonic", cv=5)}
for label, model in models.items():
    prob = model.fit(X_train, y_train).predict_proba(X_test)[:, 1]
    frac_pos, mean_pred = calibration_curve(y_test, prob, n_bins=5, strategy="quantile")
    print(f"{label}:")
    for pred, actual in zip(mean_pred, frac_pos):
        print(f"  predicted {pred:.2f}  ->  actually positive {actual:.2f}")

$ python calibration_curve.py
Uncalibrated:
  predicted 0.01  ->  actually positive 0.11
  predicted 0.10  ->  actually positive 0.23
  predicted 0.41  ->  actually positive 0.45
  predicted 0.80  ->  actually positive 0.77
  predicted 0.98  ->  actually positive 0.99
Calibrated:
  predicted 0.10  ->  actually positive 0.12
  predicted 0.24  ->  actually positive 0.23
  predicted 0.39  ->  actually positive 0.45
  predicted 0.75  ->  actually positive 0.78
  predicted 0.98  ->  actually positive 1.00

The uncalibrated model’s lowest bin is the clearest tell: it predicts 0.01 on average, but 11% of those examples are actually positive, so it is overconfident at the low end. After calibration the predicted and actual columns track each other closely in every bin, which is the same improvement the Brier score summarised as one number.

Practice

from sklearn.calibration import CalibratedClassifierCV
from sklearn.datasets import make_classification
from sklearn.metrics import brier_score_loss
from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB

X, y = make_classification(n_samples=3000, n_features=20, n_informative=15,
                           n_redundant=5, random_state=42)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)
for method in ["sigmoid", "isotonic"]:
    cal = CalibratedClassifierCV(GaussianNB(), method=method, cv=5).fit(X_train, y_train)
    print(f"{method:8s} Brier {brier_score_loss(y_test, cal.predict_proba(X_test)[:, 1]):.4f}")

$ python p1_sigmoid.py
sigmoid  Brier 0.1424
isotonic Brier 0.1410

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

X, y = make_classification(n_samples=3000, n_features=20, n_informative=15,
                           n_redundant=5, random_state=42)
_, _, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)
base = np.full(len(y_test), y_train.mean())
print(f"base rate {y_train.mean():.3f}  Brier {brier_score_loss(y_test, base):.4f}")

$ python p2_baseline.py
base rate 0.498  Brier 0.2500

from sklearn.datasets import make_classification
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import average_precision_score, roc_auc_score
from sklearn.model_selection import train_test_split

X, y = make_classification(n_samples=10_000, weights=[0.99, 0.01], random_state=42)
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, stratify=y, random_state=42)
prob = RandomForestClassifier(random_state=42).fit(X_train, y_train).predict_proba(X_test)[:, 1]
print(f"ROC AUC           {roc_auc_score(y_test, prob):.4f}")
print(f"average precision {average_precision_score(y_test, prob):.4f}")
print(f"positive rate     {y_test.mean():.4f}")

$ python p3_ap_vs_auc.py
ROC AUC           0.8287
average precision 0.3465
positive rate     0.0145