4.2 Multiple regression

With several predictors, each slope answers a new question: how does the outcome change with this predictor when the others are held fixed? That can change the answer completely.

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 scipy import stats
from sklearn import datasets
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.regression.linear_model import RegressionResults

from statly import style

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

states = sm.datasets.statecrime.load_pandas().data.drop(
    "District of Columbia"
)  # crime, poverty and families in the 50 US states
iris = datasets.load_iris(as_frame=True).frame.set_axis(
    ["sepal_length", "sepal_width", "petal_length", "petal_width", "species"], axis=1
)
iris["species"] = iris.species.map(dict(enumerate(SPECIES)))
diabetes = datasets.load_diabetes(
    as_frame=True, scaled=False
).frame  # 442 diabetes patients: blood tests and disease progression a year later

4.2.1 Holding the others fixed

Across the 50 US states, violent crime is higher where more people are poor: each extra percentage point of poverty goes with about 16 more violent crimes per 100,000 people. But poor states also have more single-parent households. Which one goes with crime?

Multiple regression puts both in one model:

$$\text{crime} = a + b_1 \cdot \text{poverty} + b_2 \cdot \text{single parents} + \varepsilon$$

Now $b_1$ answers a different question: among states with the same share of single parents, how does crime change with poverty?

alone = smf.ols("violent ~ poverty", states).fit()
both = smf.ols("violent ~ poverty + single", states).fit()


def plot_states(ax: Axes) -> None:
    """Violent crime against poverty, coloured by the share of single parents."""
    points = ax.scatter(
        states.poverty,
        states.violent,
        c=states.single,
        cmap="Blues",
        s=30,
        edgecolors=style.INK,
        linewidths=0.4,
    )
    grid = np.array([8, 23])
    ax.plot(grid, alone.params.Intercept + alone.params.poverty * grid, color=style.CORAL)
    plt.colorbar(points, ax=ax, label="single parents (%)")
    ax.set(title="States", xlabel="poverty (%)", ylabel="violent crimes per 100,000")


def plot_slopes(ax: Axes) -> None:
    """The poverty slope with its interval, in the two models."""
    for i, (name, model) in enumerate(
        (("poverty alone", alone), ("with single parents", both))
    ):
        low, high = model.conf_int().loc["poverty"]
        ax.plot([low, high], [i, i], color=style.BLUE, lw=2.5)
        ax.plot(model.params.poverty, i, "o", ms=7, color=style.BLUE)
    ax.axvline(0, color=style.INK, lw=1)
    ax.set_yticks([0, 1], ["poverty alone", "with single parents"])
    ax.invert_yaxis()
    ax.set(
        title="Slope of poverty",
        xlabel="more violent crimes per point of poverty",
        ylim=(1.6, -0.6),
    )


fig, axes = style.subplots(2, height=3.3, width=10)
plot_states(axes[0])
plot_slopes(axes[1])
fig.tight_layout()
plt.show()

Figure 4.2.1: Left: violent crime and poverty in the 50 states, coloured by the share of single-parent households (light: few, dark: many), with the line for poverty alone. Right: the slope of poverty, with its 95% interval, alone and with single parents in the model. Held fixed, the single-parent share takes over and the poverty slope drops to about zero.

With both in the model, poverty’s slope drops from 16 to about −1, with an interval that easily includes 0, while each extra point of single parents goes with about 31 more violent crimes. Poverty and single parenthood travel together across states, and once single parenthood is accounted for, poverty adds little.

Careful with the conclusion. This does not show that single parenthood causes crime. Both may stand for something else, such as unemployment or history, that the model does not contain. Multiple regression adjusts only for what is in it.

4.2.2 What “held fixed” means

There is a precise way to see a slope in a multiple regression. To get poverty’s slope with single parenthood held fixed: take out of violent crime everything single parenthood can explain (the residuals of crime on single parenthood), take out of poverty everything single parenthood can explain (the residuals of poverty on single parenthood), and draw a line through what is left. Its slope is exactly the multiple regression slope. This picture is an added-variable plot.

def residual_of(target: str, given: str) -> np.ndarray:
    """What is left of one column after a straight-line fit on another."""
    return smf.ols(f"{target} ~ {given}", states).fit().resid.to_numpy()


