6.5 Beyond straight lines

A curve can be built from simple pieces: powers, steps or short polynomials joined smoothly at knots. The fitting is still least squares; only the columns change.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from matplotlib.axes import Axes
from patsy import dmatrix
from scipy.interpolate import make_smoothing_spline

from statly import style

style.use()

co2 = sm.datasets.co2.load_pandas().data.interpolate()  # weekly CO2 at Mauna Loa, 1958 to 2001
part = co2.loc["1990":"1995"].co2
t = ((part.index - part.index[0]).days / 365.25).to_numpy()  # years since the start of 1990
ppm = part.to_numpy()
GRID = np.linspace(t.min(), t.max(), 600)

star = sm.datasets.star98.load_pandas().data  # 303 California school districts, 1998
star["passed"] = star.NABOVE / (
    star.NABOVE + star.NBELOW
)  # share above the national median in maths


def plot_data(ax: Axes, title: str) -> None:
    """The weekly CO₂ values as small grey dots."""
    ax.scatter(1990 + t, ppm, s=4, color=style.GREY, alpha=0.6)
    ax.set(title=title, xlabel="year", ylabel="CO₂ (ppm)")

6.5.1 Polynomials

Six years of weekly CO₂ at Mauna Loa: a rise plus a yearly wave. No straight line will do. The oldest fix is to add powers of x: x², x³ and so on. The model is still linear in its coefficients, so least squares fits it as before; only the columns are new.

WIDE = np.linspace(t.min() - 0.25, t.max() + 0.25, 700)


def plot_polynomial(ax: Axes, degree: int) -> None:
    """The data with a polynomial of the given degree, drawn a little beyond both ends."""
    plot_data(ax, f"Degree {degree}")
    for start, end in ((WIDE[0], t.min()), (t.max(), WIDE[-1])):
        ax.axvspan(1990 + start, 1990 + end, color=style.GREY, alpha=0.15, lw=0)
    curve = np.polynomial.Polynomial.fit(t, ppm, degree)
    ax.plot(1990 + WIDE, curve(WIDE), color=style.CORAL, lw=1.8)
    ax.set_ylim(ppm.min() - 4, ppm.max() + 4)


fig, axes = style.subplots(3, height=2.9, width=10.5, sharey=True)
for ax, degree in zip(axes, (3, 15, 25)):
    plot_polynomial(ax, degree)
fig.tight_layout()
plt.show()

Figure 6.5.1: Polynomials of increasing degree fitted to six years of weekly CO₂, drawn a little beyond the data (shaded). Degree 3 cannot bend often enough; degree 15 bends in the wrong places; degree 25 follows the wave, but just past the ends of the data it shoots off.

Polynomials are global: every coefficient affects the whole curve, so to bend in one place they must bend everywhere. It takes 26 coefficients to follow six waves, and at the ends, where no points hold it down, the curve shoots off. Use polynomials for gentle curves, a square or a cube, and not for wiggly shapes or for predicting beyond the data.

6.5.2 Steps and splines

The opposite idea is to work locally. Cut x into pieces and fit each piece on its own. With a constant in each piece, the result is a step function. With a cubic polynomial in each piece, joined so that the curve and its slope and bend match at the joins, the result is a regression spline. The joins are called knots.

KNOTS = np.arange(0.25, 6, 0.25)  # a knot every quarter year
blocks = dmatrix(
    "bs(t, knots=KNOTS, degree=3, include_intercept=True) - 1",
    {"t": t},
    return_type="dataframe",
)
grid_blocks = dmatrix(blocks.design_info, {"t": GRID})
spline = sm.OLS(ppm, blocks).fit()


def plot_steps(ax: Axes) -> None:
    """The average in each quarter year, drawn as steps."""
    plot_data(ax, "Steps")
    quarter = np.floor(t * 4) / 4
    means = pd.Series(ppm).groupby(quarter).mean()
    ax.step(1990 + means.index, means.to_numpy(), where="post", color=style.CORAL, lw=1.8)


