""" Tutorial 85 - Simpson's paradox: when every split disagrees with the total. Reads the bundled nba_league_shots.csv - 25,000 real shot-level records from one NBA season - so it runs offline. Inside it sits a genuine aggregation reversal, found by scanning the data, not constructed: Lauri Markkanen shoots a BETTER percentage than Shai Gilgeous-Alexander on two-pointers AND a better percentage on three-pointers, yet his overall field-goal percentage is LOWER. Nothing is wrong with the arithmetic. An overall percentage is a weighted average of the split percentages, weighted by each player's shot mix - and when the mix differs enough, the total can contradict every split inside it. We surface the flip, open up the weighted-average machinery that causes it, then scan every pair of high-volume shooters to show this is a standing hazard of aggregated data (53 reversal pairs), not a freak. No randomness anywhere - every number is deterministic and re-runs identically. Run: python downloads/85_simpsons_paradox_in_sports_data.py """ import itertools import os import matplotlib.pyplot as plt import pandas as pd import sdt_common as sdt sdt.init("simpsons-paradox-in-sports-data") HERE = os.path.dirname(os.path.abspath(__file__)) shots = pd.read_csv(os.path.join(HERE, "nba_league_shots.csv")) A, B = "Lauri Markkanen", "Shai Gilgeous-Alexander" TWO, THREE = "2PT Field Goal", "3PT Field Goal" # --- the aggregate answer --------------------------------------------------------- pair = shots[shots["PLAYER_NAME"].isin([A, B])] with sdt.snippet("overall"): print(f"{len(shots):,} shots, {shots['PLAYER_NAME'].nunique()} players in the bundled file\n") overall = (pair.groupby("PLAYER_NAME")["SHOT_MADE"] .agg(made="sum", attempts="count", fg_pct="mean") .sort_values("fg_pct", ascending=False)) print(overall.round(3).to_string()) print(f"\noverall field goal %: {B} ahead by " f"{overall.loc[B, 'fg_pct'] - overall.loc[A, 'fg_pct']:+.3f}") # --- the same shots, split by type ------------------------------------------------ with sdt.snippet("the-flip"): split = (pair.groupby(["SHOT_TYPE", "PLAYER_NAME"])["SHOT_MADE"] .agg(made="sum", attempts="count", fg_pct="mean")) print(split.round(3).to_string()) for t in (TWO, THREE): gap = split.loc[(t, A), "fg_pct"] - split.loc[(t, B), "fg_pct"] print(f"\n{t}: {A} ahead by {gap:+.3f}") print(f"\n{A} shoots better on twos AND on threes - yet trails overall.") print("that is Simpson's paradox, and every number above is real.") # --- why: an overall % is a mix-weighted average of the split %s ------------------- with sdt.snippet("mechanism"): print("share of each player's attempts by shot type:\n") mix = (pair.groupby(["PLAYER_NAME", "SHOT_TYPE"])["SHOT_MADE"].count() .groupby(level=0).transform(lambda s: s / s.sum()) .unstack().round(3)) print(mix.to_string()) print("\noverall = (share of twos) x (2PT%) + (share of threes) x (3PT%):\n") for p in (A, B): s2 = mix.loc[p, TWO] p2 = split.loc[(TWO, p), "fg_pct"] s3 = mix.loc[p, THREE] p3 = split.loc[(THREE, p), "fg_pct"] print(f" {p:<24} {s2:.3f} x {p2:.3f} + {s3:.3f} x {p3:.3f} = " f"{s2 * p2 + s3 * p3:.3f}") league = shots.groupby("SHOT_TYPE")["SHOT_MADE"].mean() print(f"\ntwos are simply the easier shot class: league-wide in this file,") print(f"{league[TWO]:.1%} of two-pointers go in vs {league[THREE]:.1%} of threes.") print(f"\nthe hidden confounder is the mix: {B.split()[0]} takes 83% of his") print(f"shots as (easier) twos; {A.split()[0]} takes 54% as threes. the") print("aggregate rewards the shot diet, not the shooting.") # --- how common is this? scan every high-volume pair ------------------------------- MIN_TOTAL, MIN_SPLIT = 100, 15 by_type = (shots.groupby(["PLAYER_NAME", "SHOT_TYPE"])["SHOT_MADE"] .agg(["mean", "count"]).unstack()) totals = shots.groupby("PLAYER_NAME")["SHOT_MADE"].agg(["mean", "count"]) eligible = [p for p in totals.index if totals.loc[p, "count"] >= MIN_TOTAL and by_type.loc[p, ("count", TWO)] >= MIN_SPLIT and by_type.loc[p, ("count", THREE)] >= MIN_SPLIT] reversals = [] for a, b in itertools.combinations(eligible, 2): for x, y in ((a, b), (b, a)): # test both directions of the flip wins_both = (by_type.loc[x, ("mean", TWO)] > by_type.loc[y, ("mean", TWO)] and by_type.loc[x, ("mean", THREE)] > by_type.loc[y, ("mean", THREE)]) if wins_both and totals.loc[x, "mean"] < totals.loc[y, "mean"]: margin = min(by_type.loc[x, ("mean", TWO)] - by_type.loc[y, ("mean", TWO)], by_type.loc[x, ("mean", THREE)] - by_type.loc[y, ("mean", THREE)]) reversals.append((margin, x, y)) reversals.sort(reverse=True) with sdt.snippet("scan"): n_pairs = len(eligible) * (len(eligible) - 1) // 2 print(f"players with >= {MIN_TOTAL} shots and >= {MIN_SPLIT} attempts of each " f"type: {len(eligible)}") print(f"pairs scanned: {n_pairs:,}") print(f"Simpson reversals found: {len(reversals)} ({len(reversals) / n_pairs:.1%} of pairs)\n") print("strongest five (by the smaller of the two split gaps):") for margin, x, y in reversals[:5]: print(f" {x} beats {y} on both splits (weaker gap " f"{margin:+.3f}) yet trails overall") print(f"\nthe featured pair is #1 of {len(reversals)} - the strongest reversal") print("in the file, found by exhaustive scan, not construction.") # --- the chart: the flip, then the lever that causes it ---------------------------- fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11.2, 4.9)) ORANGE = sdt.sport_color("basketball") SLATE = "#2C5E8A" cats = ["2PT%", "3PT%", "overall FG%"] vals_a = [split.loc[(TWO, A), "fg_pct"], split.loc[(THREE, A), "fg_pct"], totals.loc[A, "mean"]] vals_b = [split.loc[(TWO, B), "fg_pct"], split.loc[(THREE, B), "fg_pct"], totals.loc[B, "mean"]] xs = range(3) W = 0.36 ax1.bar([x - W / 2 for x in xs], vals_a, W, color=ORANGE, label="Markkanen") ax1.bar([x + W / 2 for x in xs], vals_b, W, color=SLATE, label="Gilgeous-Alexander") for x, (va, vb) in enumerate(zip(vals_a, vals_b)): ax1.annotate(f"{va:.3f}", xy=(x - W / 2, va), xytext=(0, 3), textcoords="offset points", ha="center", fontsize=9) ax1.annotate(f"{vb:.3f}", xy=(x + W / 2, vb), xytext=(0, 3), textcoords="offset points", ha="center", fontsize=9) winner = "Markkanen" if va > vb else "SGA" ax1.annotate(f"{winner} +{abs(va - vb):.3f}", xy=(x, max(va, vb) + 0.055), ha="center", fontsize=9, fontweight="bold", color=ORANGE if va > vb else SLATE) ax1.set_xticks(list(xs)) ax1.set_xticklabels(cats) ax1.set_ylim(0, 0.76) ax1.set_ylabel("make rate") ax1.set_title("Both splits say Markkanen; the total says SGA") ax1.legend(loc="upper right", frameon=False, fontsize=9) # right panel: overall %ages live on a line from 3PT% (mix=0% twos) to 2PT% (100% twos) for p, color, name in ((A, ORANGE, "Markkanen"), (B, SLATE, "Gilgeous-Alexander")): p2 = split.loc[(TWO, p), "fg_pct"] p3 = split.loc[(THREE, p), "fg_pct"] share2 = (pair[pair.PLAYER_NAME == p]["SHOT_TYPE"] == TWO).mean() ov = totals.loc[p, "mean"] ax2.plot([0, 1], [p3, p2], color=color, lw=2.0, alpha=0.85) ax2.scatter([share2], [ov], s=70, color=color, zorder=3) ax2.annotate(f"{name}\noverall {ov:.3f} at {share2:.0%} twos", xy=(share2, ov), xytext=(share2 - 0.03, ov + (0.05 if p == A else -0.10)), fontsize=9, color="#20242B", ha="right" if p == A else "left") ax2.annotate(f"3PT {p3:.3f}", xy=(0, p3), xytext=(-8, 0), fontsize=8, textcoords="offset points", ha="right", va="center", color=color) ax2.annotate(f"2PT {p2:.3f}", xy=(1, p2), xytext=(8, 0), fontsize=8, textcoords="offset points", ha="left", va="center", color=color) ax2.set_xlim(-0.22, 1.22) ax2.set_ylim(0.26, 0.68) ax2.set_xticks([0, 0.25, 0.5, 0.75, 1]) ax2.set_xticklabels(["0%", "25%", "50%", "75%", "100%"]) ax2.set_xlabel("share of attempts that are two-pointers") ax2.set_ylabel("make rate") ax2.set_title("The total is a lever: mix decides where it lands") fig.suptitle("Simpson's paradox, live in real NBA shot data", fontweight="bold") fig.tight_layout(rect=(0, 0.015, 1, 1)) sdt.save_fig(fig, "reversal", source="Bundled sample of 25,000 real NBA shot records") print("done")