3.6 Correlation

A correlation measures how closely two variables follow a straight line. It is easy to compute and easy to misread: it sees only straight lines, it is fragile with small samples and single outliers, and it says nothing about cause.

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

from statly import style

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

crime = sm.datasets.statecrime.load_pandas().data  # crime and poverty in 50 US states and DC
states = crime.drop("District of Columbia")


def correlated(r: float, n: int) -> tuple[np.ndarray, np.ndarray]:
    """Two normal variables with true correlation r."""
    x = rng.normal(0, 1, n)
    return x, r * x + np.sqrt(1 - r**2) * rng.normal(0, 1, n)


def plot_pair(
    ax: Axes, x: np.ndarray, y: np.ndarray, title: str, color: str = style.BLUE
) -> None:
    """A scatter plot with its correlation in the corner."""
    ax.scatter(x, y, s=10, color=color, alpha=0.6)
    r = round(np.corrcoef(x, y)[0, 1], 2) + 0.0  # + 0.0 turns −0.00 into 0.00
    style.note(ax, f"r = {r:.2f}", x=0.03, y=0.98, color=style.INK)
    ax.set_title(title)

3.6.1 What r measures

Chapter 1.2 defined the Pearson correlation as the covariance divided by both standard deviations. Equivalently, turn both variables into z-scores and average their products:

$$r = \frac{1}{n - 1} \sum_{i=1}^{n} z_{x,i}\, z_{y,i}$$

A point up and to the right of the centre, or down and to the left, has a positive product; the other two corners give negative products. If most points sit in the positive corners, r is positive. It runs from −1 (a perfect falling line) through 0 (no straight-line relation) to +1 (a perfect rising line).

fig, axes = style.subplots(6, cols=3, height=2.6, width=10.5)
for ax, r in zip(axes[:3], (0, 0.5, 0.95)):
    plot_pair(ax, *correlated(r, 100), f"Correlation {r}")
plot_pair(axes[3], crime.poverty.to_numpy(), crime.hs_grad.to_numpy(), "US states", style.GREEN)
axes[3].set(xlabel="poverty (%)", ylabel="high school graduates (%)")
x = np.linspace(-1, 1, 100)
plot_pair(axes[4], x, x**2 + rng.normal(0, 0.03, 100), "A curve", style.GOLD)
x, y = rng.normal(0, 1, 100), rng.normal(0, 1, 100)
x[0], y[0] = 9, 9  # one extreme point
plot_pair(axes[5], x, y, "One outlier", style.CORAL)
fig.tight_layout()
plt.show()

Figure 3.6.1: Top: simulated data with correlations from 0 to 0.95. Bottom left: poverty and high school graduation in the US states. Bottom middle: a perfect curve with a correlation near 0, because correlation sees only straight lines. Bottom right: one point can create a correlation on its own.

A correlation of 0.5 looks weaker than people expect. One way to read r: $r^2$ is the share of the variation in one variable that a straight line on the other explains (the R² of Chapter 4.1). An r of 0.5 explains 25%; an r of 0.3, only 9%.

Two more ways to see the same number. First, standardise both variables to z-scores and fit a least squares line: its slope is r. If y followed x perfectly, one standard deviation up in x would mean one standard deviation up in y, a slope of 1; with weaker correlation, y rises less. Second, compare the squared residuals around the line with the squared distances around a flat line at the mean of y. The line removes a share of the squares, and r is the square root of that share (with the sign of the slope).

zx = stats.zscore(states.poverty.to_numpy())
zy = stats.zscore(states.hs_grad.to_numpy())
slope = np.polyfit(zx, zy, 1)[0]
removed = 1 - ((zy - slope * zx) ** 2).sum() / (zy**2).sum()


def plot_standardised(ax: Axes) -> None:
    """Both variables as z-scores, with the least squares line."""
    ax.scatter(zx, zy, s=14, color=style.BLUE, alpha=0.7)
    grid = np.array([-2.5, 2.5])
    ax.plot(grid, slope * grid, color=style.CORAL, lw=2)
    style.note(
        ax,
        f"slope = {slope:.2f}\nr = {np.corrcoef(zx, zy)[0, 1]:.2f}",
        x=0.55,
        y=0.95,
        color=style.INK,
    )
    ax.set(title="z-scores", xlabel="poverty (z)", ylabel="graduation (z)", aspect="equal")


