1.3 Models

A model is our best guess at the systematic part of the data. The art is to be flexible enough to catch the pattern, and rigid enough to ignore the noise.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
from matplotlib.axes import Axes
from sklearn import datasets
from sklearn.cluster import KMeans
from sklearn.decomposition import PCA
from sklearn.linear_model import LogisticRegression
from sklearn.neighbors import KNeighborsClassifier, KNeighborsRegressor
from sklearn.preprocessing import StandardScaler

from statly import style

style.use()
rng = np.random.default_rng(1)
SPECIES = ["setosa", "versicolor", "virginica"]
SPECIES_COLORS = [style.BLUE, style.GREEN, style.CORAL]

engel = (
    sm.datasets.engel.load_pandas().data
)  # income and food spending, Belgian households (1857)
iris = datasets.load_iris(as_frame=True).frame  # flowers of three iris species
co2 = sm.datasets.co2.load_pandas().data.interpolate()  # weekly CO2 at Mauna Loa, 1958 to 2001
diabetes = datasets.load_diabetes(
    as_frame=True, scaled=False
).frame  # BMI and disease progression

1.3.1 What is a model?

Most questions are about how one thing relates to others. Food spending rises with income, but not exactly: two households with the same income spend different amounts. We write that as

$$Y = f(X) + \varepsilon$$

• $Y$ is what we want to understand (food spending),
• $X$ is what we know (income),
• $f$ is the systematic part: what food spending is on average for each income,
• $\varepsilon$ (“epsilon”) is everything else: the noise, the part $X$ cannot explain.

def plot_signal_and_noise(ax: Axes, shown: int = 70) -> None:
    """Food spending against income: a smooth curve for f, and the gaps to it for the noise."""
    x, y = engel.income.to_numpy(), engel.foodexp.to_numpy()
    f = sm.nonparametric.lowess(y, x, frac=0.4, return_sorted=False)
    order = np.argsort(x)
    some = rng.choice(len(x), shown, replace=False)
    ax.vlines(x[some], f[some], y[some], color=style.GREY, lw=0.8, alpha=0.7)
    ax.scatter(x, y, s=10, color=style.BLUE, alpha=0.55)
    ax.plot(x[order], f[order], color=style.CORAL, lw=2.5)

    arrow = dict(arrowstyle="->", color=style.INK)
    at_3000 = (3000, f[np.argmin(abs(x - 3000))])
    ax.annotate(
        "f(X): the average\nfood spending at each income",
        xy=at_3000,
        xytext=(3200, 700),
        color=style.CORAL,
        fontsize=9.5,
        arrowprops=arrow,
    )
    ax.annotate(
        "ε: how far one household\nis from the average",
        xy=(x[some][0], y[some][0]),
        xytext=(250, 1900),
        color=style.INK,
        fontsize=9.5,
        arrowprops=arrow,
    )
    ax.set(xlabel="household income", ylabel="food spending")


fig, ax = plt.subplots(figsize=(8, 3.8))
plot_signal_and_noise(ax)
plt.show()

Figure 1.3.1: The systematic part f (the line) and the noise ε (the gaps to the points).

We never see $f$ directly. We see points, and we estimate $f$ from them. Almost everything in this book is a way to estimate $f$ (with a straight line, a curve, a tree, a network) and a way to say how sure we are.

There are two reasons to estimate $f$. Prediction: we know $X$ for a new case and want a good guess of its $Y$; how $f$ works inside does not matter, as long as the guesses are good. Inference: we want to understand $f$ itself. Which predictors matter? Does food spending rise steadily with income, or level off? The same model can serve both, but the second needs a model we can read, and an honest statement of how sure we are about its parts.

A model is a recipe for $f$ with a few knobs, its parameters, that we set using the data. “A straight line” is a recipe with two knobs, the intercept and the slope; fitting the model means turning the knobs until the line is as close to the points as it gets.

Some recipes fix the shape and only fit the knobs (a line, a parabola). Others let the data decide the shape too. Here are three for the same data:

