↑ ↓ to move · enter to open

D03 · Statistics · Pearson r, confounders, experiments

correlation vs causation

When two things rise and fall together, they are correlated. That fact alone can't tell you whether one causes the other. This page measures together-ness with one number, then shows three ways it misleads: a hidden third thing driving both, an arrow that points the other way, and plain luck. Then it runs the tool that settles the question: random assignment. All data here is synthetic and generated in your browser.

LAB 1

measuring together-ness

Put two measurements on a scatter plot, one dot per thing measured. If dots that are high on one axis tend to be high on the other, the variables are positively correlated. The Pearson correlation r turns that into a number from −1 to +1: +1 means the dots sit exactly on a rising line, −1 on a falling line, 0 means no straight-line pattern at all. Click to add dots, drag to move them.

scatter · click to add
Data
Anscombe
0.00

How to read this: the gold dashed cross marks the average x and the average y. Dots up-right or down-left of it (green) push r up; dots in the other two corners (red) push it down. r runs from −1 to 1 and measures how tightly the dots follow a straight line.

  • pushes r up
  • pushes r down
  • best-fit line
  • amber sticks (Residuals): how far each dot misses the line

Try this: press Mystery plot, guess r with the slider, then Check.

—
Pearson r
—
r² · variance explained
—
Points · best-fit line
—
Guess game

The formula

Draw a cross through the average x (x̄) and average y (ȳ). For each dot multiply its two distances from the cross. Dots in the upper-right or lower-left quadrant give positive products (green), the other two quadrants give negative ones (red). Add them up and scale so the answer can't leave −1 to +1.

r = Σ(x−x̄)(y−ȳ) / √( Σ(x−x̄)² · Σ(y−ȳ)² )
= …

r² is variance explained

Switch to Residuals. The gold line is the least-squares fit, the line that makes the total squared vertical miss as small as possible. Without x, your best guess for every y is ȳ, and the squared misses add to SStot. Using the line, they add to SSres. The share removed is r².

r² = 1 − SSres / SStot
= …

Same r, different stories

In 1973 the statistician Francis Anscombe published four small data sets, now called Anscombe's quartet. All four have 11 points, the same averages, the same best-fit line y = 3.00 + 0.500x and the same r = 0.816 (set IV rounds to 0.817). Set II is a smooth curve, set III is an almost perfect line plus one outlier, and in set IV a single far point creates the whole correlation. r only measures how well a straight line fits, so plot the data before trusting the number.

CodePython

from math import sqrt

def pearson(xs, ys):
    n = len(xs)
    mx, my = sum(xs) / n, sum(ys) / n
    sxy = sum((x - mx) * (y - my) for x, y in zip(xs, ys))
    sxx = sum((x - mx) ** 2 for x in xs)
    syy = sum((y - my) ** 2 for y in ys)
    return sxy / sqrt(sxx * syy)

# Anscombe's quartet, set I
x = [10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5]
y = [8.04, 6.95, 7.58, 8.81, 8.33, 9.96,
     7.24, 4.26, 10.84, 4.82, 5.68]
r = pearson(x, y)
print(round(r, 3), round(r * r, 3))   # 0.816 0.667
LAB 2

the hidden third thing

On hot days people buy more ice cream, and more people swim, so more swimming accidents happen. Plot ice cream sales against accidents for 200 days and you get a clear positive r, even though ice cream causes no accidents. Weather here is a confounder: a third variable that drives both. The numbers are hypothetical. Click the arrows in the diagram to change what causes what, then hold temperature fixed and watch the link fade.

no direct link hot weather ice cream sales swim accidents causal model
cooler dayhotter day
  • gold arrow: this cause is on
  • dashed arrow: off
  • line: best fit (one per temperature band when bands are on)
Direct link
0.80 off
—
Raw r, all days
—
Within temperature bands
—
Partial r, temperature held fixed
—
True direct link in this model

Now

Controlling for a variable

Slide Temperature bands up. Days are split into bands of similar temperature, and r is measured inside each band, so hot days are only compared with hot days. More bands leave less temperature variation inside each one, and the within-band r approaches the partial correlation, which removes the straight-line effect of temperature from both variables at once:

rxy·z = (rxy − rxz·ryz) / √( (1 − rxz²)(1 − ryz²) )
= …

Here x is ice cream, y is accidents and z is temperature. Partial r only removes confounders you measured. A confounder nobody recorded stays hidden in the raw r.

Which way does the arrow point?

Switch the direct link between Ice cream → accidents and Accidents → ice cream. The scatter and r barely change, because r treats x and y the same: r(x, y) = r(y, x). Data on two variables can't reveal the direction of cause. This trap is reverse causation: a town with many police officers may have them because it has a lot of crime.

CodePython

import random
from math import sqrt

def pearson(xs, ys):
    n = len(xs)
    mx, my = sum(xs) / n, sum(ys) / n
    sxy = sum((x - mx) * (y - my) for x, y in zip(xs, ys))
    sxx = sum((x - mx) ** 2 for x in xs)
    syy = sum((y - my) ** 2 for y in ys)
    return sxy / sqrt(sxx * syy)

def partial_r(rxy, rxz, ryz):
    return (rxy - rxz * ryz) / sqrt((1 - rxz**2) * (1 - ryz**2))

random.seed(1)
n, w = 200, 0.8                 # days, weather effect
s = sqrt(1 - w * w)
temp = [random.gauss(0, 1) for _ in range(n)]
ice = [400 + 120 * (w * t + s * random.gauss(0, 1)) for t in temp]
acc = [12 + 3 * (w * t + s * random.gauss(0, 1)) for t in temp]