def plot_blocks(ax: Axes) -> None:
    """Every spline building block over the six years."""
    ax.plot(1990 + GRID, np.asarray(grid_blocks), lw=1.2)
    for knot in KNOTS:
        ax.axvline(1990 + knot, color=style.INK, lw=0.5, alpha=0.4)
    ax.set(title="Building blocks", xlabel="year", ylabel="block value")


def plot_spline(ax: Axes) -> None:
    """The fitted spline, with the knots marked."""
    plot_data(ax, "Spline")
    ax.plot(
        1990 + GRID,
        np.asarray(grid_blocks) @ spline.params.to_numpy(),
        color=style.CORAL,
        lw=1.8,
    )
    for knot in KNOTS:
        ax.axvline(1990 + knot, color=style.INK, lw=0.5, alpha=0.4)


fig, axes = style.subplots(3, height=2.9, width=10.5)
for ax, plot in zip(axes, [plot_steps, plot_blocks, plot_spline]):
    plot(ax)
fig.tight_layout()
plt.show()

Figure 6.5.2: Left: a step function with one step per quarter year. Middle: the building blocks of a cubic spline with knots every quarter year (each block is one column of the model). Right: the fitted spline, a weighted sum of the blocks, follows the wave without swinging at the ends.

Each block is a little hump that lives only between a few knots, so a change in one place stays in that place. The model is still a regression: the blocks are its columns, least squares finds their weights, and everything from Part 4, intervals and tests included, works unchanged.

The knots set the flexibility. More knots, more wiggle; too many, and the spline chases noise. The number of knots is a knob for cross-validation (Chapter 6.3):

DFS = np.r_[1, np.arange(3, 16)]  # a straight line, then cubic splines (which need 3 or more)
folds = np.array_split(np.random.default_rng(1).permutation(len(star)), 10)


def cv_spline(df: int) -> float:
    """10-fold cross-validated RMSE of a spline (or a line, for df = 1) in LOWINC."""
    formula = (
        "passed ~ LOWINC"
        if df == 1
        else f"passed ~ bs(LOWINC, df={df}, lower_bound=0, upper_bound=100)"
    )
    errors = []
    for held in folds:
        train = star.drop(index=star.index[held])
        fit = smf.ols(formula, train).fit()
        errors.append(np.mean((star.passed.iloc[held] - fit.predict(star.iloc[held])) ** 2))
    return float(np.sqrt(np.mean(errors)))


def plot_spline_cv(ax: Axes) -> None:
    """Cross-validated error against degrees of freedom."""
    ax.plot(DFS, [cv_spline(df) for df in DFS], "o-", color=style.BLUE, ms=4)
    ax.set(
        title="Cross-validation",
        xlabel="degrees of freedom",
        ylabel="error (RMSE)",
        xticks=[1, 3, 5, 7, 9, 11, 13, 15],
    )


def plot_spline_fits(ax: Axes) -> None:
    """The data with three spline fits of increasing flexibility."""
    ax.scatter(star.LOWINC, star.passed, s=6, color=style.GREY, alpha=0.5)
    grid = pd.DataFrame({"LOWINC": np.linspace(star.LOWINC.min(), star.LOWINC.max(), 300)})
    for df, color in ((1, style.GOLD), (4, style.CORAL), (15, style.PURPLE)):
        formula = (
            "passed ~ LOWINC"
            if df == 1
            else f"passed ~ bs(LOWINC, df={df}, lower_bound=0, upper_bound=100)"
        )
        ax.plot(
            grid.LOWINC,
            smf.ols(formula, star).fit().predict(grid),
            color=color,
            lw=1.6,
            label=f"df = {df}",
        )
    ax.yaxis.set_major_formatter(style.PERCENT)
    ax.legend(fontsize=8)
    ax.set(title="Three fits", xlabel="low-income pupils (%)", ylabel="share above the median")


