2.4 Several variables at once

Most questions involve more than one variable, and then what matters is how they depend on each other. Conditioning, slicing the data by one variable to look at another, is the single most useful idea in statistics.

import matplotlib.pyplot as plt
import numpy as np
import statsmodels.api as sm
from matplotlib.axes import Axes
from scipy import stats
from sklearn import datasets
from sklearn.mixture import GaussianMixture

from statly import style

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

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
SPECIES_COLORS = np.array([style.BLUE, style.GREEN, style.CORAL])

2.4.1 Joint, marginal and conditional

Two variables measured on the same units have a joint distribution: how common each combination is. The scatter plot shows it. Each variable on its own has a marginal distribution, the histogram along its edge. And for each value of one variable, the other has a conditional distribution: the distribution in a thin slice.

income, food = engel.income.to_numpy(), engel.foodexp.to_numpy()
SLICES = [(300, 600), (600, 900), (900, 1300), (1300, 2600)]


def plot_scatter_with_margins(fig: plt.Figure, spec) -> None:
    """A scatter plot with histograms of each variable along its edges."""
    grid = spec.subgridspec(
        2, 2, width_ratios=[4, 1], height_ratios=[1, 4], hspace=0.05, wspace=0.05
    )
    main = fig.add_subplot(grid[1, 0])
    top = fig.add_subplot(grid[0, 0], sharex=main)
    side = fig.add_subplot(grid[1, 1], sharey=main)
    main.scatter(income, food, s=8, color=style.BLUE, alpha=0.5)
    for low, high in SLICES:
        main.axvspan(
            low,
            high,
            color=style.GREY,
            alpha=0.06 if SLICES.index((low, high)) % 2 else 0.14,
            lw=0,
        )
    top.hist(income, bins=40, color=style.BLUE, alpha=0.55)
    side.hist(food, bins=40, color=style.BLUE, alpha=0.55, orientation="horizontal")
    for ax in (top, side):
        ax.set_axis_off()
    main.set(
        xlabel="household income", ylabel="food spending", xlim=(250, 2700), ylim=(150, 1700)
    )
    top.set_title("Joint and marginal", loc="left")


def plot_slices(ax: Axes) -> None:
    """Food spending within each income slice, as boxplots."""
    groups = [food[(income >= low) & (income < high)] for low, high in SLICES]
    ax.boxplot(
        groups,
        widths=0.5,
        patch_artist=True,
        boxprops=dict(facecolor=style.BLUE, alpha=0.3, color=style.INK),
        medianprops=dict(color=style.CORAL, lw=2),
        whiskerprops=dict(color=style.INK),
        capprops=dict(color=style.INK),
        flierprops=dict(markeredgecolor=style.INK),
    )
    ax.set_xticks(range(1, 5), [f"{low} to {high}" for low, high in SLICES], fontsize=8.5)
    ax.set(title="Conditional on income", xlabel="income slice", ylabel="food spending")


fig = plt.figure(figsize=(10, 3.9))
outer = fig.add_gridspec(1, 2, width_ratios=[1, 1.1], wspace=0.3)
plot_scatter_with_margins(fig, outer[0])
plot_slices(fig.add_subplot(outer[1]))
plt.show()

Figure 2.4.1: Left: income and food spending of the 235 households (the joint distribution), with each variable’s own histogram along the edges (the marginal distributions). Right: food spending within four income slices (conditional distributions). The slices move up with income, and they widen.

The conditional mean, the average food spending at each income, is exactly what the regression function $f(x)$ of Chapter 1.3 describes: $f(x) = E(Y \mid X = x)$, read “the expected value of Y given that X equals x”. The vertical bar means “given”. Regression is the study of conditional distributions.

2.4.2 Independence

Two variables are independent when knowing one tells you nothing about the other: every conditional distribution is the same as the marginal. Then the joint probability is just the product of the two:

$$P(X = x \text{ and } Y = y) = P(X = x)\, P(Y = y)$$

Correlation measures only straight-line dependence, so independence is stronger than zero correlation:

def plot_cloud(ax: Axes, x: np.ndarray, y: np.ndarray, title: str, color: str) -> None:
    """A scatter plot with its correlation."""
    ax.scatter(x, y, s=8, color=color, alpha=0.5)
    r = round(np.corrcoef(x, y)[0, 1], 2) + 0.0
    style.note(ax, f"r = {r:.2f}", x=0.03, y=0.98, color=style.INK)
    ax.set(title=title)


fig, axes = style.subplots(3, height=3, width=10.5)
plot_cloud(axes[0], income, rng.permutation(food), "Shuffled: independent", style.GREY)
plot_cloud(axes[1], income, food, "Real: dependent", style.BLUE)
x = rng.uniform(-1, 1, 300)
plot_cloud(axes[2], x, x**2 + rng.normal(0, 0.03, 300), "Dependent, r ≈ 0", style.GOLD)
for ax in axes[:2]:
    ax.set(xlabel="household income", ylabel="food spending")
axes[2].set(xlabel="x", ylabel="y")
fig.tight_layout()
plt.show()

Figure 2.4.2: Left: the incomes paired with the food spending of a randomly chosen other household: independent by construction, a shapeless cloud. Middle: the real pairs: dependent. Right: dependent but uncorrelated: y is fully determined by x (a parabola plus a little noise), yet the correlation is about 0.

Independent implies uncorrelated, but not the other way round. The shuffled panel is also the picture behind the permutation tests of Chapter 3.1: shuffling one column makes the two independent, and shows what “no relation” looks like.

2.4.3 Conditional probability and Bayes’ rule

The probability of A given B is the share of B cases in which A also happens:

$$P(A \mid B) = \frac{P(A \text{ and } B)}{P(B)}$$

The order matters, and mixing it up is one of the most common errors in reasoning about chance. A screening test catches 90% of people with a disease and wrongly alarms 5% of healthy people. One person in a hundred has the disease. Someone tests positive. How likely is it that they have it?

PEOPLE, PREVALENCE, SENSITIVITY, FALSE_ALARMS = 1000, 0.01, 0.9, 0.05
sick = np.zeros(PEOPLE, bool)
sick[rng.choice(PEOPLE, int(PEOPLE * PREVALENCE), replace=False)] = True
flagged = np.zeros(PEOPLE, bool)
flagged[np.where(sick)[0][: int(sick.sum() * SENSITIVITY)]] = True
flagged[
    rng.choice(np.where(~sick)[0], int(round((~sick).sum() * FALSE_ALARMS)), replace=False)
] = True


def plot_people(ax: Axes) -> None:
    """A grid of 1,000 people: sick in red, flagged ones circled."""
    rows, cols = 25, 40
    x, y = np.meshgrid(np.arange(cols), np.arange(rows))
    x, y = x.ravel(), y.ravel()
    ax.scatter(
        x, y, s=14, c=np.where(sick, style.CORAL, style.GREY), alpha=np.where(sick, 1.0, 0.35)
    )
    ax.scatter(
        x[flagged], y[flagged], s=55, facecolors="none", edgecolors=style.INK, linewidths=1
    )
    ax.set(xlim=(-1, cols), ylim=(-1, rows), aspect="equal")
    ax.set_axis_off()


fig, ax = plt.subplots(figsize=(8, 5.2))
plot_people(ax)
plt.show()

Figure 2.4.3: 1,000 people, 10 of whom have the disease (red). The test flags 9 of the 10 sick people and 50 of the 990 healthy ones (circled). Of the 59 flagged people, only 9 are sick: a positive test means a 15% chance of disease, not 90%.

Counting is the safest way to get it right. Written as a formula, it is Bayes’ rule, which turns $P(B \mid A)$ around into $P(A \mid B)$:

$$P(\text{sick} \mid +) = \frac{P(+ \mid \text{sick})\,P(\text{sick})}{P(+ \mid \text{sick})\,P(\text{sick}) + P(+ \mid \text{healthy})\,P(\text{healthy})} = \frac{0.9 \times 0.01}{0.9 \times 0.01 + 0.05 \times 0.99} = 0.15$$

The base rate, how common the disease is to begin with, dominates. The same logic explains why a p-value is not the probability that the null hypothesis is true (Chapter 3.1), and why precision can be low for a rare class even with a good classifier (Chapter 6.1).

2.4.4 The multivariate normal