def plot_line(ax: Axes, x: np.ndarray, y: np.ndarray, grid: np.ndarray) -> None:
    """A straight line: two knobs, and a sentence you can say out loud."""
    slope, intercept = np.polyfit(x, y, 1)
    ax.plot(grid, intercept + slope * grid, color=style.CORAL)
    style.note(ax, f"+{slope:.0f} progression\nper BMI point", y=0.9)
    ax.set_title("Straight line")


def plot_smooth(ax: Axes, x: np.ndarray, y: np.ndarray, grid: np.ndarray) -> None:
    """A smooth curve (lowess): the data choose the shape."""
    curve = sm.nonparametric.lowess(y, x, frac=0.5)
    ax.plot(curve[:, 0], curve[:, 1], color=style.CORAL)
    ax.set_title("Smooth curve")


def plot_every_point(ax: Axes, x: np.ndarray, y: np.ndarray, grid: np.ndarray) -> None:
    """One nearest neighbour: a curve that goes through every training point."""
    model = KNeighborsRegressor(n_neighbors=1).fit(x.reshape(-1, 1), y)
    ax.plot(grid, model.predict(grid.reshape(-1, 1)), color=style.CORAL, lw=1)
    ax.set_title("Follows every point")


x, y = diabetes.bmi.to_numpy(), diabetes.target.to_numpy()
grid = np.linspace(x.min(), x.max(), 400)
fig, ax = style.subplots(3, height=3)
for a, plot in zip(ax, [plot_line, plot_smooth, plot_every_point]):
    a.scatter(x, y, s=6, color=style.BLUE, alpha=0.3)
    plot(a, x, y, grid)
    a.set(xlabel="body mass index", ylabel="disease progression")
fig.tight_layout()
plt.show()

Figure 1.3.2: Three estimates of f for the same data: from rigid to very flexible.

The line gives a sentence you can say out loud. The curve shows that the effect is real but not quite straight. The last one passes through every training point, and would make terrible predictions for new people, because most of its wiggles are noise.

Recipes with a fixed shape (parametric models, like the line) need little data and are easy to read, but are wrong if the true shape is different. Recipes that let the data choose the shape (non-parametric, like the other two) can fit anything, and so can also fit noise.

1.3.2 Model accuracy versus interpretability

Flexible models can be more accurate. Rigid models are easier to read: each extra unit of X adds b to Y is a sentence you can act on. Methods differ in both:

# (method, flexibility, interpretability), placed by hand on a scale from 0 to 10
METHODS = [
    ("Lasso", 0.8, 9.6),
    ("Linear regression", 1.2, 9.2),
    ("Logistic regression", 1.6, 8.6),
    ("Decision trees", 4.6, 7.4),
    ("Generalised additive models", 4.0, 6.6),
    ("Splines", 3.8, 6.0),
    ("K-nearest neighbours", 6.0, 4.0),
    ("Random forests", 7.2, 3.0),
    ("Boosting", 8.1, 2.3),
    ("Support vector machines", 6.8, 2.0),
    ("Neural networks", 9.1, 1.0),
]


def plot_method_map(ax: Axes) -> None:
    """Methods placed by how flexible they are and how easy they are to read."""
    for name, flexibility, interpretability in METHODS:
        ax.scatter(flexibility, interpretability, s=40, color=style.BLUE)
        ax.text(flexibility + 0.15, interpretability + 0.12, name, fontsize=9.5)
    ax.set(xlim=(0, 11.5), ylim=(0, 10.5), xticks=[], yticks=[])
    ax.set(xlabel="flexibility →", ylabel="interpretability →")
    ax.grid(False)


fig, ax = plt.subplots(figsize=(8, 4.2))
plot_method_map(ax)
plt.show()

Figure 1.3.3: Methods in this book, from easy to read and rigid (top left) to flexible and opaque (bottom right).

There is no best place on this map. For inference, stay up and to the left: you need to read the model. For prediction, move right as far as the data allow. The data allow less than you think, because more flexibility needs more data to tell signal from noise.

1.3.3 Supervised models

