Correlation vs causation: Pearson and Spearman

#Why it matters

Correlation vs causation is the distinction every stakeholder forgets and every interviewer checks. Analysts compute correlations daily (does delivery time relate to review score, does discount depth relate to repeat purchase), and the number is easy to get. The hard part is knowing what it measures, when one bad row has faked it, and what you are allowed to say about it. Statistics appears in 58% of Data Analyst and 82% of Data Scientist job descriptions we read; "correlation does not imply causation, so what would you do instead?" shows up in almost every product analytics interview.

#The concept

#Three correlation coefficients

Spearman is Pearson computed on ranks. Kendall's τ has a direct reading ("of all pairs of trips, τ more agree than disagree on ordering"), behaves better with small samples and many ties, and is usually smaller in size than Spearman for the same data. Use Pearson when the relationship looks linear and there are no wild values, Spearman for monotonic or skewed data, Kendall for small samples or ordinal scales.

#Range and strength language

All three run from -1 (perfect negative) through 0 (no association of that kind) to +1 (perfect positive). Rough language that most teams accept for |r|:

These cut-offs are conventions, not laws. In noisy behavioural data an r of 0.2 can be commercially important; in physical measurements an r of 0.9 can be disappointing. And $r^2$ is the share of variance in one variable that a straight line on the other explains: r = 0.5 explains only 25%.

#Zero correlation does not mean no relationship

Pearson only sees straight lines. A perfect U-shape (y = x², with x symmetric around zero) has r = 0. Always plot.

#Correlation is not causation: three reasons

If X and Y move together, any of these can be true:

1. Confounding: a third variable Z drives both. Ice cream sales and drownings both rise in hot weather. Temperature is the confounder.

2. Reverse causation: Y causes X. Areas with more police often have more crime, because police are sent where crime is.

3. Coincidence: with enough variables, some pairs line up by chance. Two series that both trend upward over time correlate strongly whether or not they are connected.

Only after ruling these out, usually with an experiment or a careful design (later lessons on A/B tests and causal inference), may you say "X causes Y".

#Simpson's paradox

A trend that appears in every group can reverse when the groups are combined. It happens when the grouping variable (a confounder) is related to both variables. You will reproduce one below with penguin bills, where the combined correlation is negative and every species shows a positive one.

#Restricted range

If you only observe a narrow slice of X, the correlation shrinks even when the true relationship is strong. Correlating exam scores with job performance among hired candidates (who all scored high) understates how well the exam predicts. You will see the same thing with short taxi trips.

#Correlation matrices and multiple comparisons

A matrix of 20 variables has 190 pairs. At a 5% significance level, about 9 or 10 of those will look "significant" even if every variable is random noise. Treat a large correlation matrix as a map for hypotheses, not a list of findings.

#Worked example: penguins, taxis and Anscombe

You need pandas, SciPy and DuckDB (pip install pandas scipy duckdb). The Palmer Penguins data (CC0, by Allison Horst, Alison Hill and Kristen Gorman) holds body measurements for 344 penguins of three species in Antarctica.

#Step 1: flipper length and body mass, overall and by species

import numpy as np
import pandas as pd
from scipy import stats

url = "https://raw.githubusercontent.com/allisonhorst/palmerpenguins/main/inst/extdata/penguins.csv"
cols = ["bill_length_mm", "bill_depth_mm", "flipper_length_mm", "body_mass_g"]
penguins = pd.read_csv(url).dropna(subset=cols)
print("n =", len(penguins))

r, p = stats.pearsonr(penguins["flipper_length_mm"], penguins["body_mass_g"])
print(f"overall: r = {r:.3f}, p = {p:.2g}")
for species, grp in penguins.groupby("species"):
    r = grp["flipper_length_mm"].corr(grp["body_mass_g"])
    rho = grp["flipper_length_mm"].corr(grp["body_mass_g"], method="spearman")
    print(f"{species:<10} n = {len(grp):>3}  pearson {r:.3f}  spearman {rho:.3f}")

n = 342
overall: r = 0.871, p = 4.4e-107
Adelie     n = 151  pearson 0.468  spearman 0.475
Chinstrap  n =  68  pearson 0.642  spearman 0.670
Gentoo     n = 123  pearson 0.703  spearman 0.717

Overall r = 0.87 is "very strong", but inside each species it is 0.47 to 0.70. Part of the overall correlation comes from the fact that Gentoo penguins are bigger in every way. That is not wrong, but it answers a different question ("across all penguins") from "within a species, do longer-flippered birds weigh more?".

