6.1 Classification

A classifier gives each case a score, and a threshold turns the score into a decision. Where to put the threshold depends on which mistake is worse.

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.axes import Axes
from matplotlib.colors import ListedColormap
from scipy import stats
from sklearn import datasets
from sklearn.base import BaseEstimator
from sklearn.calibration import calibration_curve
from sklearn.discriminant_analysis import (
    LinearDiscriminantAnalysis,
    QuadraticDiscriminantAnalysis,
)
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import confusion_matrix, precision_recall_curve, roc_auc_score, roc_curve
from sklearn.model_selection import cross_val_predict, train_test_split
from sklearn.naive_bayes import GaussianNB
from sklearn.neighbors import KNeighborsClassifier
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

from statly import style

style.use()
rng = np.random.default_rng(1)

cancer = datasets.load_breast_cancer(
    as_frame=True
)  # 569 breast tumour samples, 30 measurements
X = cancer.data
y = (cancer.target == 0).astype(int).to_numpy()  # 1 = malignant, 0 = benign
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, stratify=y, random_state=1
)
CLASS_COLORS = [style.GREEN, style.CORAL]  # benign, malignant

6.1.1 Drawing a boundary

Doctors take a sample from a breast lump and measure the cell nuclei under a microscope: their size, shape and texture, 30 numbers in all. For 569 samples we also know the answer: benign or malignant. Can the measurements tell the two apart?

A classifier learns from labelled cases and draws a boundary between the classes. Here are three, using just two of the measurements, so we can see the boundary. Each is trained on 70% of the samples; the dots are the other 30%, which the classifiers have not seen.

FEATURES = ["worst radius", "worst concave points"]
gx, gy = np.meshgrid(np.linspace(7, 37, 300), np.linspace(0, 0.3, 300))
grid = np.c_[gx.ravel(), gy.ravel()]
CLASSIFIERS_2D = [
    ("Logistic regression", LogisticRegression()),
    ("1 nearest neighbour", KNeighborsClassifier(1)),
    ("25 nearest neighbours", KNeighborsClassifier(25)),
]


def plot_boundary(ax: Axes, name: str, classifier: BaseEstimator) -> None:
    """The zones a classifier assigns to each class, with the test cases on top."""
    model = make_pipeline(StandardScaler(), classifier).fit(
        X_train[FEATURES].to_numpy(), y_train
    )
    zones = model.predict(grid).reshape(gx.shape)
    ax.contourf(
        gx, gy, zones, levels=[-0.5, 0.5, 1.5], cmap=ListedColormap(CLASS_COLORS), alpha=0.15
    )
    test = X_test[FEATURES].to_numpy()
    ax.scatter(test[:, 0], test[:, 1], c=np.array(CLASS_COLORS)[y_test], s=10, alpha=0.8)
    errors = (model.predict(test) != y_test).mean()
    style.note(ax, f"wrong on new cases: {errors:.0%}", x=0.03, color=style.INK)
    ax.set(title=name, xlabel="worst radius", ylabel="worst concave points")


fig, axes = style.subplots(3, height=3.1, width=10.5)
for ax, (name, classifier) in zip(axes, CLASSIFIERS_2D):
    plot_boundary(ax, name, classifier)
fig.tight_layout()
plt.show()

Figure 6.1.1: Benign (green) and malignant (red) tumours in the test set, by two measurements. The shading shows what each classifier would say for a new tumour. Logistic regression draws a straight boundary; one nearest neighbour copies the closest training case and follows every quirk; 25 neighbours vote and give a smoother boundary.

The trade-off of Chapter 1.3 is back. One neighbour is the most flexible classifier possible and draws islands around single training cases, which do not help on new ones. Twenty-five neighbours and the logistic line are smoother and make fewer errors on new tumours.

6.1.2 The confusion matrix

With all 30 measurements, logistic regression does even better. “Errors” is too blunt a summary, though, because the two kinds of error are not equally bad: calling a malignant tumour benign delays treatment, calling a benign one malignant leads to more tests. The confusion matrix counts all four outcomes on the test set:

full_model = make_pipeline(StandardScaler(), LogisticRegression(max_iter=1000)).fit(
    X_train, y_train
)
scores = full_model.predict_proba(X_test)[:, 1]  # predicted chance of malignant
matrix = confusion_matrix(y_test, scores >= 0.5)
NAMES = ["benign", "malignant"]


def plot_matrix(ax: Axes) -> None:
    """The four counts as coloured squares, correct ones in blue."""
    for truth in range(2):
        for call in range(2):
            color = style.BLUE if truth == call else style.CORAL
            ax.add_patch(
                plt.Rectangle(
                    (call, 1 - truth),
                    1,
                    1,
                    color=color,
                    alpha=0.15 + 0.6 * matrix[truth, call] / matrix.max(),
                )
            )
            ax.text(
                call + 0.5,
                1.5 - truth,
                str(matrix[truth, call]),
                ha="center",
                va="center",
                fontsize=18,
            )
    ax.set_xticks([0.5, 1.5], [f"called {n}" for n in NAMES])
    ax.set_yticks([1.5, 0.5], [f"truly {n}" for n in NAMES])
    ax.set(xlim=(0, 2), ylim=(0, 2))
    ax.grid(False)
    ax.tick_params(length=0)


fig, ax = plt.subplots(figsize=(4.4, 3))
plot_matrix(ax)
plt.show()

Figure 6.1.2: Logistic regression with all 30 measurements, on the 171 test tumours it has not seen. The diagonal holds the correct calls.

From the four counts come the numbers everyone uses. Call the malignant class “positive” and write TP for true positives (malignant, called malignant), FN for false negatives (malignant, called benign), TN and FP likewise:

Measure	Formula	Here
accuracy	(TP + TN) / all	164 of 171 = 96%
sensitivity, recall, true positive rate	TP / (TP + FN)	58 of 64 = 91%
specificity, true negative rate	TN / (TN + FP)	106 of 107 = 99%
precision, positive predictive value	TP / (TP + FP)	58 of 59 = 98%
false positive rate	FP / (TN + FP) = 1 − specificity	1 of 107 = 1%
F1 score	2 · precision · recall / (precision + recall)	0.94

The F1 score is the harmonic mean of precision and recall: high only if both are.

Which one matters depends on the use. A screening test needs high sensitivity: missing a cancer is worse than an extra biopsy. A test that decides on surgery needs high precision.

6.1.3 Moving the threshold

Logistic regression does not really say “malignant”. It gives a probability, and we call everything above 50% malignant. Nothing makes 50% special. A lower threshold catches more cancers and alarms more healthy people; a higher one the reverse. The ROC curve shows every threshold at once: the share of sick cases caught against the share of healthy cases flagged.

A useless test, a coin toss, follows the diagonal. A perfect test goes straight up to the top left corner. The area under the curve (AUC) sums up a test across all thresholds: 0.5 for a coin, 1 for perfect.

Building the ROC curve and its area, step by step

With real data the ROC curve is a staircase. Sort the cases by their score, from highest to lowest, and lower the threshold past them one at a time. Each sick case passed moves the curve up by 1/(number of sick cases); each healthy case passed moves it right by 1/(number of healthy cases). Twenty test tumours show it:

one_feature = make_pipeline(StandardScaler(), LogisticRegression()).fit(
    X_train[["mean texture"]], y_train
)
pick = np.r_[np.where(y_test == 1)[0][:10], np.where(y_test == 0)[0][:10]]
few_scores = one_feature.predict_proba(X_test[["mean texture"]].iloc[pick])[:, 1]
few_labels = y_test[pick]
order = np.argsort(-few_scores)


def plot_sorted(ax: Axes) -> None:
    """The scores from high to low, coloured by the true class."""
    for rank, i in enumerate(order):
        ax.plot(few_scores[i], -rank, "o", ms=8, color=CLASS_COLORS[few_labels[i]])
    ax.set(
        title="Sorted scores",
        xlabel="predicted chance of malignant",
        yticks=[],
        xlim=(-0.05, 1.05),
    )