The normal distribution has a version for several variables at once. For two, its density looks like a hill, and its contours are ellipses. The means set the centre; the covariance matrix sets the shape:

$$\Sigma = \begin{pmatrix} \sigma_x^2 & \text{cov}(x, y) \\ \text{cov}(x, y) & \sigma_y^2 \end{pmatrix}$$

The variances on the diagonal stretch the hill along each axis; the covariance tilts it.

def plot_bivariate(ax: Axes, r: float) -> None:
    """A sample, density contours and the axes of the ellipse for correlation r."""
    cov = np.array([[1, r], [r, 1]])
    sample = rng.multivariate_normal([0, 0], cov, 500)
    ax.scatter(sample[:, 0], sample[:, 1], s=5, color=style.BLUE, alpha=0.35)
    gx, gy = np.meshgrid(np.linspace(-3.5, 3.5, 200), np.linspace(-3.5, 3.5, 200))
    density = stats.multivariate_normal([0, 0], cov).pdf(np.dstack([gx, gy]))
    ax.contour(gx, gy, density, levels=4, colors=style.INK, linewidths=0.8)
    values, vectors = np.linalg.eigh(cov)
    for value, vector in zip(values, vectors.T):
        ax.annotate(
            "",
            xy=vector * 2 * np.sqrt(value),
            xytext=(0, 0),
            arrowprops=dict(arrowstyle="->", color=style.CORAL, lw=2),
        )
    ax.set(
        title=f"correlation {r}",
        aspect="equal",
        xlim=(-3.5, 3.5),
        ylim=(-3.5, 3.5),
        xlabel="x",
        ylabel="y",
    )


fig, axes = style.subplots(3, height=3.3, width=10.5)
for ax, r in zip(axes, (0, 0.6, -0.9)):
    plot_bivariate(ax, r)
fig.tight_layout()
plt.show()

Figure 2.4.4: Samples from two-variable normal distributions with standard deviations 1 and 1, and correlations 0, 0.6 and −0.9, with contours of the density. The arrows are the ellipse’s axes: the directions of largest and smallest spread, which principal component analysis finds (Chapter 7.1).

In a multivariate normal, every slice is normal, and the conditional mean of y given x is a straight line: linear regression is exactly right. With more variables the ellipse becomes an ellipsoid in more dimensions, still described by a mean vector and a covariance matrix.

2.4.5 Latent variables

Sometimes the most important variable is one we never measured. A latent variable is hidden, but leaves its mark on the variables we see. The two-humped petal lengths of Chapter 1.1 come from a hidden variable, the species. Suppose we did not know it. A mixture model assumes the data come from a few hidden groups, each with its own normal distribution, and finds the groups and each flower’s chance of belonging to them:

petals = iris[["petal length (cm)", "petal width (cm)"]].to_numpy()
mixture = GaussianMixture(3, n_init=5, random_state=0).fit(petals)
found = mixture.predict(petals)
order = np.argsort(mixture.means_[:, 0])  # name the groups from short to long petals
relabel = np.empty(3, int)
relabel[order] = np.arange(3)
found = relabel[found]


def draw_ellipses(ax: Axes) -> None:
    """Two-standard-deviation ellipses of each fitted group."""
    angles = np.linspace(0, 2 * np.pi, 100)
    for g in range(3):
        values, vectors = np.linalg.eigh(mixture.covariances_[g])
        circle = np.c_[np.cos(angles), np.sin(angles)] * 2 * np.sqrt(values)
        ellipse = circle @ vectors.T + mixture.means_[g]
        ax.plot(ellipse[:, 0], ellipse[:, 1], color=SPECIES_COLORS[relabel[g]], lw=1.5)


fig, axes = style.subplots(3, height=3.2, width=10.5, sharey=True)
axes[0].scatter(petals[:, 0], petals[:, 1], s=10, color=style.GREY, alpha=0.7)
axes[0].set_title("What we see")
axes[1].scatter(petals[:, 0], petals[:, 1], s=10, c=SPECIES_COLORS[found], alpha=0.7)
draw_ellipses(axes[1])
axes[1].set_title("Groups found")
axes[2].scatter(petals[:, 0], petals[:, 1], s=10, c=SPECIES_COLORS[iris.target], alpha=0.7)
axes[2].set_title("True species")
for ax in axes:
    ax.set_xlabel("petal length (cm)")