In supervised learning we have an outcome $Y$ to learn from: every row comes with the right answer, and the model learns to give it. What kind of $Y$ it is decides the kind of problem.

def plot_regression(ax: Axes) -> None:
    """Regression: predict a number (disease progression from body mass index)."""
    x, y = diabetes.bmi.to_numpy(), diabetes.target.to_numpy()
    slope, intercept = np.polyfit(x, y, 1)
    ends = np.array([x.min(), x.max()])
    ax.scatter(x, y, s=7, color=style.BLUE, alpha=0.4)
    ax.plot(ends, intercept + slope * ends, color=style.CORAL)
    ax.set(title="Regression: a number", xlabel="body mass index", ylabel="disease progression")


def plot_classification(ax: Axes) -> None:
    """Classification: predict a group (the species, from two petal measurements)."""
    petals = iris[["petal length (cm)", "petal width (cm)"]].to_numpy()
    model = LogisticRegression(max_iter=1000).fit(petals, iris.target)
    gx, gy = np.meshgrid(np.linspace(0.5, 7.2, 300), np.linspace(0, 2.7, 200))
    zones = model.predict(np.c_[gx.ravel(), gy.ravel()]).reshape(gx.shape)
    ax.contourf(gx, gy, zones, levels=[-0.5, 0.5, 1.5, 2.5], colors=SPECIES_COLORS, alpha=0.12)
    for k, (name, color) in enumerate(zip(SPECIES, SPECIES_COLORS)):
        group = petals[iris.target == k]
        ax.scatter(group[:, 0], group[:, 1], s=9, color=color, label=name)
    ax.legend(fontsize=7.5, loc="upper left")
    ax.set(
        title="Classification: a group", xlabel="petal length (cm)", ylabel="petal width (cm)"
    )


def trend_and_season(t: np.ndarray, months: pd.Index) -> np.ndarray:
    """The design matrix: a curved trend, plus one column per month of the year."""
    month_columns = pd.get_dummies(months).astype(float).to_numpy()[:, 1:]
    return np.column_stack([np.ones(len(t)), t, t**2, month_columns])


def plot_forecast(ax: Axes, months_ahead: int = 36) -> None:
    """Time series: fit trend and season to monthly CO₂, then carry them forward."""
    monthly = co2.co2.resample("MS").mean().loc["1985":]
    t = np.arange(len(monthly))
    coef, *_ = np.linalg.lstsq(
        trend_and_season(t, monthly.index.month), monthly.to_numpy(), rcond=None
    )
    start = monthly.index[-1] + pd.offsets.MonthBegin()
    future = pd.date_range(start, periods=months_ahead, freq="MS")
    t_future = np.arange(len(monthly), len(monthly) + months_ahead)
    ax.plot(monthly.index, monthly.to_numpy(), color=style.BLUE, lw=1.2)
    ax.plot(
        future,
        trend_and_season(t_future, future.month) @ coef,
        color=style.CORAL,
        lw=1.5,
        ls="--",
    )
    ax.set(title="Time series: what comes next", xlabel="year", ylabel="CO₂ (ppm)")


fig, ax = style.subplots(3, height=3.1, width=10)
for a, plot in zip(ax, [plot_regression, plot_classification, plot_forecast]):
    plot(a)
fig.tight_layout()
plt.show()

Figure 1.3.4: Three kinds of supervised problems, each on real data.

Kind	What $Y$ is	Typical question	Methods
Regression	a number	How much? How many?	linear regression, splines, trees
Classification	a group	Which one? Yes or no?	logistic regression, trees, SVMs
Time series	a number, in order	What comes next? Is it changing?	trend + season, ARIMA

Part 4 (regression) and Part 6 (supervised learning) are about these models; Part 9 covers time series.

1.3.4 Unsupervised models

Sometimes there is no $Y$: no right answers to learn from. Unsupervised methods look for structure in the variables themselves: groups of similar rows (clustering), or a few directions that capture most of the variation (principal components).