def plot_staircase(ax: Axes) -> None:
    """The ROC staircase traced through the sorted cases."""
    x, y = [0.0], [0.0]
    for i in order:
        if few_labels[i] == 1:
            x.append(x[-1])
            y.append(y[-1] + 1 / 10)
        else:
            x.append(x[-1] + 1 / 10)
            y.append(y[-1])
    ax.fill_between(x, y, step=None, color=style.BLUE, alpha=0.15, lw=0)
    ax.plot(x, y, color=style.BLUE, lw=2, marker="o", ms=3)
    ax.plot([0, 1], [0, 1], color=style.INK, ls=":", lw=0.8)
    style.note(
        ax, f"AUC = {roc_auc_score(few_labels, few_scores):.2f}", x=0.45, y=0.2, color=style.INK
    )
    ax.set(
        title="The staircase",
        xlabel="benign flagged",
        ylabel="malignant caught",
        aspect="equal",
        xlim=(0, 1),
        ylim=(0, 1.02),
    )


def plot_pairs(ax: Axes) -> None:
    """Every malignant-benign pair: blue if ranked correctly."""
    sick = np.sort(few_scores[few_labels == 1])[::-1]
    healthy = np.sort(few_scores[few_labels == 0])[::-1]
    correct = sick[:, None] > healthy[None, :]
    ax.imshow(
        correct, cmap=ListedColormap([style.CORAL, style.BLUE]), alpha=0.6, vmin=0, vmax=1
    )
    style.note(ax, f"correct pairs: {correct.mean():.0%}", x=0.03, y=-0.08, color=style.INK)
    ax.set(
        title="All 100 pairs",
        xlabel="benign, highest score first",
        ylabel="malignant, highest first",
        xticks=[],
        yticks=[],
    )
    ax.grid(False)


fig, axes = style.subplots(3, height=3.3, width=10.5)
plot_sorted(axes[0])
plot_staircase(axes[1])
plot_pairs(axes[2])
fig.tight_layout()
plt.show()

Figure 6.1.3: Twenty test tumours, scored by logistic regression on a single, weaker measurement (mean texture), so the classes overlap. Left: the scores, from highest (top) to lowest; red are malignant, green benign. Middle: lowering the threshold past each tumour in turn traces the ROC staircase: up for a malignant tumour, right for a benign one; the shaded area is the AUC. Right: every pair of one malignant and one benign tumour; a cell is blue if the malignant one scored higher. The share of blue cells is exactly the AUC.

Each column of the staircase area belongs to one benign tumour, and its height is the share of malignant tumours scored above it. So the area is the share of malignant-benign pairs in the right order:

$$\text{AUC} = P(\text{score of a random positive} > \text{score of a random negative})$$

the probability of superiority of Chapter 3.4. An AUC of 0.9 means that nine times out of ten the model ranks a random sick case above a random healthy one. It says nothing about any particular threshold, nor about whether the probabilities themselves are right (calibration, Section 6.1.6).

6.1.4 Comparing classifiers

ROC curves make classifiers easy to compare, independent of any threshold:

CLASSIFIERS = [
    ("logistic regression", LogisticRegression(max_iter=1000), style.CORAL),
    ("linear discriminant", LinearDiscriminantAnalysis(), style.BLUE),
    ("15 nearest neighbours", KNeighborsClassifier(15), style.GREEN),
    ("naive Bayes", GaussianNB(), style.GOLD),
]


def plot_roc(
    ax: Axes,
    model: BaseEstimator,
    features: list[str] | None,
    label: str,
    color: str,
    ls: str = "-",
) -> None:
    """Fit on the training set; draw the ROC curve on the test set, with the AUC in the label."""
    columns = features or list(X.columns)
    fitted = make_pipeline(StandardScaler(), model).fit(X_train[columns], y_train)
    probability = fitted.predict_proba(X_test[columns])[:, 1]
    fpr, tpr, _ = roc_curve(y_test, probability)
    ax.plot(
        fpr,
        tpr,
        color=color,
        ls=ls,
        lw=1.8,
        label=f"{label}: {roc_auc_score(y_test, probability):.3f}",
    )


