Chapter 91 · Solutions
Correlation vs. Causation — Worked Solutions ✅
Five challenges, each verified in code.
⚙️ Setup¶
In [1]:
import numpy as np, pandas as pd
from scipy import stats
import statsmodels.api as sm
BASE="https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"
rng=np.random.default_rng(910)
CHALLENGE 1
Manufacture a confounded correlation
Make z drive both x and y; show x and y correlate even though neither causes the other.
In [2]:
z=rng.normal(0,1,400); x=2*z+rng.normal(0,0.6,400); y=1.5*z+rng.normal(0,0.6,400)
print(f"corr(x,y) = {stats.pearsonr(x,y)[0]:.2f} (spurious)")
corr(x,y) = 0.88 (spurious)
CHALLENGE 2
Control for the confounder
Compute the partial correlation of x and y given z and show it collapses.
In [3]:
def resid(v,w): return sm.OLS(v,sm.add_constant(w)).fit().resid
print(f"partial(x,y|z) = {stats.pearsonr(resid(x,z),resid(y,z))[0]:.2f}")
partial(x,y|z) = -0.07
CHALLENGE 3
Spurious by chance
Across 1000 unrelated 30-point pairs, find the largest |r| that appears by luck.
In [4]:
best=max(abs(np.corrcoef(rng.normal(0,1,30),rng.normal(0,1,30))[0,1]) for _ in range(1000))
print(f"largest |r| from pure noise = {best:.2f}")
largest |r| from pure noise = 0.55
CHALLENGE 4
Partial-correlation formula
Verify the residual method matches the formula r_xy.z = (r_xy - r_xz r_yz)/sqrt((1-r_xz^2)(1-r_yz^2)).
In [5]:
rxy=stats.pearsonr(x,y)[0]; rxz=stats.pearsonr(x,z)[0]; ryz=stats.pearsonr(y,z)[0]
formula=(rxy-rxz*ryz)/np.sqrt((1-rxz**2)*(1-ryz**2))
print(f"residual method = {stats.pearsonr(resid(x,z),resid(y,z))[0]:.3f}, formula = {formula:.3f}")
residual method = -0.067, formula = -0.067
CHALLENGE 5
Real data: debunk the ice-cream link
Load correlation-vs-causation--confounding.xlsx; show the raw correlation and the partial correlation controlling for temperature.
In [6]:
try: d = pd.read_excel("../../data/correlation-vs-causation--confounding.xlsx", sheet_name="Days")
except FileNotFoundError: d = pd.read_excel(BASE+"correlation-vs-causation--confounding.xlsx", sheet_name="Days")
def res(v): return sm.OLS(d[v],sm.add_constant(d.temperature_f)).fit().resid
print(f"raw = {d.ice_cream_sales.corr(d.drownings):.2f}, partial|temp = {stats.pearsonr(res("ice_cream_sales"),res("drownings"))[0]:.2f}")
raw = 0.58, partial|temp = -0.03
Statistics, Data Science and AI: A Visual Handbook · © 2026 John Fisher