rxy = pearson(ice, acc)         # about 0.64
rxz, ryz = pearson(ice, temp), pearson(acc, temp)
print(rxy, partial_r(rxy, rxz, ryz))   # partial is near 0
Part 2 · correlation by coincidence

The third trap needs no cause at all. Generate many series of pure random numbers, unrelated by construction, and compare every pair. A few pairs will match closely by luck. A standard test calls a correlation "significant" when unrelated data would produce one that strong less than 5% of the time, so among thousands of pairs about 5% pass anyway. Searching many comparisons and reporting the winner is called the multiple comparisons problem.

all pairs · best match
40 15

How to read this: top: one bar per range of r, counting all pairs of random, unrelated series. Red bars passed the usual 5% test for "significant" purely by luck. Bottom: the luckiest pair, both rescaled to the same height.

  • "significant" by chance
  • the 5% cutoff (±r_crit)
  • "flipped": drawn upside down so a negative match lines up
—
Pairs compared
—
Strongest r found
—
Pairs "significant" (p < 0.05)
—
Expected if unrelated noise

Now

The null distribution

For n unrelated, normally distributed points, r scatters around 0 with a known shape (gold curve). The cutoff rcrit leaves 5% of that shape in the two tails:

f(r) ∝ (1 − r²)(n−4)/2 P(|r| > rcrit) = 0.05
= …

Switch to Random trends: each series is a random walk, so neighbouring points are linked. The test assumes independent points, and with trending data far more than 5% of pairs pass. Many striking correlations between yearly statistics have this shape.

CodePython

import random
from math import sqrt
from itertools import combinations

def pearson(xs, ys):
    n = len(xs)
    mx, my = sum(xs) / n, sum(ys) / n
    sxy = sum((x - mx) * (y - my) for x, y in zip(xs, ys))
    sxx = sum((x - mx) ** 2 for x in xs)
    syy = sum((y - my) ** 2 for y in ys)
    return sxy / sqrt(sxx * syy)

M, n = 40, 15
series = [[random.gauss(0, 1) for _ in range(n)]
          for _ in range(M)]
rs = [pearson(a, b) for a, b in combinations(series, 2)]
r_crit = 0.514     # |r| for p < 0.05 when n = 15
print(len(rs), max(map(abs, rs)))
print(sum(abs(r) > r_crit for r in rs) / len(rs))   # about 0.05
LAB 3

how we find out what causes what

Suppose a study app raises test scores by a known amount, the true effect. In an observational study students choose for themselves, and motivated students are both more likely to use the app and likely to score higher anyway. Motivation is a confounder nobody measured. In a randomized experiment a coin flip decides who gets the app, which cuts the arrow from motivation to app use. Both studies below use the same kind of students. Each dot is one student, colored by motivation. All numbers are hypothetical.

one study of each kind
less motivatedmore motivated
+3.0 0.80
Students per study

How to read this: each dot is a student; each column is students without or with the app. The pale bar is a column's average score. The colored vertical bar between the averages is the study's estimate of the app's effect, and the gold dashed bar beside it shows the true effect for comparison.

300 repeats
observational estimatesrandomized estimatestrue effect
—
True effect (points)
—
Observational estimate
—
Randomized estimate · 95% CI
—
Average over repeats

Now

Why the coin flip works

estimate = ȳapp − ȳno app
observational = true effect + 8c·(ūapp − ūno app) + noise
randomized: E[ūapp − ūno app] = 0

The model: score = 70 + 8c·u + effect·(app) + noise, where u is motivation and c is the confounding slider. Students choose the app with probability 1/(1+e−2cu), so app users have higher average u and the observational difference mixes their motivation into the estimate. A coin flip is unrelated to u, so both groups get the same mix of motivation on average. Repeat the study and the randomized estimates center on the true effect. Chance still moves any single estimate, and the 95% confidence interval ±1.96·SE shows how far. It narrows as the sample grows, like the margin of error in D00 · Polls. A poll randomizes who gets asked; an experiment randomizes who gets treated.

A confounder can do more than inflate an effect. It can reverse its sign, which is D02 · Simpson's paradox.

When you can't flip the coin

Many questions can't be randomized: nobody can assign smoking or poverty. Researchers then look for natural experiments, where something close to random already split people into groups, such as a lottery for school places or a rule that applies only above an age or score cutoff. Comparing people just either side of a cutoff is called regression discontinuity. Other designs compare changes over time between a group hit by a new policy and one that wasn't (difference-in-differences). Each one depends on an assumption that the split is unrelated to everything else that matters, and that assumption has to be argued from knowledge of the setting.

Without an experiment, evidence builds from several directions at once: a plausible mechanism, a dose-response pattern, the right order in time, and results that repeat across different populations and designs.

CodePython

import random
from math import exp

def study(N, effect, c, randomized):
    app, no_app = [], []
    for _ in range(N):
        u = random.gauss(0, 1)              # motivation, unmeasured
        p = 0.5 if randomized else 1 / (1 + exp(-2 * c * u))
        t = random.random() < p
        y = 70 + 8 * c * u + effect * t + random.gauss(0, 5)
        (app if t else no_app).append(y)
    return sum(app) / len(app) - sum(no_app) / len(no_app)

obs = [study(1000, 3, 0.8, False) for _ in range(300)]
rnd = [study(1000, 3, 0.8, True) for _ in range(300)]
print(sum(obs) / 300, sum(rnd) / 300)   # far above 3, about 3