def plot_added(ax: Axes, predictor: str, other: str, label: str) -> None:
    """Residual of crime against residual of the predictor, both cleaned of the other predictor."""
    x, y = residual_of(predictor, other), residual_of("violent", other)
    slope = np.polyfit(x, y, 1)[0]
    ax.scatter(x, y, s=14, color=style.PURPLE, alpha=0.7)
    grid = np.array([x.min(), x.max()])
    ax.plot(grid, slope * grid, color=style.CORAL)
    style.note(
        ax,
        f"slope {slope:.1f}\nmodel's slope {both.params[predictor]:.1f}",
        x=0.03,
        y=0.98,
        color=style.INK,
    )
    ax.set(
        title=f"{label}, {other.replace('single', 'single parents')} removed",
        xlabel=f"{label}, residual",
        ylabel="violent crime, residual",
    )


fig, axes = style.subplots(2, height=3.2, width=9.5)
plot_added(axes[0], "poverty", "single", "poverty")
plot_added(axes[1], "single", "poverty", "single parents")
fig.tight_layout()
plt.show()

Figure 4.2.2: Added-variable plots for the model of violent crime on poverty and single parenthood. Left: poverty, after removing what single parenthood explains from both axes: the line through the residuals has exactly the multiple regression slope for poverty (about −1). Right: single parenthood, after removing poverty: a steep slope (about 31).

The added-variable plot also shows which states drive a slope, and whether the relation that is left is straight. It is one of the most useful checks of a multiple regression (Chapter 4.3).

4.2.3 The general model and its table

With p predictors the model is

$$Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \dots + \beta_p X_p + \varepsilon$$

Stack the predictors as columns of a table X (with a first column of ones for the intercept) and the outcomes as a column y. Least squares then has one compact solution, which every program computes:

$$\hat{\boldsymbol\beta} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y}$$

The details of the matrix algebra matter less than the result: one coefficient per predictor, each with a standard error, a t and a p-value, read exactly as in Chapter 4.1, except that each one now means “with the others held fixed”. For the diabetes patients:

PREDICTORS = ["age", "bmi", "bp", "s1", "s2", "s3", "s4", "s5", "s6"]
full = smf.ols("target ~ " + " + ".join(PREDICTORS), diabetes).fit()
table = full.summary2().tables[1]
table.columns = ["estimate", "standard error", "t", "p-value", "95% from", "95% to"]
table = table.round(2).astype(object)
table["p-value"] = ["< 0.001" if p < 0.001 else f"{p:.3f}" for p in full.pvalues]
print(table.to_markdown(colalign=["left"] + ["right"] * 6))

Table 4.2.1: Disease progression a year later, regressed on age, body mass index, blood pressure and the six blood measurements (s1 to s6) of 442 diabetes patients.

	estimate	standard error	t	p-value	95% from	95% to
Intercept	-363.9	68.14	-5.34	< 0.001	-497.83	-229.97
age	-0.12	0.22	-0.55	0.583	-0.55	0.31
bmi	6	0.72	8.32	< 0.001	4.59	7.42
bp	0.95	0.22	4.23	< 0.001	0.51	1.39
s1	-0.98	0.58	-1.68	0.093	-2.12	0.16
s2	0.66	0.54	1.22	0.223	-0.4	1.72
s3	0.51	0.79	0.65	0.518	-1.05	2.08
s4	4.66	6.04	0.77	0.441	-7.21	16.53
s5	68.95	15.93	4.33	< 0.001	37.64	100.25
s6	0.2	0.28	0.73	0.465	-0.34	0.75

Body mass index, blood pressure and s5 are clearly related to progression with the others held fixed. Most blood measurements are not individually significant, yet they are strongly related to each other (Section 4.2.9), which is one reason single p-values can mislead.

4.2.4 Testing several coefficients at once: the F-test

The table tests one coefficient at a time. Two other questions need a test of several at once. Is the model any better than no model at all (are all slopes zero)? And do the six blood measurements, as a group, add anything to age, body mass index and blood pressure?

Both are answered by comparing the residual sum of squares of a smaller model (RSS₀) with a larger one that contains it (RSS₁). The F statistic asks whether the drop in RSS is more than the extra q coefficients would buy by chance:

$$F = \frac{(\text{RSS}_0 - \text{RSS}_1) / q}{\text{RSS}_1 / (n - p - 1)}$$