def plot_clusters(ax: Axes) -> None:
    """Clustering: k-means finds three groups in the petal measurements, with no labels given."""
    petals = iris[["petal length (cm)", "petal width (cm)"]].to_numpy()
    groups = KMeans(3, n_init=10, random_state=0).fit_predict(petals)
    colors = np.array([style.GOLD, style.PURPLE, style.GREY])[groups]
    ax.scatter(petals[:, 0], petals[:, 1], s=10, c=colors)
    ax.set(
        title="Clustering: groups nobody named",
        xlabel="petal length (cm)",
        ylabel="petal width (cm)",
    )


def plot_components(ax: Axes) -> None:
    """Principal components: all four measurements squeezed into two directions."""
    measurements = StandardScaler().fit_transform(iris.drop(columns="target"))
    pca = PCA(2).fit(measurements)
    z = pca.transform(measurements)
    for k, (name, color) in enumerate(zip(SPECIES, SPECIES_COLORS)):
        group = z[iris.target == k]
        ax.scatter(group[:, 0], group[:, 1], s=10, color=color, label=name)
    first, second = pca.explained_variance_ratio_
    ax.legend(fontsize=8)
    ax.set(
        title="Principal components: 4 measurements → 2",
        xlabel=f"component 1 ({first:.0%} of the variation)",
        ylabel=f"component 2 ({second:.0%})",
    )


fig, (left, right) = style.subplots(2, height=3.3)
plot_clusters(left)
plot_components(right)
fig.tight_layout()
plt.show()

Figure 1.3.5: No labels given. Clustering finds groups in the petal measurements (left); principal components squeeze all four measurements into two (right). The species, never shown to the methods, line up with what they found.

Unsupervised results are harder to judge, because there is no right answer to compare with. They are often a first step: find the groups, then study them. Part 7 covers unsupervised learning. Generative models (Chapter 8.4), which learn to produce new data like the old, are unsupervised too.

1.3.5 The bias and variance trade-off

How flexible should a model be? Too rigid, and it misses the shape: its error is bias. Too flexible, and it follows the noise of this particular sample: fit it to another sample and you get a different curve. That is variance. To see both we need a case where we know the truth, so we simulate one: a known curve, plus noise, sampled 30 times.

NOISE = 0.35  # spread of the noise around the truth
N_TRAIN = 60  # points in each training sample
GRID = np.linspace(0, 1, 200)


def truth(x: np.ndarray) -> np.ndarray:
    """The true f, which in real life we never get to see."""
    return np.sin(2.2 * np.pi * x) * 0.9 + x


def training_sample(size: int = N_TRAIN) -> tuple[np.ndarray, np.ndarray]:
    """One training sample: random x values, and the truth plus noise."""
    x = rng.uniform(0, 1, size)
    return x, truth(x) + rng.normal(0, NOISE, size)


def knn_fits(k: int, repeats: int = 30) -> np.ndarray:
    """Fit k-nearest neighbours to `repeats` fresh samples: one row of predictions per fit."""
    fits = []
    for _ in range(repeats):
        x, y = training_sample()
        model = KNeighborsRegressor(n_neighbors=k).fit(x[:, None], y)
        fits.append(model.predict(GRID[:, None]))
    return np.array(fits)


def plot_fits(ax: Axes, k: int, title: str) -> None:
    """Thirty fits of the same model, each on its own sample, against the truth."""
    for curve in knn_fits(k):
        ax.plot(GRID, curve, color=style.BLUE, lw=0.8, alpha=0.35)
    ax.plot(GRID, truth(GRID), color=style.GREY, lw=2, ls="--")
    ax.set(title=title, ylim=(-1.2, 2.4), xlabel="x", ylabel="y")