#Step 2: Simpson's paradox with bill length and depth

r_all = penguins["bill_length_mm"].corr(penguins["bill_depth_mm"])
print(f"all penguins: r = {r_all:.3f}")
for species, grp in penguins.groupby("species"):
    print(f"{species:<10} r = {grp['bill_length_mm'].corr(grp['bill_depth_mm']):.3f}")

means = penguins.groupby("species")[["bill_length_mm", "bill_depth_mm"]].mean().round(1)
print(means)

all penguins: r = -0.235
Adelie     r = 0.391
Chinstrap  r = 0.654
Gentoo     r = 0.643
           bill_length_mm  bill_depth_mm
species                                 
Adelie               38.8           18.3
Chinstrap            48.8           18.4
Gentoo               47.5           15.0

Across all 342 penguins, longer bills go with shallower bills (r = -0.24, and significant). Within every species, longer bills go with deeper bills (r between 0.39 and 0.65). The reversal happens because Gentoo penguins have long, shallow bills and Adelie penguins short, deep ones: species is related to both variables. Ignore it and you get the wrong sign.

In business data the same pattern shows up all the time: a discount looks like it lowers revenue per order overall because it is used mostly on cheap categories, while within every category it raises it. Whenever you correlate across groups that differ, check within groups too.

#Step 3: Anscombe's quartet

Francis Anscombe built four small datasets in 1973 to make one point. The seaborn-data repository has them:

anscombe = pd.read_csv("https://raw.githubusercontent.com/mwaskom/seaborn-data/master/anscombe.csv")
summary = anscombe.groupby("dataset").apply(lambda d: pd.Series({
    "mean_x": d["x"].mean(), "mean_y": d["y"].mean(), "sd_y": d["y"].std(),
    "pearson": d["x"].corr(d["y"]), "spearman": d["x"].corr(d["y"], method="spearman"),
}), include_groups=False)
print(summary.round(3))

         mean_x  mean_y   sd_y  pearson  spearman
dataset                                          
I           9.0   7.501  2.032    0.816     0.818
II          9.0   7.501  2.032    0.816     0.691
III         9.0   7.500  2.030    0.816     0.991
IV          9.0   7.501  2.031    0.817     0.500

Same means, same SD, same Pearson r of 0.816. Plotted, they are a noisy line (I), a smooth curve (II), a perfect line with one outlier (III) and a vertical stack with one far point that creates the whole correlation (IV). Spearman already hints that they differ (0.69, 0.99, 0.50), but only a scatter plot tells you how. Draw one before you quote any correlation:

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 4, figsize=(14, 3.2), sharey=True)
for ax, (name, d) in zip(axes, anscombe.groupby("dataset")):
    ax.scatter(d["x"], d["y"])
    ax.set_title(f"{name}: r = {d['x'].corr(d['y']):.2f}")
fig.savefig("anscombe.png", dpi=120, bbox_inches="tight")

#Step 4: one bad row against 7.6 million good ones

Now the January 2019 NYC TLC yellow taxi file from Mean vs median: central tendency and outliers (download link there). Trip distance and fare should be almost perfectly related, since the meter charges by distance.

import duckdb

taxi = duckdb.sql("""
    select trip_distance, fare_amount from 'yellow_tripdata_2019-01.csv.gz'
    where tpep_pickup_datetime >= '2019-01-01' and tpep_pickup_datetime < '2019-02-01'
""").df()

def both(d):
    pearson = d["trip_distance"].corr(d["fare_amount"])
    spearman = d["trip_distance"].corr(d["fare_amount"], method="spearman")
    return f"pearson {pearson:.3f}, spearman {spearman:.3f}"

print(f"raw     n = {len(taxi):,}  {both(taxi)}")
clean = taxi[taxi["fare_amount"].between(0.01, 500) & taxi["trip_distance"].between(0.01, 100)]
print(f"clean   n = {len(clean):,}  {both(clean)}")

raw     n = 7,667,255  pearson 0.038, spearman 0.921
clean   n = 7,605,890  pearson 0.945, spearman 0.935

Pearson on the raw data says distance and fare are unrelated (0.04). After removing about 61,000 impossible records (0.8% of trips, including the 623,259.86 fare), it is 0.95. Spearman was 0.92 on the raw data all along, because a fare of 623,259.86 is just "rank 1" to it. When Pearson and Spearman disagree sharply, look for outliers before you look for meaning.