with q and n − p − 1 degrees of freedom, p being the number of predictors in the larger model. It is the ANOVA of Chapter 3.5 in regression clothes.

small = smf.ols("target ~ age + bmi + bp", diabetes).fit()
comparison = sm.stats.anova_lm(small, full)
f_value, f_p = comparison.F.iloc[1], comparison["Pr(>F)"].iloc[1]


def plot_rss(ax: Axes) -> None:
    """The two residual sums of squares as bars."""
    values = [small.ssr / 1e6, full.ssr / 1e6]
    ax.bar(
        ["age, bmi, bp", "+ six blood\nmeasurements"],
        values,
        color=[style.GREY, style.BLUE],
        alpha=0.7,
    )
    for i, v in enumerate(values):
        ax.text(i, v + 0.01, f"{v:.2f} million", ha="center", fontsize=9)
    ax.set(
        title="What is left",
        ylabel="residual sum of squares (millions)",
        ylim=(0, max(values) * 1.15),
    )


def plot_f_curve(ax: Axes) -> None:
    """The F distribution under 'the six add nothing', with the observed F."""
    x = np.linspace(0.01, 12, 400)
    ax.plot(x, stats.f.pdf(x, 6, full.df_resid), color=style.BLUE)
    ax.axvline(f_value, color=style.CORAL, lw=1.6)
    style.note(ax, f"F = {f_value:.1f}\np = {f_p:.0e}", x=0.6, color=style.INK)
    ax.set(title="If the six added nothing", xlabel="F", ylabel="how common", yticks=[])


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

Figure 4.2.3: The F-test for adding the six blood measurements to a model with age, body mass index and blood pressure. Left: residual sums of squares of the two models. Right: the F distribution with 6 and 432 degrees of freedom, with the observed F far in its tail.

The six together clearly add something, even though most of them are not significant alone: they share their information, and the F-test counts it once. The overall F-test at the bottom of every regression output is the same comparison against the model with only an intercept.

4.2.5 R² and adjusted R²

R² (Chapter 4.1) can only go up when a predictor is added, even a useless one: least squares will always find some use for it in the sample. The adjusted R² charges for each predictor:

$$R^2_{\text{adj}} = 1 - \frac{\text{RSS} / (n - p - 1)}{\text{TSS} / (n - 1)}$$

noise = rng.normal(size=(len(diabetes), 60))
base = diabetes[PREDICTORS].to_numpy()


def r2_pair(extra: int) -> tuple[float, float]:
    """R² and adjusted R² with a number of noise columns added."""
    X = sm.add_constant(np.c_[base, noise[:, :extra]])
    model = sm.OLS(diabetes.target.to_numpy(), X).fit()
    return model.rsquared, model.rsquared_adj


extras = np.arange(0, 61, 3)
pairs = np.array([r2_pair(k) for k in extras])

fig, ax = plt.subplots(figsize=(8, 3.1))
ax.plot(extras, pairs[:, 0], "o-", ms=3, color=style.CORAL, label="R²")
ax.plot(extras, pairs[:, 1], "o-", ms=3, color=style.BLUE, label="adjusted R²")
ax.legend(fontsize=8)
ax.set(xlabel="predictors of pure noise added", ylabel="share of variation explained")
plt.show()

Figure 4.2.4: The diabetes model with 0 to 60 extra predictors of pure noise added one by one. R² creeps up with every one of them; adjusted R² wobbles by chance and goes nowhere. Neither is a substitute for checking on new data (Chapter 6.3).

4.2.6 Confounding and Simpson’s paradox

A third variable that is related to both the predictor and the outcome is a confounder. Leave it out and the slope you see mixes its effect with the one you want. In extreme cases the slope even flips its sign:

overall = smf.ols("sepal_width ~ sepal_length", iris).fit()
with_species = smf.ols("sepal_width ~ sepal_length + species", iris).fit()
GRID = np.linspace(4.2, 8, 50)


def plot_flowers(ax: Axes, by_species: bool) -> None:
    """The flowers, grey or coloured by species."""
    for name, color in zip(SPECIES, SPECIES_COLORS):
        group = iris[iris.species == name]
        ax.scatter(
            group.sepal_length,
            group.sepal_width,
            s=12,
            color=color if by_species else style.GREY,
            alpha=0.6,
            label=name if by_species else None,
        )
    ax.set(xlabel="sepal length (cm)", ylabel="sepal width (cm)")