def plot_trade_off(ax: Axes, max_k: int = 40) -> None:
    """The error on new data split into bias², variance and noise, from rigid to flexible."""
    ks = np.arange(1, max_k + 1)
    bias2, variance = [], []
    for k in ks:
        fits = knn_fits(k, repeats=60)
        bias2.append(((fits.mean(axis=0) - truth(GRID)) ** 2).mean())
        variance.append(fits.var(axis=0).mean())
    bias2, variance = np.array(bias2), np.array(variance)
    flexibility = 1 / ks
    ax.plot(flexibility, bias2, color=style.GOLD, label="bias²")
    ax.plot(flexibility, variance, color=style.PURPLE, label="variance")
    ax.plot(
        flexibility,
        bias2 + variance + NOISE**2,
        color=style.CORAL,
        lw=2.5,
        label="error on new data",
    )
    ax.axhline(NOISE**2, color=style.GREY, ls=":", lw=1)
    ax.text(
        0.98, NOISE**2 + 0.005, "noise", ha="right", va="bottom", color=style.GREY, fontsize=8.5
    )
    ax.set_xscale("log")
    ax.set(
        xticks=[], xlabel="flexibility →", ylabel="mean squared error", title="The trade-off"
    )
    ax.minorticks_off()
    ax.legend(fontsize=8)


fig, ax = style.subplots(3, height=3.2, width=10.5)
plot_fits(ax[0], k=30, title="Rigid: 30 neighbours")
plot_fits(ax[1], k=2, title="Flexible: 2 neighbours")
plot_trade_off(ax[2])
fig.tight_layout()
plt.show()

Figure 1.3.6: The same model fitted to 30 different samples from the same truth (grey dashes). Rigid (left): the fits agree with each other but miss the shape: that is bias. Flexible (middle): right on average, but every sample gives a different curve: that is variance. Right: the error on new data is lowest in between.

Error is measured as the mean squared error: the average squared gap between what the model predicts, $\hat f(x)$, and what happens:

$$\text{MSE} = \frac{1}{n} \sum_{i=1}^{n} \big(y_i - \hat f(x_i)\big)^2$$

The hat on $\hat f$ marks an estimate. For a new case at $x$, averaged over all the training samples we might have had, the expected squared error splits exactly into three parts:

$$E\big(y - \hat f(x)\big)^2 = \underbrace{\text{Var}\big(\hat f(x)\big)}_{\text{variance}} + \underbrace{\big(E \hat f(x) - f(x)\big)^2}_{\text{bias}^2} + \underbrace{\text{Var}(\varepsilon)}_{\text{noise}}$$

Making a model more flexible lowers the bias and raises the variance, so the total makes a U. The noise part never goes away: no model can predict what is truly random.

The catch: the error we can easily measure, on the data we fitted, keeps falling as the model gets more flexible. It cannot see the variance. Move the slider and compare:

That is why a model is always judged on data it has not seen: hold some data back, or use cross-validation (Chapter 6.3). Training error rewards memorising; only new data rewards learning.

1.3.6 Classification: the best possible classifier

For a group outcome the error is the error rate, the share of cases put in the wrong class. If we knew, for every $x$, the probability of each class, the best possible rule would be to pick the most probable class. That rule is the Bayes classifier; its boundary is where the two classes are equally likely, and its error, the Bayes error rate, is the smallest any classifier can reach, like the noise in regression.

Real classifiers must estimate those probabilities from data. K nearest neighbours does it the simplest way: look at the K training cases closest to $x$ and take the share of each class among them. A simulation, where the true probabilities are known, shows how K plays the role of flexibility:

def chance_orange(x1: np.ndarray, x2: np.ndarray) -> np.ndarray:
    """The true probability of the orange class at each point (known only in a simulation)."""
    return 1 / (1 + np.exp(-(3 * np.sin(2.2 * x1) + 2.2 * x2 - 0.8)))


def draw_cases(n: int) -> tuple[np.ndarray, np.ndarray]:
    """Random points with classes drawn from the true probabilities."""
    x = rng.uniform(-2, 2, (n, 2))
    return x, (rng.uniform(size=n) < chance_orange(x[:, 0], x[:, 1])).astype(int)


X_train, y_train = draw_cases(200)
X_test, y_test = draw_cases(5_000)
gx, gy = np.meshgrid(np.linspace(-2, 2, 250), np.linspace(-2, 2, 250))
CLASS_COLORS = np.array([style.BLUE, style.GOLD])