def plot_squares(ax: Axes, fitted: np.ndarray, title: str, color: str) -> None:
    """Each point's vertical distance to `fitted` drawn as a square."""
    for x0, y0, f0 in zip(zx, zy, fitted):
        side = y0 - f0
        corner = (x0, min(y0, f0))
        ax.add_patch(plt.Rectangle(corner, abs(side), abs(side), color=color, alpha=0.18, lw=0))
    ax.scatter(zx, zy, s=10, color=style.INK)
    order = np.argsort(zx)
    ax.plot(zx[order], fitted[order], color=color, lw=2)
    style.note(
        ax, f"sum of squares {((zy - fitted) ** 2).sum():.0f}", x=0.03, y=0.08, color=style.INK
    )
    ax.set(title=title, xlabel="poverty (z)", aspect="equal", xlim=(-3, 3), ylim=(-3, 3))


fig, axes = style.subplots(3, height=3.3, width=10.5)
plot_standardised(axes[0])
plot_squares(axes[1], np.zeros_like(zy), "Around the mean", style.GREY)
plot_squares(axes[2], slope * zx, f"Around the line: {removed:.0%} less", style.CORAL)
fig.tight_layout()
plt.show()

Figure 3.6.2: One correlation, three readings, for poverty and high school graduation in the 50 states. Left: in z-scores the least squares line has slope r. Middle: squared distances to a flat line at the mean. Right: squared residuals around the regression line. The line removes the share r² of the squares; its square root, with the sign of the slope, is r.

3.6.2 Testing a correlation

Is a correlation further from 0 than chance allows? With no true correlation, sample correlations scatter around 0, and the scatter shrinks with the sample size. The test statistic

$$t = \frac{r\sqrt{n - 2}}{\sqrt{1 - r^2}}$$

follows a t distribution with n − 2 degrees of freedom.

def plot_null_r(ax: Axes, n: int) -> None:
    """Sample correlations of unrelated variables, with the 5% cut-offs."""
    x, y = rng.normal(size=(2, 10_000, n))
    xs, ys = x - x.mean(axis=1, keepdims=True), y - y.mean(axis=1, keepdims=True)
    r = (xs * ys).sum(axis=1) / np.sqrt((xs**2).sum(axis=1) * (ys**2).sum(axis=1))
    style.hist(ax, r, bins=np.linspace(-1, 1, 61))
    t_cut = stats.t.ppf(0.975, n - 2)
    r_cut = t_cut / np.sqrt(n - 2 + t_cut**2)
    for edge in (-r_cut, r_cut):
        ax.axvline(edge, color=style.CORAL, lw=1.4)
    style.note(ax, f"±{r_cut:.2f}", x=0.03, color=style.CORAL)
    ax.set(title=f"{n} pairs", xlabel="correlation")


fig, axes = style.subplots(3, height=2.8, width=10.5, sharex=True)
for ax, n in zip(axes, (10, 30, 100)):
    plot_null_r(ax, n)
axes[0].set_ylabel("samples")
fig.tight_layout()
plt.show()

Figure 3.6.3: Correlations in 10,000 samples of truly unrelated variables, for samples of 10, 30 and 100. The red lines cut off the outer 5%: the smallest correlation that would be significant. With 10 pairs, r must be beyond ±0.63; with 100, ±0.20 is enough.

So the same r can be noise or a clear finding, depending on n. For the 51 states, poverty and high school graduation have r = −0.74: t = −7.7, p < 0.0001. A significant correlation can still be tiny: with 10,000 pairs, r = 0.02 is significant and explains 0.04% of the variation.

3.6.3 An interval for a correlation

The sampling distribution of r is lopsided when the true correlation is far from 0, because r cannot pass 1. Fisher’s z-transformation fixes that:

$$z = \tfrac{1}{2} \ln\frac{1 + r}{1 - r} \quad (= \operatorname{artanh} r), \qquad \text{SE}(z) \approx \frac{1}{\sqrt{n - 3}}$$

On the z scale the sampling distribution is close to normal, so take z ± 1.96 standard errors and turn the ends back into correlations.

TRUE_R, N = 0.8, 20
samples = [correlated(TRUE_R, N) for _ in range(10_000)]
rs = np.array([np.corrcoef(x, y)[0, 1] for x, y in samples])