fig, axes = style.subplots(2, height=3.1, width=9.5)
plot_spline_cv(axes[0])
plot_spline_fits(axes[1])
fig.tight_layout()
plt.show()

Figure 6.5.3: Choosing a spline’s flexibility for the share of pupils above the national median in maths against the share of low-income pupils (303 California districts). Left: 10-fold cross-validated error by degrees of freedom (number of spline columns); 3 or 4 are as good as anything. Right: the fits with 1 (a straight line), 4 and 15 degrees of freedom.

Cubic splines have one more weakness: beyond the outer knots they are still cubics, free to swing. A natural spline adds the constraint that the curve is a straight line beyond the boundary knots. That costs nothing inside the data and makes the edges, and any careful extrapolation, much calmer:

EDGE = np.linspace(t.min() - 0.5, t.max() + 0.5, 700)
EDGE_KNOTS = np.quantile(
    t, np.linspace(0, 1, 24)
)  # 24 knots, the outer two at the ends of the data


def cubic_basis(x: np.ndarray) -> np.ndarray:
    """Cubic spline columns: 1, x, x², x³ and (x − knot)³ for x past each inner knot."""
    powers = [np.ones_like(x), x, x**2, x**3]
    return np.column_stack(powers + [np.clip(x - k, 0, None) ** 3 for k in EDGE_KNOTS[1:-1]])


def natural_basis(x: np.ndarray) -> np.ndarray:
    """Natural cubic spline columns: cubic between the knots, straight beyond the outer ones."""
    last = EDGE_KNOTS[-1]

    def d(k: float) -> np.ndarray:
        return (np.clip(x - k, 0, None) ** 3 - np.clip(x - last, 0, None) ** 3) / (last - k)

    return np.column_stack(
        [np.ones_like(x), x] + [d(k) - d(EDGE_KNOTS[-2]) for k in EDGE_KNOTS[:-2]]
    )


cubic_coef = np.linalg.lstsq(cubic_basis(t), ppm, rcond=None)[0]
natural_coef = np.linalg.lstsq(natural_basis(t), ppm, rcond=None)[0]


def plot_edges(ax: Axes) -> None:
    """Both splines over the data and a little beyond it."""
    plot_data(ax, "Cubic and natural splines")
    for start, end in ((EDGE[0], t.min()), (t.max(), EDGE[-1])):
        ax.axvspan(1990 + start, 1990 + end, color=style.GREY, alpha=0.15, lw=0)
    ax.plot(
        1990 + EDGE,
        cubic_basis(EDGE) @ cubic_coef,
        color=style.BLUE,
        lw=1.5,
        label="cubic spline",
    )
    ax.plot(
        1990 + EDGE,
        natural_basis(EDGE) @ natural_coef,
        color=style.CORAL,
        lw=1.5,
        label="natural spline",
    )
    ax.set_ylim(ppm.min() - 5, ppm.max() + 5)
    ax.legend(fontsize=8, loc="upper left")


fig, ax = plt.subplots(figsize=(8, 3.2))
plot_edges(ax)
plt.show()

Figure 6.5.4: The six years of CO₂ fitted with a cubic spline (blue) and a natural spline (red) with the same number of columns, drawn half a year beyond the data on each side (shaded). Inside, they agree. Outside, the cubic spline bends away; the natural spline continues in a straight line.

6.5.3 Smoothing splines

A smoothing spline avoids choosing knots. It puts a knot at every data point, which would fit every point exactly, and then adds a penalty for bending, like the ridge penalty of Chapter 5.3 but on the curve’s wiggliness. One knob, λ, sets the balance between following the points and staying smooth.

def plot_smoothing(ax: Axes, lam: float, title: str) -> None:
    """A smoothing spline with penalty lam."""
    plot_data(ax, title)
    ax.plot(
        1990 + GRID, make_smoothing_spline(t, ppm, lam=lam)(GRID), color=style.CORAL, lw=1.6
    )