def plot_knn(ax: Axes, k: int) -> None:
    """The K-nearest-neighbour boundary against the Bayes boundary, with the training cases."""
    model = KNeighborsClassifier(k).fit(X_train, y_train)
    zones = model.predict(np.c_[gx.ravel(), gy.ravel()]).reshape(gx.shape)
    ax.contourf(gx, gy, zones, levels=[-0.5, 0.5, 1.5], colors=list(CLASS_COLORS), alpha=0.12)
    ax.contour(gx, gy, zones, levels=[0.5], colors=style.INK, linewidths=1.2)
    ax.contour(
        gx,
        gy,
        chance_orange(gx, gy),
        levels=[0.5],
        colors=style.CORAL,
        linewidths=1.6,
        linestyles="--",
    )
    ax.scatter(X_train[:, 0], X_train[:, 1], c=CLASS_COLORS[y_train], s=10)
    error = (model.predict(X_test) != y_test).mean()
    ax.set(title=f"K = {k}: {error:.0%} wrong", xticks=[], yticks=[], aspect="equal")


fig, axes = style.subplots(3, height=3.4, width=10.5)
for ax, k in zip(axes, (1, 15, 100)):
    plot_knn(ax, k)
fig.tight_layout()
plt.show()

Figure 1.3.7: Simulated data with two classes (blue and orange) whose true probabilities are known. The dashed line is the Bayes boundary, the best possible. Left: K = 1 follows every training point. Middle: K = 15 is close to the Bayes boundary. Right: K = 100 is too smooth.

ks = np.unique(np.geomspace(1, 200, 30).astype(int))
bayes_error = np.mean(
    np.minimum(
        chance_orange(X_test[:, 0], X_test[:, 1]), 1 - chance_orange(X_test[:, 0], X_test[:, 1])
    )
)


def plot_error_curves(ax: Axes) -> None:
    """Training and test error rates against flexibility (1/K)."""
    train, test = [], []
    for k in ks:
        model = KNeighborsClassifier(k).fit(X_train, y_train)
        train.append((model.predict(X_train) != y_train).mean())
        test.append((model.predict(X_test) != y_test).mean())
    ax.plot(1 / ks, train, color=style.GREY, label="training")
    ax.plot(1 / ks, test, color=style.CORAL, label="new cases")
    ax.axhline(bayes_error, color=style.INK, ls="--", lw=1.2, label="Bayes error rate")
    ax.set_xscale("log")
    ax.set_xticks(1 / np.array([100, 30, 10, 3, 1]), ["100", "30", "10", "3", "1"])
    ax.minorticks_off()
    ax.yaxis.set_major_formatter(style.PERCENT)
    ax.legend(fontsize=8)
    ax.set(xlabel="neighbours K (flexible to the right)", ylabel="error rate")


fig, ax = plt.subplots(figsize=(8, 3.3))
plot_error_curves(ax)
plt.show()

Figure 1.3.8: Error rates of K nearest neighbours on the training cases and on 5,000 new cases, as K goes from 200 (rigid) to 1 (flexible). The dashed line is the Bayes error rate. Training error falls to zero; test error makes the familiar U and never goes below the Bayes rate.

The picture is the same as for regression: too rigid misses the shape, too flexible follows the noise, and only error on new data shows where the balance lies.

1.3.7 In one picture

• A model is a recipe for $f$ with knobs we fit to data: $Y = f(X) + \varepsilon$.
• Parametric models fix the shape; non-parametric ones let the data choose it.
• Flexible models can be more accurate, rigid ones are easier to read.
• Supervised: learn from right answers (numbers, groups, series). Unsupervised: find structure with no answers given.
• Error on new data = bias² + variance + noise. Training error always looks better than it is.
• Prediction needs good guesses; inference needs a model we can read.
• The Bayes classifier picks the most probable class; no classifier can beat its error rate. K nearest neighbours estimates it, with K as the flexibility knob.

Next: graphs that show data honestly, and the ways they can mislead.