def plot_overall(ax: Axes) -> None:
    """All flowers with one line."""
    plot_flowers(ax, by_species=False)
    ax.plot(GRID, overall.predict(pd.DataFrame({"sepal_length": GRID})), color=style.INK)
    style.note(ax, f"slope {overall.params.sepal_length:+.2f}", x=0.62, color=style.INK)
    ax.set_title("Ignoring species")


def plot_within(ax: Axes) -> None:
    """Each species with its own line from the model with species."""
    plot_flowers(ax, by_species=True)
    for name, color in zip(SPECIES, SPECIES_COLORS):
        group = iris[iris.species == name]
        grid = np.linspace(group.sepal_length.min(), group.sepal_length.max(), 20)
        ax.plot(
            grid,
            with_species.predict(pd.DataFrame({"sepal_length": grid, "species": name})),
            color=color,
        )
    style.note(ax, f"slope {with_species.params.sepal_length:+.2f}", x=0.62, color=style.INK)
    ax.legend(fontsize=8, loc="lower right")
    ax.set_title("Within species")


fig, axes = style.subplots(2, height=3.3, width=9.5, sharey=True)
plot_overall(axes[0])
plot_within(axes[1])
fig.tight_layout()
plt.show()

Figure 4.2.5: Sepal width against sepal length for the 150 iris flowers. Left: ignoring species, the line falls slightly. Right: within each species, wider goes with longer; species is the confounder. The lines are from one model with species as a predictor.

Setosa flowers have short but wide sepals; the other two species have long but narrow ones. Lumped together, the species differences hide the relation inside each species. This reversal is Simpson’s paradox, and it is not rare: a treatment can look harmful overall because the sickest patients get it most often.

Which slope is right depends on the question. To predict the width of a flower of unknown species, the overall line is fine. To understand how width grows with length, hold species fixed.

4.2.7 Groups as predictors

Species is not a number, so how does it go into a regression? As indicator columns: one column that is 1 for versicolor and 0 otherwise, one for virginica. Setosa, with zeros in both, is the reference group. Each indicator’s coefficient is how far that species’ line sits above or below setosa’s, at the same sepal length. That is what drew the three parallel lines above. A few rows of the table the model actually sees:

flower	intercept	sepal length	versicolor	virginica
a setosa	1	5.1	0	0
a versicolor	1	7.0	1	0
a virginica	1	6.3	0	1

Why not a third column for setosa? It would always equal 1 minus the other two, so it would carry no new information, and the model could not separate it from the intercept. With k groups, use k − 1 indicators. Testing whether species matters at all means testing both indicators together: the F-test of Section 4.2.4.

4.2.8 Interactions

Parallel lines assume the slope is the same in every group. An interaction lets it differ: the effect of one predictor depends on the value of another.

interaction = smf.ols("sepal_width ~ sepal_length * species", iris).fit()


def plot_model_lines(ax: Axes, model: RegressionResults, title: str) -> None:
    """The flowers with one fitted line per species from the given model."""
    plot_flowers(ax, by_species=True)
    for name, color in zip(SPECIES, SPECIES_COLORS):
        group = iris[iris.species == name]
        grid = np.linspace(group.sepal_length.min(), group.sepal_length.max(), 20)
        ax.plot(
            grid,
            model.predict(pd.DataFrame({"sepal_length": grid, "species": name})),
            color=color,
            lw=2.2,
        )
    ax.set_title(title)


fig, axes = style.subplots(2, height=3.3, width=9.5, sharey=True)
plot_model_lines(axes[0], with_species, "Same slope")
plot_model_lines(axes[1], interaction, "Own slope per species")
axes[1].legend(fontsize=8, loc="lower right")
fig.tight_layout()
plt.show()

Figure 4.2.6: Left: one slope for all species, so the lines are parallel. Right: with an interaction, each species gets its own slope; setosa’s is much steeper.

In the model, the interaction adds columns that are sepal length times each indicator. Their coefficients are how much each species’ slope differs from setosa’s: setosa gains 0.80 cm of width per cm of length, versicolor about 0.32 and virginica about 0.23. With an interaction in the model, never read a main slope alone: it is the slope for the reference group only.