fig, axes = style.subplots(3, height=2.9, width=10.5, sharey=True)
for ax, (lam, title) in zip(
    axes, ((1e-7, "Tiny penalty"), (1e-3, "Middle penalty"), (10, "Huge penalty"))
):
    plot_smoothing(ax, lam, title)
fig.tight_layout()
plt.show()

Figure 6.5.5: Smoothing splines with three penalties. A tiny penalty follows the weekly noise, a huge one flattens the wave into the trend, a middle one catches the wave and ignores the noise.

In symbols, the smoothing spline minimises

$$\sum_i \big(y_i - g(x_i)\big)^2 + \lambda \int g''(x)^2\, dx$$

where $g''$ is the curve’s second derivative, its bending. λ = 0 interpolates every point; λ → ∞ leaves only the straight line, which never bends. The amount of flexibility a λ allows is usually reported as effective degrees of freedom, comparable with the number of columns of a regression spline.

6.5.4 Local regression

Local regression (loess, or lowess, as used in many figures of this book) takes a different route. To estimate the curve at one point $x_0$, it takes the share s of the data closest to $x_0$ (the span), weights them so that the nearest count most (usually with the tricube weight $(1 - d^3)^3$, where d is the distance scaled to the window), and fits a weighted straight line through them. The curve’s height at $x_0$ is that line’s height there. Then it moves to the next point and starts again.

x_loc, y_loc = star.LOWINC.to_numpy(), star.passed.to_numpy()
X0, SPAN = 60, 0.3


def local_weights(x0: float, span: float) -> np.ndarray:
    """Tricube weights for the share `span` of points nearest to x0 (zero for the rest)."""
    distance = abs(x_loc - x0)
    radius = np.sort(distance)[int(span * len(x_loc)) - 1]
    scaled = np.clip(distance / radius, 0, 1)
    return (1 - scaled**3) ** 3


def plot_one_point(ax: Axes) -> None:
    """The local window, its weights and the weighted line at one point."""
    w = local_weights(X0, SPAN)
    inside = w > 0
    ax.scatter(x_loc[~inside], y_loc[~inside], s=5, color=style.GREY, alpha=0.3)
    ax.scatter(x_loc[inside], y_loc[inside], s=4 + 60 * w[inside], color=style.BLUE, alpha=0.6)
    slope, intercept = np.polyfit(x_loc[inside], y_loc[inside], 1, w=np.sqrt(w[inside]))
    window = np.array([x_loc[inside].min(), x_loc[inside].max()])
    ax.plot(window, intercept + slope * window, color=style.CORAL, lw=2)
    ax.plot(X0, intercept + slope * X0, "o", ms=9, color=style.CORAL)
    ax.axvline(X0, color=style.INK, ls="--", lw=1)
    ax.yaxis.set_major_formatter(style.PERCENT)
    ax.set(title="One point", xlabel="low-income pupils (%)", ylabel="share above the median")


def plot_spans(ax: Axes) -> None:
    """Loess curves for three spans."""
    ax.scatter(x_loc, y_loc, s=5, color=style.GREY, alpha=0.4)
    for span, color in ((0.1, style.PURPLE), (0.3, style.CORAL), (0.8, style.GOLD)):
        curve = sm.nonparametric.lowess(y_loc, x_loc, frac=span)
        ax.plot(curve[:, 0], curve[:, 1], color=color, lw=1.8, label=f"span {span:.0%}")
    ax.yaxis.set_major_formatter(style.PERCENT)
    ax.legend(fontsize=8)
    ax.set(title="The whole curve", xlabel="low-income pupils (%)")


fig, axes = style.subplots(2, height=3.2, width=9.5, sharey=True)
plot_one_point(axes[0])
plot_spans(axes[1])
fig.tight_layout()
plt.show()

