Correlation, causation and Simpson's paradox

Before you start

• Level: 2 (intermediate).
• Time: about 2 hours.
• Builds on Mean, median and skew (the mean and the standard deviation), Outliers, box plots and robust summaries (why one extreme value can move a summary) and Hypothesis testing by permutation (why a randomised test lets you talk about cause). You need Python 3.9 or later.

The problem

Kettlebrook Fitness, a fictional workout app, pulled 28 days of data for 600 active users. Leah, the head of engagement, made a scatter plot: users who got more push notifications logged more workouts. She wants to raise the daily notification cap from one to three in the next release, "because notifications drive workouts", and asks you to confirm the relationship before Friday's planning meeting. You will measure the relationship with a correlation you compute by hand, then split users by how long they have been on the app and see what happens to it.

The file data/users.csv is synthetic, generated by make_data.py with a fixed seed. Each row is one made-up user. 10 users turned workout tracking off, so their workout count is blank.

In plain words

Towns with more ice cream sales have more drownings. Ice cream does not drown anyone: hot weather makes people buy ice cream and go swimming. The two numbers move together because a third thing pushes both. That third thing is called a confounder.

Kettlebrook has a likely confounder too. Users who have been on the app for a year are the committed ones. They work out more, and they have also switched on more reminder types over time, so they get more notifications. Put everyone on one scatter plot and notifications and workouts rise together, even if notifications do nothing, or even if they hurt.

The fix is to compare like with like: look at new users only, then regular users only, then veterans only. If the relationship flips direction inside every group, you have Simpson's paradox: a trend in the combined data that reverses in each part.

Words you will use:

• Correlation (Pearson's r): a number from -1 to 1 for how closely two variables follow a straight line. 1 is a perfect rising line, -1 a perfect falling line, 0 no straight-line relationship.
• Deviation: a value minus its mean. Positive means above average.
• Confounder: a variable that affects both things you are comparing, and so creates or hides a relationship between them.
• Stratify: split the data by a variable (here tenure, how long a user has been on the app) and analyse each group (each stratum) separately.
• Within-group correlation: the correlation that is left after removing the differences between groups' averages.
• Causation: changing X would change Y. Correlation alone cannot show this; a randomised experiment can.

What you will build from scratch

• read_users(path): users with a known workout count, skipping the ones who turned tracking off.
• pearson_r(xs, ys): the correlation written out from deviations, with the cases where it is undefined.
• r_by_group(users): r inside each tenure group.
• within_group_r(users): one correlation with the between-group differences removed.
• is_simpson(overall_r, group_rs): does every group point the other way?

The idea

Pearson's r, from deviations. For each user take the deviation of x (notifications) from its mean and the deviation of y (workouts) from its mean, and multiply them. If a user is above average on both, or below average on both, the product is positive. If they are above on one and below on the other, it is negative. Add the products up, then divide by the size of the deviations so the answer does not depend on units:

r = sum((x - mean_x) * (y - mean_y)) / sqrt( sum((x - mean_x)^2) * sum((y - mean_y)^2) )

r is undefined when one variable never varies (a denominator of 0) and with fewer than 2 pairs. r only measures straight lines: y = x squared on x from -2 to 2 is a perfect relationship with r = 0.

Why a confounder fools r. The causal story built into this data (which Leah cannot see) looks like this:

flowchart LR
    T["Tenure<br/>(new, regular, veteran)"] -->|"veterans opt in to more reminders"| N["Notifications in 28 days"]
    T -->|"committed users train more"| W["Workouts in 28 days"]
    N -->|"each extra notification:<br/>a small DROP (fatigue)"| W

Pooled across tenure, the two strong arrows from tenure dominate, so notifications and workouts rise together. Inside one tenure group, tenure is held fixed and only the arrow from notifications to workouts is left, and it points down.

Within-group correlation. Subtract each group's own mean from its users' x and y values, then compute r on the leftovers. Every group is now centred on (0, 0), so the gaps between groups no longer count. If this r has the opposite sign to the pooled r, the pooled relationship was mostly about who is in which group.

Correlation is not causation, even within groups. Stratifying removes the confounders you thought of. It cannot remove the ones you did not measure. To learn what sending more notifications does, randomly assign users to one or three notifications a day and compare, as in lessons 03 and 08.

Walk through

Run each block from this folder. r by hand on 5 made-up users, so you can see every step:

import math
xs = [2, 4, 6, 8, 10]    # notifications
ys = [3, 7, 5, 11, 9]    # workouts
mx, my = sum(xs) // 5, sum(ys) // 5   # both sums divide by 5 exactly
dx = [x - mx for x in xs]
dy = [y - my for y in ys]
products = [a * b for a, b in zip(dx, dy)]
print("means", mx, my)
print("dx", dx)
print("dy", dy)
print("products", products, "sum", sum(products))
r = sum(products) / math.sqrt(sum(a * a for a in dx) * sum(b * b for b in dy))
print("r", round(r, 4))

means 6 7
dx [-4, -2, 0, 2, 4]
dy [-4, 0, -2, 4, 2]
products [16, 0, 0, 8, 8] sum 32
r 0.8

No product is negative: users above average on notifications tend to be above average on workouts. The products add up to 32, the two sums of squared deviations are 40 and 40 (the square root of 40 x 40 is 40), so r = 32 / 40 = 0.8.

Now the real data, all tenures together:

import csv, math
users = []
with open("data/users.csv", newline="") as f:
    for row in csv.DictReader(f):
        if row["workouts_28d"] == "":
            continue
        users.append((row["tenure"], int(row["notifications_28d"]), int(row["workouts_28d"])))

def r(pairs):
    xs, ys = [p[0] for p in pairs], [p[1] for p in pairs]
    mx, my = sum(xs) / len(xs), sum(ys) / len(ys)
    sxy = sum((x - mx) * (y - my) for x, y in zip(xs, ys))
    return sxy / math.sqrt(sum((x - mx) ** 2 for x in xs) * sum((y - my) ** 2 for y in ys))

print(len(users), "users, pooled r =", round(r([(u[1], u[2]) for u in users]), 4))

590 users, pooled r = 0.6661

r = 0.67 is a strong positive correlation. This is Leah's scatter plot in one number. Now look at the groups:

import csv, math
users = []
with open("data/users.csv", newline="") as f:
    for row in csv.DictReader(f):
        if row["workouts_28d"] != "":
            users.append((row["tenure"], int(row["notifications_28d"]), int(row["workouts_28d"])))

def r(pairs):
    xs, ys = [p[0] for p in pairs], [p[1] for p in pairs]
    mx, my = sum(xs) / len(xs), sum(ys) / len(ys)
    sxy = sum((x - mx) * (y - my) for x, y in zip(xs, ys))
    return sxy / math.sqrt(sum((x - mx) ** 2 for x in xs) * sum((y - my) ** 2 for y in ys))

for t in ("new", "regular", "veteran"):
    g = [(u[1], u[2]) for u in users if u[0] == t]
    print("%-8s n=%d  mean notifications %5.1f  mean workouts %5.1f  r %7.4f" % (
        t, len(g), sum(p[0] for p in g) / len(g), sum(p[1] for p in g) / len(g), r(g)))

new      n=217  mean notifications  12.0  mean workouts   5.2  r -0.2226
regular  n=197  mean notifications  23.1  mean workouts  11.1  r -0.2780
veteran  n=176  mean notifications  37.9  mean workouts  17.2  r -0.3246

Veterans get about three times the notifications of new users and do about three times the workouts. But inside every group the correlation is negative: among users with the same tenure, the ones who get more notifications work out a little less. That is Simpson's paradox.

One number for "inside the groups": centre each group on its own means, then correlate:

import csv, math
users = []
with open("data/users.csv", newline="") as f:
    for row in csv.DictReader(f):
        if row["workouts_28d"] != "":
            users.append((row["tenure"], int(row["notifications_28d"]), int(row["workouts_28d"])))
dx, dy = [], []
for t in ("new", "regular", "veteran"):
    g = [u for u in users if u[0] == t]
    mx = sum(u[1] for u in g) / len(g)
    my = sum(u[2] for u in g) / len(g)
    dx += [u[1] - mx for u in g]
    dy += [u[2] - my for u in g]
sxy = sum(a * b for a, b in zip(dx, dy))
print("within-group r", round(sxy / math.sqrt(sum(a * a for a in dx) * sum(b * b for b in dy)), 4))

within-group r -0.2715

The pooled r was +0.67. With the gaps between tenure groups removed it is -0.27.

On the lesson data, python3 solution/solution.py prints:

n 590
overall_r 0.6661
by_tenure {'new': {'n': 217, 'r': -0.2226}, 'regular': {'n': 197, 'r': -0.278}, 'veteran': {'n': 176, 'r': -0.3246}}
within_r -0.2715
simpson True

One more fact, which Leah cannot see but make_data.py can: within a tenure group, each extra notification was built to cost 0.15 workouts in 28 days. Raising the cap would most likely lower workouts, the opposite of the plan.

Common mistakes

1. Turning the pooled correlation into a forecast (the trap). The slope of the best straight line through the points is r x (SD of y) / (SD of x). Pooled, it promises a gain for every notification:

import csv, math
users = []
with open("data/users.csv", newline="") as f:
    for row in csv.DictReader(f):
        if row["workouts_28d"] != "":
            users.append((row["tenure"], int(row["notifications_28d"]), int(row["workouts_28d"])))

def slope(pairs):
    xs, ys = [p[0] for p in pairs], [p[1] for p in pairs]
    mx, my = sum(xs) / len(xs), sum(ys) / len(ys)
    return sum((x - mx) * (y - my) for x, y in zip(xs, ys)) / sum((x - mx) ** 2 for x in xs)

print("pooled: %+.3f workouts per extra notification" % slope([(u[1], u[2]) for u in users]))
for t in ("new", "regular", "veteran"):
    print("%-7s %+.3f" % (t, slope([(u[1], u[2]) for u in users if u[0] == t])))

pooled: +0.314 workouts per extra notification
new     -0.109
regular -0.144
veteran -0.161

"Two more notifications a day is about 56 more in 28 days, so about 18 more workouts per user" sounds like a business case. Inside every tenure group the slope is negative. The pooled number measures who the heavy-notification users are, not what notifications do to them.

2. Trusting r without looking at the points. One unusual user can make or break a correlation. Five users with no relationship, then one test account that got 60 notifications and logged 40 workouts:

import math
def r(xs, ys):
    mx, my = sum(xs) / len(xs), sum(ys) / len(ys)
    sxy = sum((x - mx) * (y - my) for x, y in zip(xs, ys))
    return sxy / math.sqrt(sum((x - mx) ** 2 for x in xs) * sum((y - my) ** 2 for y in ys))
xs, ys = [10, 12, 14, 16, 18], [9, 6, 10, 5, 8]
print("five users:", round(r(xs, ys), 3))
print("plus one test account:", round(r(xs + [60], ys + [40]), 3))

five users: -0.229
plus one test account: 0.975

From slightly negative to nearly perfect, from one row. Plot the data, and check extreme rows before you report r.

3. Reading r = 0 as "no relationship". r measures straight lines only. Workouts that rise with notifications up to a point and then fall (a hill shape) can give r near 0 while the relationship is strong.

4. Treating the stratified answer as proof. The within-group r of -0.27 removes tenure, the confounder we measured. Users who get more notifications within a band might still differ in ways we did not record. Stratifying is how you challenge a correlation; a randomised test is how you settle what a change causes. Also, only stratify by things that are fixed before the notifications (tenure is). Splitting by something the notifications themselves change, like "opened the app this week", can create a fake relationship of its own.

Your turn

Open exercise.py and fill in the six functions. Then run, from this folder:

python3 test_lesson.py

The test checks pearson_r against a different formula for r (the one-pass sums version), on a toy example where r is exactly 0.8, on a curve where r is 0, and on inputs where r is undefined. It checks the groups on a hand-made set where the pooled r is positive and each group's r is -1, and runs your pipeline on the committed data and on a second dataset built with another seed. To see the reference solution pass: CHECK_SOLUTION=1 python3 test_lesson.py.

Stretch: use the permutation idea from lesson 08 to ask whether the new-user r of -0.22 could be luck: shuffle workouts within the new users 2,000 times and count how often |r| reaches 0.22.

Check yourself

1. What does r = 0.67 tell you, and what does it not tell you?
2. Explain Simpson's paradox to Leah using the three tenure groups.
3. Why is it fine to stratify by tenure but risky to stratify by "opened the app this week"?
4. What would you run to find out whether more notifications cause more workouts?

What to tell the stakeholder

"Across all users, more notifications go with more workouts (r = 0.67), but that is because long-time users both get more notifications and train more. Among users with the same tenure, more notifications go with slightly fewer workouts (within-group r = -0.27), so raising the cap could backfire." Caveat: this is still a correlation, so before changing the cap for everyone we recommend a randomised test of one against three notifications a day.

In real job postings

In the Real Analyst Jobs data drop of 25 October 2026, 996 of the 2,056 US analyst postings that list a tool (48.4%) ask for Python and 324 (15.8%) ask for R, the two languages where an analysis like this is usually run. See Python data analyst jobs and product analyst jobs.

Next

11: Linear regression from scratch: the slope from Common mistakes, done properly.