Curves enter a regression the same way, as extra columns. Adding $x^2$ lets the line bend: the model $y = \beta_0 + \beta_1 x + \beta_2 x^2$ is still linear in its coefficients, so everything in this chapter applies (Chapter 6.5 goes further). The hierarchy principle says: if a model contains $x^2$ or an interaction $x \cdot z$, keep $x$ (and $z$) in it too, even if their own p-values are large; without them, the curve or the interaction is forced through an arbitrary point.

4.2.9 When predictors carry the same information

The diabetes data have two cholesterol measures, total (s1) and the “bad” LDL kind (s2). They are almost the same thing: correlation 0.90. Put both in a model and the model cannot tell which one deserves the credit.

FORMULA = "target ~ bmi + bp + s1 + s2 + s5"
resamples = [
    diabetes.sample(len(diabetes), replace=True, random_state=rng) for _ in range(1_000)
]
slopes = np.array([smf.ols(FORMULA, r).fit().params[["s1", "s2"]] for r in resamples])


def plot_measures(ax: Axes) -> None:
    """Total against LDL cholesterol."""
    ax.scatter(diabetes.s1, diabetes.s2, s=8, color=style.BLUE, alpha=0.5)
    ax.set(title="Two measures", xlabel="total cholesterol (s1)", ylabel="LDL cholesterol (s2)")


def plot_slopes(ax: Axes) -> None:
    """The pair of slopes in every resample."""
    ax.scatter(slopes[:, 0], slopes[:, 1], s=6, color=style.CORAL, alpha=0.4)
    ax.axhline(0, color=style.INK, lw=0.8)
    ax.axvline(0, color=style.INK, lw=0.8)
    ax.set(title="Their slopes", xlabel="slope of s1", ylabel="slope of s2")


fig, axes = style.subplots(2, height=3.3, width=9.5)
plot_measures(axes[0])
plot_slopes(axes[1])
fig.tight_layout()
plt.show()

Figure 4.2.7: Left: the two cholesterol measures of 442 diabetes patients move together. Right: their slopes in 1,000 bootstrap resamples of a model for disease progression that has both. The estimates swing widely and trade off against each other: when one is high, the other is low.

This is collinearity. The model as a whole predicts fine, but the single slopes are unstable, and here they even have opposite signs: more total cholesterol “lowers” progression while more LDL “raises” it, which is nonsense as a statement about either one. The standard errors are large, honestly reporting that the data cannot separate the two.

What to do: keep only one of them, combine them, or accept that only their joint effect is known. The variance inflation factor (VIF) measures how much a predictor’s variance (its squared standard error) is inflated by the others. Regress predictor j on all the other predictors; if that regression has $R^2_j$, then

$$\text{VIF}_j = \frac{1}{1 - R^2_j}$$

A predictor unrelated to the others has VIF 1. One that the others predict with $R^2_j = 0.9$ has VIF 10: its standard error is $\sqrt{10} \approx 3$ times what it would be on its own.

design = sm.add_constant(diabetes[PREDICTORS])
vif = {
    name: variance_inflation_factor(design.to_numpy(), i)
    for i, name in enumerate(design.columns)
    if name != "const"
}
print(
    "| predictor | VIF |\n|---|---:|\n"
    + "\n".join(f"| {name} | {value:.1f} |" for name, value in vif.items())
)

Table 4.2.2: Variance inflation factors in the model of Table 1.4.2.1. The cholesterol measurements s1 to s4 and s5 are largely predictable from each other; values above 5 or 10 are the usual warning.

predictor	VIF
age	1.2
bmi	1.5
bp	1.4
s1	59.1
s2	39.1
s3	15.4
s4	8.8
s5	10.1
s6	1.5

4.2.10 In one picture

• Each slope in a multiple regression is the effect of its predictor with the others held fixed; the added-variable plot shows it.
• $\hat{\boldsymbol\beta} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}$; each coefficient has its own t-test.
• The F-test compares nested models: are several coefficients zero together?
• Adjusted R² charges for predictors; plain R² never goes down.
• Leaving out a confounder biases a slope, and can even flip it: Simpson’s paradox.
• Groups enter as indicator columns against a reference group.
• An interaction lets one predictor’s slope depend on another.
• Predictors that carry the same information make single slopes unstable (collinearity); the VIF $= 1/(1 - R^2_j)$ measures it.
• Adjusting for variables in the model is not the same as showing a cause.

Next: how to check that a regression’s assumptions hold.