def plot_r_scale(ax: Axes) -> None:
    """The sampling distribution of r."""
    style.hist(ax, rs, bins=50)
    ax.axvline(TRUE_R, color=style.INK, ls="--", lw=1.2)
    ax.set(title="r", xlabel="sample correlation", ylabel="samples")


def plot_z_scale(ax: Axes) -> None:
    """The same correlations after Fisher's z-transformation, with the normal curve."""
    z = np.arctanh(rs)
    width = style.hist(ax, z, bins=50, color=style.GREEN)
    grid = np.linspace(z.min(), z.max(), 200)
    ax.plot(
        grid,
        stats.norm.pdf(grid, np.arctanh(TRUE_R), 1 / np.sqrt(N - 3)) * len(z) * width,
        color=style.CORAL,
    )
    ax.axvline(np.arctanh(TRUE_R), color=style.INK, ls="--", lw=1.2)
    ax.set(title="z = artanh(r)", xlabel="Fisher's z", ylabel="samples")


fig, axes = style.subplots(2, height=3, width=9.5)
plot_r_scale(axes[0])
plot_z_scale(axes[1])
fig.tight_layout()
plt.show()

Figure 3.6.4: Correlations in 10,000 samples of 20 pairs from a population with true correlation 0.8 (dashed). Left: on the r scale the distribution is squeezed against 1 and lopsided. Right: after Fisher’s z-transformation it is close to the normal curve (red) with standard error 1/√17.

For the states, r = −0.74 with n = 51 gives a 95% interval from −0.85 to −0.59. The bootstrap (Chapter 3.3) is the alternative when the data are far from normal.

3.6.4 Rank correlations

Two correlations work on order instead of values, so they ignore how extreme a value is and catch any relation that keeps rising (or falling), straight or not.

• Spearman’s ρ is Pearson’s r computed on the ranks.
• Kendall’s τ looks at every pair of points. A pair is concordant if the point that is further right is also higher, discordant if it is lower. τ is (concordant − discordant) / all pairs.

few = crime.sample(12, random_state=3)


def plot_pairs(ax: Axes) -> None:
    """Concordant and discordant pairs among a few states."""
    x, y = few.poverty.to_numpy(), few.murder.to_numpy()
    for i in range(len(x)):
        for j in range(i + 1, len(x)):
            same = (x[i] - x[j]) * (y[i] - y[j]) > 0
            ax.plot(
                [x[i], x[j]],
                [y[i], y[j]],
                color=style.BLUE if same else style.CORAL,
                lw=0.7,
                alpha=0.5,
            )
    ax.scatter(x, y, s=24, color=style.INK, zorder=3)
    tau = stats.kendalltau(x, y).statistic
    style.note(ax, f"τ = {tau:.2f}", x=0.03, y=0.98, color=style.INK)
    ax.set(title="Pairs", xlabel="poverty (%)", ylabel="murders per 100,000")


def plot_three(ax: Axes) -> None:
    """Pearson, Spearman and Kendall, with and without DC."""
    names = ["Pearson r", "Spearman ρ", "Kendall τ"]
    for shift, data, color, label in (
        (-0.15, crime, style.CORAL, "with DC"),
        (0.15, states, style.BLUE, "without DC"),
    ):
        values = [
            stats.pearsonr(data.poverty, data.murder).statistic,
            stats.spearmanr(data.poverty, data.murder).statistic,
            stats.kendalltau(data.poverty, data.murder).statistic,
        ]
        ax.barh(np.arange(3) + shift, values, height=0.3, color=color, alpha=0.7, label=label)
    ax.set_yticks(range(3), names)
    ax.invert_yaxis()
    ax.legend(fontsize=8, loc="lower right")
    ax.set(title="Poverty and murder, 51 states", xlabel="correlation", xlim=(0, 0.7))


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

Figure 3.6.5: Left: Kendall’s τ for 12 states, poverty against murder rate. Every pair of points is joined by a line: blue if the poorer state also has more murders (concordant), red if fewer (discordant). Right: the three correlations for all 51 states, with and without Washington DC. Pearson’s r moves the most when one extreme point is removed; the rank correlations hardly move.

Kendall’s τ is usually smaller in size than Spearman’s ρ for the same data; they are on different scales, so compare like with like. Both have their own tests, and both are the right choice for ordered categories (such as “how religious, 1 to 4”).