fig, ax = plt.subplots(figsize=(6.2, 4.4))
for name, model, color in CLASSIFIERS:
    plot_roc(ax, model, None, name, color)
plot_roc(ax, LogisticRegression(), FEATURES, "logistic, 2 measurements", style.GREY, "--")
ax.plot([0, 1], [0, 1], color=style.INK, lw=0.8, ls=":")
ax.legend(title="area under the curve", fontsize=8, title_fontsize=8, loc="lower right")
ax.set(
    xlabel="benign tumours flagged (false positive rate)",
    ylabel="malignant tumours caught",
    aspect="equal",
)
plt.show()

Figure 6.1.4: ROC curves on the test tumours for four classifiers using all 30 measurements, and for logistic regression with only two (dashed). All are excellent: the curves almost reach the top left corner, and two well-chosen measurements do nearly as well as thirty.

The methods differ in how they draw the boundary:

• Logistic regression models the log-odds as a straight line (Chapter 4.4).
• Linear discriminant analysis assumes each class is a normal cloud with the same shape and finds the line between the clouds. It often gives nearly the same boundary as logistic regression.
• K nearest neighbours looks up the k most similar training cases and lets them vote. It needs no shape at all, but it needs many cases and suffers when there are many measurements.
• Naive Bayes pretends the measurements are independent within each class. That is almost never true, yet it often works surprisingly well, especially with many features.

On this test set of 171 tumours, the differences between the first three are a handful of cases, well within chance. The next chapter shows how to compare methods more reliably than with a single split.

6.1.5 How the generative classifiers work

Logistic regression models the probability of each class given x directly. Linear and quadratic discriminant analysis and naive Bayes go the other way round. They model what x looks like within each class, a density $f_k(x)$ for class k, and turn it around with Bayes’ rule (Chapter 2.4), using the share of each class, $\pi_k$ (the prior):

$$P(\text{class } k \mid x) = \frac{\pi_k\, f_k(x)}{\sum_l \pi_l\, f_l(x)}$$

• LDA assumes each class is a normal cloud, all with the same covariance matrix. The boundary between two classes is then a straight line (a flat plane in more dimensions).
• QDA lets each class have its own covariance matrix; the boundaries become curves.
• Naive Bayes assumes the measurements are independent within each class, so each $f_k$ is a product of one-dimensional densities: crude, but cheap and robust with many measurements.

radius = X_train["worst radius"].to_numpy()
grid_1d = np.linspace(radius.min() - 2, radius.max() + 2, 400)
priors = np.bincount(y_train) / len(y_train)


def plot_bayes_1d(ax: Axes) -> None:
    """Class densities times priors, and the posterior probability they imply."""
    weighted = []
    pooled_sd = np.sqrt(
        sum(((radius[y_train == k] - radius[y_train == k].mean()) ** 2).sum() for k in (0, 1))
        / (len(radius) - 2)
    )
    for k, (color, name) in enumerate(zip(CLASS_COLORS, ("benign", "malignant"))):
        values = radius[y_train == k]
        density = priors[k] * stats.norm.pdf(
            grid_1d, values.mean(), pooled_sd
        )  # LDA: one shared spread
        weighted.append(density)
        ax.fill_between(grid_1d, density, color=color, alpha=0.3, lw=0, label=name)
    ax.set(
        title="Bayes' rule in one dimension",
        xlabel="worst radius",
        ylabel="prior × density",
        yticks=[],
    )
    ax.legend(fontsize=8, loc="upper right")
    twin = ax.twinx()
    twin.plot(grid_1d, weighted[1] / (weighted[0] + weighted[1]), color=style.INK, ls="--")
    twin.yaxis.set_major_formatter(style.PERCENT)
    twin.set_ylim(0, 1.05)