axes[0].set_ylabel("petal width (cm)")
fig.tight_layout()
plt.show()

Figure 2.4.5: Petal length and width of the 150 flowers. Left: what we see, with no labels. Middle: a mixture of three two-variable normals fitted without labels; the ellipses are the groups it found, and each flower is coloured by its most likely group. Right: the true species. The hidden groups match the species closely.

A latent variable can also be a hidden amount rather than a group. Several test scores are correlated because they all depend on one underlying ability; several measurements of a flower are correlated because they all depend on its overall size. Factor analysis looks for a few such hidden amounts that explain the correlations among many measured ones:

size = rng.normal(0, 1, 400)
measured = np.c_[0.9 * size, 0.8 * size, 0.7 * size] + rng.normal(0, 0.5, (400, 3))
NAMES = ["sepal length", "petal length", "petal width"]


def plot_structure(ax: Axes) -> None:
    """The hidden variable with arrows to the three measured ones."""
    ax.add_patch(plt.Circle((0.5, 0.78), 0.13, fill=False, ls="--", ec=style.CORAL, lw=1.5))
    ax.text(
        0.5, 0.78, "size\n(hidden)", ha="center", va="center", fontsize=9, color=style.CORAL
    )
    for x, name in zip((0.15, 0.5, 0.85), NAMES):
        ax.add_patch(plt.Rectangle((x - 0.13, 0.12), 0.26, 0.16, fill=False, ec=style.INK))
        ax.text(x, 0.2, name.replace(" ", "\n"), ha="center", va="center", fontsize=8.5)
        ax.annotate(
            "",
            xy=(x, 0.29),
            xytext=(0.5, 0.65),
            arrowprops=dict(arrowstyle="->", color=style.INK, lw=1),
        )
    ax.set(xlim=(0, 1), ylim=(0, 1), aspect="equal")
    ax.set_axis_off()


def plot_correlations(ax: Axes, values: np.ndarray, title: str) -> None:
    """The correlation matrix of the three measurements."""
    matrix = np.corrcoef(values.T)
    ax.imshow(matrix, cmap="RdBu_r", vmin=-1, vmax=1)
    for i in range(3):
        for j in range(3):
            ax.text(
                j,
                i,
                f"{matrix[i, j]:.2f}",
                ha="center",
                va="center",
                fontsize=10,
                color="white" if abs(matrix[i, j]) > 0.6 else "black",
            )
    short = [n.split()[0][0] + ". " + n.split()[1] for n in NAMES]
    ax.set_xticks(range(3), short, fontsize=8)
    ax.set_yticks(range(3), short, fontsize=8)
    ax.grid(False)
    ax.set_title(title)


residuals = measured - np.outer(
    size, np.linalg.lstsq(size[:, None], measured, rcond=None)[0][0]
)
fig, axes = style.subplots(3, height=3.2, width=10.5)
plot_structure(axes[0])
plot_correlations(axes[1], measured, "Measured")
plot_correlations(axes[2], residuals, "Size held fixed")
fig.tight_layout()
plt.show()

Figure 2.4.6: Simulated: one hidden quantity (a flower’s overall size) drives three measurements, each with its own noise. Left: the structure, with the hidden variable in the dashed circle. Middle: the three measurements are correlated with each other. Right: once the hidden size is known, the correlations vanish: the measurements only ever shared the size.

That is the signature of a latent variable: correlations that disappear once it is held fixed. Principal components (Chapter 7.1) are a close relative, and many methods in machine learning, from mixture models to the hidden layers of a neural network (Chapter 8.1), are ways of finding useful hidden variables.

2.4.6 In one picture

• The joint distribution shows combinations; marginals show each variable alone; conditional distributions show one variable within a slice of another.
• Regression is the conditional mean: $f(x) = E(Y \mid X = x)$.
• Independent means every slice looks the same; independence implies zero correlation, not the reverse.
• Bayes’ rule turns $P(B \mid A)$ into $P(A \mid B)$; the base rate matters.
• A multivariate normal is set by a mean and a covariance matrix; its contours are ellipses.
• Latent variables are hidden causes; mixture models and factor analysis look for them.

Next: extreme values, the laws of the largest observations.