3.6.5 Partial correlation

Two variables can correlate because a third drives both. The partial correlation of x and y given z removes z from both: regress x on z and keep the residuals, regress y on z and keep the residuals, and correlate the two sets of residuals.

def residuals(target: np.ndarray, given: np.ndarray) -> np.ndarray:
    """What is left of target after a straight-line fit on given."""
    fit = sm.OLS(target, sm.add_constant(given)).fit()
    return fit.resid


poverty, violent, single = (
    states.poverty.to_numpy(),
    states.violent.to_numpy(),
    states.single.to_numpy(),
)
fig, axes = style.subplots(2, height=3.2, width=9.5)
plot_pair(axes[0], poverty, violent, "Plain correlation")
axes[0].set(xlabel="poverty (%)", ylabel="violent crimes per 100,000")
plot_pair(
    axes[1],
    residuals(poverty, single),
    residuals(violent, single),
    "Single parenthood removed",
    style.PURPLE,
)
axes[1].set(xlabel="poverty, residual", ylabel="violent crime, residual")
fig.tight_layout()
plt.show()

Figure 3.6.6: Left: violent crime and poverty in the 50 states: r = 0.31. Right: the same after removing the share of single-parent households from both (each point is a state’s residual): the partial correlation is about zero. Poverty and violent crime were linked through single parenthood (compare Chapter 4.2).

This is exactly what multiple regression does for each predictor (Chapter 4.2). It removes only the variables you name: a partial correlation cannot adjust for a confounder nobody measured.

3.6.6 Traps

Correlation is the most misused number in statistics. Besides the curves and outliers of the first figure:

score, grade = correlated(0.6, 1_000)
admitted = score > np.quantile(score, 0.8)
temperature = rng.uniform(5, 30, 200)
ice_cream = 20 * temperature + rng.normal(0, 80, 200)
drownings = rng.poisson(0.4 * temperature)  # whole numbers, never negative


def plot_range(ax: Axes) -> None:
    """All applicants, with the admitted ones highlighted and both correlations."""
    ax.scatter(score[~admitted], grade[~admitted], s=6, color=style.GREY, alpha=0.4)
    ax.scatter(score[admitted], grade[admitted], s=8, color=style.BLUE, alpha=0.7)
    all_r = np.corrcoef(score, grade)[0, 1]
    top_r = np.corrcoef(score[admitted], grade[admitted])[0, 1]
    style.note(
        ax, f"all: r = {all_r:.2f}\nadmitted: r = {top_r:.2f}", x=0.03, y=0.98, color=style.INK
    )
    ax.set(title="Restriction of range", xlabel="entrance test score", ylabel="later grade")


fig, axes = style.subplots(2, height=3.2, width=9.5)
plot_range(axes[0])
plot_pair(axes[1], ice_cream, drownings, "A third variable", style.GOLD)
axes[1].set(xlabel="ice cream sold", ylabel="drownings")
fig.tight_layout()
plt.show()

Figure 3.6.7: Left: restriction of range. Simulated entrance test scores and later grades for 1,000 applicants (r = 0.6). Among the 20% admitted (blue), the correlation is much weaker, and the test looks useless. Right: a third variable. Ice cream sales and drownings both rise with temperature (simulated); they correlate without one causing the other.

• Restriction of range: correlations shrink when only part of the range is observed. Studies of selected groups (admitted students, hired workers, treated patients) understate relationships.
• Third variables: correlation is not causation. Hot days cause both.
• Averages: correlations between group averages (countries, states) are usually much stronger than between the individuals inside them, the ecological fallacy.

3.6.7 In one picture

Correlation	Measures	Robust to outliers?	Test
Pearson’s r	straight-line relation	no	$t = r\sqrt{n-2}/\sqrt{1-r^2}$, t with n − 2 df
Spearman’s ρ	any rising or falling relation, via ranks	yes	its own, or the same t
Kendall’s τ	concordant minus discordant pairs	yes	its own
partial correlation	relation after removing named variables	no	t with n − 2 − (number removed) df

• r is the average product of z-scores; r² is the share of variation explained.
• Whether an r is significant depends on n; whether it matters depends on its size.
• Use Fisher’s z for intervals; ranks for outliers and curves that keep rising.
• Always look at the scatter plot, and remember the third variable.

Next: agreement and reliability, when two measurements should give the same answer.