def plot_discriminant(ax: Axes, model: BaseEstimator, title: str) -> None:
    """Class ellipses and decision zones for a discriminant model."""
    fitted = make_pipeline(StandardScaler(), model).fit(X_train[FEATURES].to_numpy(), y_train)
    zones = fitted.predict(grid).reshape(gx.shape)
    ax.contourf(
        gx, gy, zones, levels=[-0.5, 0.5, 1.5], cmap=ListedColormap(CLASS_COLORS), alpha=0.15
    )
    test = X_test[FEATURES].to_numpy()
    ax.scatter(test[:, 0], test[:, 1], c=np.array(CLASS_COLORS)[y_test], s=8, alpha=0.7)
    ax.set(title=title, xlabel="worst radius", ylabel="worst concave points")


fig, axes = style.subplots(3, height=3.2, width=10.5)
plot_bayes_1d(axes[0])
plot_discriminant(axes[1], LinearDiscriminantAnalysis(), "LDA: a straight boundary")
plot_discriminant(axes[2], QuadraticDiscriminantAnalysis(), "QDA: a curved boundary")
fig.tight_layout()
plt.show()

Figure 6.1.5: Left: one measurement (worst radius), the normal density of each class scaled by its share of the training tumours, and the resulting probability of malignant (dashed, right axis): Bayes’ rule turns two bells into an S-curve. Middle and right: LDA and QDA on two measurements; LDA’s equal-shape ellipses give a straight boundary, QDA’s own-shape ellipses a curved one.

When the classes really are normal clouds, LDA and QDA are the best possible classifiers, and they need less data than logistic regression. When they are far from normal, logistic regression is safer. QDA needs a covariance matrix per class, many parameters, so it pays off only with enough data or truly different shapes.

6.1.6 When one class is rare, and calibration

If 1 in 100 people has a disease, a classifier that always says “healthy” is 99% accurate and completely useless. With rare classes, accuracy hides everything. The ROC curve does not: it does not depend on how common the classes are. But that is also its blind spot. The precision, what share of the alarms are real, falls with the base rate (Chapter 2.4), and the ROC curve never shows it. The precision-recall curve does:

SEPARATION, CASES = 2.0, 200_000


def simulated_scores(prevalence: float) -> tuple[np.ndarray, np.ndarray]:
    """Scores for healthy (normal around 0) and sick (normal around SEPARATION) cases."""
    sick = rng.uniform(size=CASES) < prevalence
    return rng.normal(np.where(sick, SEPARATION, 0), 1), sick.astype(int)


fig, axes = style.subplots(2, height=3.4, width=9.5)
for prevalence, color in ((0.5, style.BLUE), (0.1, style.GOLD), (0.02, style.CORAL)):
    scores, labels = simulated_scores(prevalence)
    fpr, tpr, _ = roc_curve(labels, scores)
    axes[0].plot(fpr, tpr, color=color, lw=1.8, label=f"{prevalence:.0%} sick")
    precision, recall, _ = precision_recall_curve(labels, scores)
    axes[1].plot(recall, precision, color=color, lw=1.8, label=f"{prevalence:.0%} sick")
axes[0].set(
    title="ROC", xlabel="false positive rate", ylabel="sensitivity (recall)", aspect="equal"
)
axes[1].set(
    title="Precision-recall", xlabel="recall (sensitivity)", ylabel="precision", ylim=(0, 1.02)
)
for ax in axes:
    ax.legend(fontsize=8)
fig.tight_layout()
plt.show()

Figure 6.1.6: The same simulated test (scores of sick cases two standard deviations higher than those of healthy ones) applied where the disease is common (50%), uncommon (10%) and rare (2%). Left: the ROC curves are identical. Right: the precision-recall curves are not; when the disease is rare, most alarms are false at every threshold.

With rare classes, look at precision and recall, think about the costs of each error, and set the threshold accordingly: often far below 50% when missing a case is expensive.

A last property of a classifier’s scores: are they honest probabilities? A model is calibrated if, among cases it gives 30%, about 30% are positive. Logistic regression is usually close; naive Bayes and many machine learning methods are not, and their scores need recalibrating before they can be read as probabilities:

def plot_calibration(ax: Axes) -> None:
    """Observed share against predicted probability, per bin, for two classifiers."""
    ax.plot([0, 1], [0, 1], color=style.INK, ls="--", lw=1)
    for name, model, color in (
        ("logistic regression", LogisticRegression(max_iter=1000), style.CORAL),
        ("naive Bayes", GaussianNB(), style.GOLD),
    ):
        scores = cross_val_predict(
            make_pipeline(StandardScaler(), model), X, y, cv=10, method="predict_proba"
        )[:, 1]
        observed, predicted = calibration_curve(y, scores, n_bins=8, strategy="quantile")
        ax.plot(predicted, observed, "o-", color=color, label=name)
    ax.xaxis.set_major_formatter(style.PERCENT)
    ax.yaxis.set_major_formatter(style.PERCENT)
    ax.legend(fontsize=8)
    ax.set(
        xlabel="predicted chance of malignant", ylabel="share truly malignant", aspect="equal"
    )


fig, ax = plt.subplots(figsize=(4.6, 4.2))
plot_calibration(ax)
plt.show()

Figure 6.1.7: Calibration for all 569 tumours, using cross-validated predictions (each tumour scored by a model that did not see it, Chapter 6.3) pooled into eight equal-sized bins: the share that was truly malignant against the average predicted probability. Logistic regression (red) lies close to the diagonal; naive Bayes (gold) is overconfident, pushing its probabilities towards 0 and 1.

6.1.7 More than two classes

Everything carries over to more classes: logistic regression becomes multinomial (Chapter 4.4), LDA and QDA work as before, and the confusion matrix grows to one row and one column per class. Its off-diagonal cells show which classes get mixed up:

digits = datasets.load_digits()
d_train, d_test, t_train, t_test = train_test_split(
    digits.data, digits.target, test_size=0.3, random_state=1
)
digit_model = make_pipeline(StandardScaler(), LogisticRegression(max_iter=3000)).fit(
    d_train, t_train
)
digit_matrix = confusion_matrix(t_test, digit_model.predict(d_test))


def plot_digit_matrix(ax: Axes) -> None:
    """The 10 × 10 confusion matrix with counts."""
    off = digit_matrix.copy().astype(float)
    np.fill_diagonal(off, np.nan)
    ax.imshow(np.where(np.eye(10, dtype=bool), 0.15, 0), cmap="Blues", vmin=0, vmax=1)
    for i in range(10):
        for j in range(10):
            count = digit_matrix[i, j]
            if count:
                ax.text(
                    j,
                    i,
                    str(count),
                    ha="center",
                    va="center",
                    fontsize=8.5,
                    color=style.CORAL if i != j else style.INK,
                    fontweight="bold" if i != j else "normal",
                )
    ax.set_xticks(range(10))
    ax.set_yticks(range(10))
    ax.set(xlabel="predicted digit", ylabel="true digit")
    ax.grid(False)
    accuracy = np.trace(digit_matrix) / digit_matrix.sum()
    ax.set_title(f"Accuracy {accuracy:.1%}")


fig, ax = plt.subplots(figsize=(4.8, 4.6))
plot_digit_matrix(ax)
plt.show()

Figure 6.1.8: Confusion matrix for handwritten digits (8 × 8 pixels), classified by logistic regression on 30% held-out images. Most digits land on the diagonal; the few errors (red) are scattered, such as two 3s taken for 2s.

6.1.8 In one picture

• A classifier draws a boundary; more flexible boundaries can overfit, as in regression.
• Judge it on new cases, with the confusion matrix: sensitivity TP/(TP + FN), specificity TN/(TN + FP), precision TP/(TP + FP), F1.
• The threshold decides the balance between the two errors. The ROC curve shows all thresholds; the AUC is the chance a random positive outscores a random negative.
• With rare classes, use precision-recall; check that probabilities are calibrated.
• Logistic regression models P(class | x) directly; LDA, QDA and naive Bayes model each class’s shape and use Bayes’ rule; nearest neighbours votes.
• With rare classes, accuracy misleads.

Next: how to split the data so that a model is judged honestly.