#Step 5: restricted range and small samples

for limit in [1, 2, 100]:
    short = clean[clean["trip_distance"] < limit]
    print(f"trips under {limit:>3} miles: n = {len(short):>9,}  pearson = {short['trip_distance'].corr(short['fare_amount']):.3f}")

small = clean.sample(20, random_state=1)
print("n = 20 sample:",
      "pearson", round(stats.pearsonr(small["trip_distance"], small["fare_amount"]).statistic, 3),
      "spearman", round(stats.spearmanr(small["trip_distance"], small["fare_amount"]).statistic, 3),
      "kendall", round(stats.kendalltau(small["trip_distance"], small["fare_amount"]).statistic, 3))

trips under   1 miles: n = 2,071,854  pearson = 0.107
trips under   2 miles: n = 4,688,935  pearson = 0.489
trips under 100 miles: n = 7,605,890  pearson = 0.945
n = 20 sample: pearson 0.949 spearman 0.945 kendall 0.824

The same meter, the same pricing rule, and the correlation falls from 0.95 to 0.11 when you only look at trips under a mile. On short trips, time in traffic and the fixed starting charge dominate the fare, and distance varies too little to show its effect. If your dataset only covers a narrow slice (one city, one price band, only customers who converted), expect weaker correlations than the true relationship.

#Step 6: how many "significant" correlations does noise give you?

rng = np.random.default_rng(42)
noise = pd.DataFrame(rng.normal(size=(342, 20)), columns=[f"v{i}" for i in range(20)])
pvals = []
for i in range(20):
    for j in range(i + 1, 20):
        pvals.append(stats.pearsonr(noise.iloc[:, i], noise.iloc[:, j]).pvalue)
pvals = np.array(pvals)
print(f"pairs tested: {len(pvals)}, p < 0.05: {(pvals < 0.05).sum()}, expected by chance: {0.05 * len(pvals):.1f}")

pairs tested: 190, p < 0.05: 10, expected by chance: 9.5

Twenty columns of pure random numbers, the same size as the penguin data, and a handful of pairs come out "significant". A big heatmap of correlations will always contain some. The multiple testing lesson in the inference module covers how to correct for this; for now, treat surprising cells in a correlation matrix as questions.

#Common mistakes

• Quoting Pearson on skewed or dirty data. One bad row turned 0.95 into 0.04. Compute Spearman next to it and plot.

• Reading r as a slope. r = 0.9 does not mean "fare rises 0.9 per mile". It measures how tightly points hug a line, not how steep the line is. Use regression for the slope.

• Correlating across groups that differ. Species reversed the sign of the bill correlation. Check within the obvious groups (city, category, device, plan).

• Concluding "no relationship" from r near 0. The relationship may be curved, or your range may be restricted.

• Saying "drives", "impacts" or "leads to" when you only have a correlation. Say "is associated with" and name the possible confounders.

• Correlating two time series that both trend. Two growing series correlate. Correlate the period-to-period changes instead.

#Interview angle

#Practice task

1. Compute the full Pearson and Spearman correlation matrices for the four penguin measurements, overall and for Gentoo only. Which pairs change the most?

2. Find another Simpson's paradox in the penguin data: try flipper_length_mm against bill_depth_mm.

3. In the taxi data, correlate tip_amount with fare_amount for card payments (payment_type = 1) and cash payments (2) separately. Explain the cash result from how tips are recorded.

4. Write one paragraph for a product manager explaining why a correlation of 0.6 between "used the app's wishlist" and "made a purchase" does not prove that the wishlist drives purchases, naming two confounders.

Done when each answer has a number, a plot or table, and a sentence about what you can and cannot conclude.

print(penguins[cols].corr().round(2))
print(penguins[cols].corr(method="spearman").round(2))
print(penguins.loc[penguins["species"] == "Gentoo", cols].corr().round(2))

print("all:", round(penguins["flipper_length_mm"].corr(penguins["bill_depth_mm"]), 3))
for s, g in penguins.groupby("species"):
    print(s, round(g["flipper_length_mm"].corr(g["bill_depth_mm"]), 3))

#Further reading

• SciPy pearsonr, spearmanr and kendalltau: signatures and p-value methods.

• Anscombe, Graphs in Statistical Analysis (The American Statistician, 1973): the original quartet paper.

• Palmer Penguins: the dataset's documentation and licence.

• Stanford Encyclopedia of Philosophy: Simpson's Paradox: a careful explanation of why it happens and what it means for causal claims.