Figure 6.5.6: Local regression of the share of pupils above the median on the share of low-income pupils. Left: at one point (60% low-income pupils, dashed), the nearest 30% of districts get weights (dot size) that fall with distance, and a weighted line (red) is fitted through them; its height at 60% is one point of the curve. Right: repeating this at every point traces the curve, shown for spans of 10%, 30% and 80%.

The span is the flexibility knob: small spans follow every bump, large ones approach a straight line. Local regression needs no knots and adapts well, but it gives no simple formula, and in many dimensions the “nearby” points become too few (the curse of dimensionality).

6.5.5 Several curves at once: additive models

Real questions have several predictors, each with its own curve. A generalised additive model (GAM) adds one smooth curve per predictor:

$$Y = a + f_1(X_1) + f_2(X_2) + \dots + \varepsilon$$

Each $f$ is a spline, and the whole model is fitted at once. It keeps what makes linear regression readable, one effect per predictor with the others held fixed, but lets each effect bend.

In 303 California school districts, what goes with the share of pupils scoring above the national median in maths?

PREDICTORS = {
    "LOWINC": "low-income pupils (%)",
    "PCTAF": "in university-prep courses (%)",
    "PERMINTE": "minority teachers (%)",
    "AVYRSEXP": "teacher experience (years)",
}
formula = "passed ~ " + " + ".join(f"bs({name}, df=5)" for name in PREDICTORS)
gam = smf.ols(formula, star).fit()


def plot_effect(ax: Axes, name: str, label: str) -> None:
    """One predictor's curve, the others held at their median, with its band."""
    grid = pd.DataFrame({other: np.full(100, star[other].median()) for other in PREDICTORS})
    grid[name] = np.linspace(star[name].quantile(0.01), star[name].quantile(0.99), 100)
    band = gam.get_prediction(grid).summary_frame(alpha=0.05)
    ax.fill_between(
        grid[name], band.mean_ci_lower, band.mean_ci_upper, color=style.BLUE, alpha=0.2, lw=0
    )
    ax.plot(grid[name], band["mean"], color=style.BLUE)
    ax.plot(star[name], np.full(len(star), 0.05), "|", ms=5, color=style.INK, alpha=0.4)
    ax.yaxis.set_major_formatter(style.PERCENT)
    ax.set(xlabel=label, ylim=(0, 1))


fig, axes = style.subplots(4, cols=4, height=2.9, width=11, sharey=True)
for ax, (name, label) in zip(axes, PREDICTORS.items()):
    plot_effect(ax, name, label)
axes[0].set_ylabel("share above the median")
fig.tight_layout()
plt.show()

Figure 6.5.7: An additive model for the share of pupils above the national median in maths, in 303 California school districts. Each panel shows one predictor’s curve, with the others held at their median, and a 95% band. Low income matters most and bends: the drop is steepest between 0 and 40%.

Each panel reads like a slope in Chapter 4.2, but bent. The share of low-income pupils dominates, and its effect is strongest at the low end. Above about 15%, more pupils in university-prep courses go with higher scores. The other two add little once these are known. The ticks along the bottom show where the districts are; the bands widen where they are few.

GAMs are a good default when you suspect curves and still want to read the model. What they cannot do on their own is interactions, where the curve of one predictor depends on another; trees, in the next chapter, find those automatically.

6.5.6 In one picture

• Polynomials add powers of x; they are global and swing at the ends.
• Step functions and splines work locally; splines join cubic pieces smoothly at knots.
• Natural splines are straight beyond the boundary knots, which calms the edges.
• Smoothing splines put a knot everywhere and penalise bending; one knob sets the smoothness.
• Local regression fits a weighted line in a moving window; the span sets the smoothness.
• Choose knots, λ or span by cross-validation.
• All of these are still least squares on new columns, so intervals and tests carry over.
• A GAM gives every predictor its own curve, added together.

Next: trees, which cut the data into boxes.