C10 · Who draws fouls in the NBA? — Bayesian Modelling in Practice

The brief

Close NBA games are decided at the free-throw line. A team's analytics group is preparing for the trade deadline and asks two questions:

Which players are best at drawing fouls late in close games, and which defenders are best at not being whistled?
Both answers must be separated from the matchup: a guard who is always defended by foul-prone opponents looks good for free, and a centre who only ever guards stars looks bad for free.

Their intern ranked players by raw foul rate. The top of the list is a set of players nobody has heard of, each at 100%. They would like something better - and they want to know how far to trust any ranking at all.

Since 2015 the league has reviewed every material play in the last two minutes of close games and published its verdict. Each row of the data is one reviewed play.

How this notebook works

Each task states what to deliver, not how. Write your code in the YOUR CODE HERE cells.
Stuck? h.hint("task2") reveals hints one level at a time: nudge → approach → code skeleton. Try to get by on nudges.
h.check("task2", sigma_theta=...) compares your numbers with the reference solution.
A full worked solution lives in notebooks/solutions/. Open it only when you are done (or truly stuck).

The models here have about 1,600 parameters and 47,000 observations. pm.sample(nuts_sampler="nutpie") fits each in roughly a minute; budget accordingly.

import time

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
from scipy import stats
from scipy.special import expit

from pymc_challenges import Hints, data

RANDOM_SEED = 2021
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-variat")

h = Hints("C10")
h.tasks()

C10 - Who draws the whistle? An item-response model of NBA foul calls

task0 - Define the outcome, look at raw rates (3 hints, 2 check(s))
task1 - A colleague's Rasch model and what is wrong with it (3 hints, 1 check(s))
task2 - Identification and centred vs non-centred (3 hints, 3 check(s))
task3 - Position as a group-level predictor (3 hints, 2 check(s))
task4 - Posterior predictive check by position pair (3 hints)
task5 - Shrinkage - raw ranking vs posterior ranking (3 hints, 1 check(s))
task6 - Rank uncertainty and a specific matchup (3 hints, 2 check(s))

The data

column	meaning
`committing`	the player who committed the foul - or would have, had one been called
`disadvantaged`	the player on the receiving end
`decision`	the league's review: `CC` correct call, `IC` incorrect call, `CNC` correct no-call, `INC` incorrect no-call
`committing_position`, `disadvantaged_position`	G(uard), F(orward), C(entre) and hybrids such as `G-F`

data.describe("nba_fouls")
fouls = data.load("nba_fouls")
fouls.head()

nba_fouls: Foul calls / non-calls in close games, with committing and disadvantaged player.
  source: NBA Last Two Minute reports 2015-2021, via pymc-examples.
  url:    https://raw.githubusercontent.com/pymc-devs/pymc-examples/main/examples/data/item_response_nba.csv

	committing	disadvantaged	decision	committing_position	disadvantaged_position
play_id
1	Avery Bradley	Stephen Curry	CNC	G	G
3	Isaiah Thomas	Andre Iguodala	CNC	G	G-F
4	Jae Crowder	Harrison Barnes	INC	F	F
6	Draymond Green	Isaiah Thomas	CC	F	G
7	Avery Bradley	Stephen Curry	CC	G	G

fouls.decision.value_counts()

decision
CNC    33401
CC     10498
INC     2540
IC       422
Name: count, dtype: int64

Task 0 · Define the outcome, then look at raw rates

Deliver

A binary outcome built from decision, with a justification: the four codes allow more than one defensible definition - which one answers the team's question?
The intern's table: the ten players with the highest raw rate of drawing fouls, with the number of plays each rate is based on.
A plot of raw foul-drawing rate against number of plays for every player, and one sentence on why ranking by raw rate is a bad idea.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

decision crosses two things: was a whistle blown and was that correct. Which of the two produces free throws? Then, before trusting any rate, ask how many plays it is based on - what is the raw rate of a player with one play?

2 Approach

"A foul was called" is decision in {CC, IC}; "a foul occurred" would be {CC, INC}. fouls.groupby("disadvantaged").foul_called.agg(plays="size", raw_rate="mean") gives the table. Plot raw_rate against plays on a log x axis; adding the binomial 95% range around the league-wide rate (scipy.stats.binom.interval) turns it into a funnel plot that shows how much of the spread is pure noise.

3 Skeleton

fouls["foul_called"] = fouls.decision.isin([___, ___]).astype(int)
drawn = fouls.groupby("disadvantaged").foul_called.agg(plays="size", raw_rate="mean")
drawn.sort_values(["raw_rate", "plays"], ascending=False).head(10)
(drawn.plays <= 5).sum()

Self-checks in the notebook (h.check("task0", …)):

call_rate — overall share of reviewed plays on which a foul was CALLED (correctly or not), about 0.23
players_with_5_or_fewer — number of players with at most 5 plays as the disadvantaged player, a bit over 200

Task 1 · A colleague's Rasch model

Item-response theory was built for exams: the probability that student $i$ answers question $j$ correctly depends on the student's ability minus the question's difficulty. Here the "student" is the disadvantaged player and the "question" is the committing player:

$$\text{logit}\,P(\text{foul called}) = \mu + \theta_{\text{disadvantaged}} - b_{\text{committing}}$$

A colleague proposes the natural hierarchical version, "with weakly-informative priors":

$$\theta_i \sim \text{Normal}(\mu_\theta, \sigma_\theta), \quad b_j \sim \text{Normal}(\mu_b, \sigma_b), \quad \mu \sim \text{Normal}(0, 1.5), \quad \mu_\theta, \mu_b \sim \text{Normal}(0, 5), \quad \sigma_\theta, \sigma_b \sim \text{HalfNormal}(1)$$

Deliver

This model, fitted to all plays, with named dimensions for the two sets of players.
Diagnostics for the five top-level parameters. Something is badly wrong. Show with a plot which direction in parameter space the sampler cannot resolve.
Find a function of $\mu, \mu_\theta, \mu_b$ that is well determined - report its posterior mean and r_hat - and explain in two sentences what the data can and cannot tell you in this model.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Look at r_hat and ESS for mu, mu_theta and mu_b, not just at the player effects. Then think about the likelihood: it only ever sees mu + theta_i - b_j. What happens to that sum if you add 1 to every theta and subtract 1 from mu?

2 Approach

Use pd.factorize(column, sort=True) to get an integer index array and the labels for coords; index the parameter vectors with those arrays (theta[d_idx]). Sample with nuts_sampler="nutpie".

Scatter the draws of mu_theta against mu_b coloured by chain: each chain sits on a different piece of a long ridge. Only one linear combination of the three location parameters enters the likelihood for an average matchup. Build it from idata.posterior (xarray arithmetic keeps the chain/draw dims) and pass combo.to_dataset() to az.summary to get its own r_hat.

3 Skeleton

d_idx, d_names = pd.factorize(fouls.disadvantaged, sort=True)
c_idx, c_names = pd.factorize(fouls.committing, sort=True)
coords = {"disadvantaged": d_names, "committing": c_names}
with pm.Model(coords=coords) as colleague_model:
    ...
    theta = pm.Normal("theta", mu_theta, sigma_theta, dims="disadvantaged")
    b = pm.Normal("b", mu_b, sigma_b, dims="committing")
    pm.Bernoulli("foul_called", logit_p=mu + theta[d_idx] - b[c_idx], observed=y)
    idata = pm.sample(nuts_sampler="nutpie", random_seed=RANDOM_SEED)

post = idata.posterior
combo = (post["mu"] ___ post["mu_theta"] ___ post["mu_b"]).rename("combo")
az.summary(combo.to_dataset())

Self-checks in the notebook (h.check("task1", …)):

identified_combination — posterior mean of mu + mu_theta - mu_b, the log-odds of a call for an average player against an average opponent, about -1.14

Task 2 · Identify it, then worry about geometry

Deliver

An identified version of the model, with a justification of the constraint you chose (there is more than one), and a prior predictive check on the probability scale.
Fits of both a centred and a non-centred parameterisation of the player effects. Compare them on divergences and on the effective sample size of $\mu$, $\sigma_\theta$ and $\sigma_b$, and choose one.
An explanation of the difference in terms of the data: which players favour which parameterisation, and how many of each are there? Show the relevant posterior geometry for one player of each kind.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Three location parameters, one degree of freedom in the data: you must remove two. What is the least arbitrary way to do that, given that you care about players relative to each other?

For the geometry: a player with 600 plays has a posterior dominated by his likelihood; a player with one play has a posterior that is essentially his prior Normal(0, sigma). What does the joint posterior of (effect, log sigma) look like in the second case? Count how many players are of each kind.

2 Approach

Identification: either fix the population means at zero (theta ~ Normal(0, sigma_theta)) or use pm.ZeroSumNormal("theta", sigma, dims=...).

Non-centred: z ~ Normal(0, 1) and theta = pm.Deterministic("theta", sigma * z). Build the model in a function with a centred flag, fit both with nutpie and compare idata.sample_stats["diverging"].sum() and the ess_bulk column of az.summary(idata, var_names=["mu", "sigma_theta", "sigma_b"]).

To see the geometry, scatter log(sigma_theta) against the effect of a one-play player and against that of the busiest player, in whichever coordinates the sampler actually works in (theta for centred, z_theta for non-centred).

3 Skeleton

def rasch(centred):
    with pm.Model(coords=coords) as model:
        mu = pm.Normal("mu", -1, 1)
        sigma_theta = pm.HalfNormal("sigma_theta", ___)
        sigma_b = pm.HalfNormal("sigma_b", ___)
        if centred:
            theta = pm.Normal("theta", 0, sigma_theta, dims="disadvantaged")
            b = ___
        else:
            z_theta = pm.Normal("z_theta", 0, 1, dims="disadvantaged")
            theta = pm.Deterministic("theta", ___, dims="disadvantaged")
            ...
        pm.Bernoulli("foul_called", logit_p=mu + theta[d_idx] - b[c_idx], observed=y)
    return model

Self-checks in the notebook (h.check("task2", …)):

mu — posterior mean of the intercept (log-odds of a call, average player vs average opponent), about -1.14
sigma_theta — posterior mean sd of foul-drawing skill across players, about 0.30 on the logit scale
sigma_b — posterior mean sd of the committing-player effect across players, about 0.41

Task 3 · Does position explain it?

Guards handle the ball and attack the rim; centres set screens and protect it. Give both skills a group-level predictor: a player's effect is drawn around the mean of his position.

Deliver

The model with position means for both $\theta$ and $b$, identified, with clean diagnostics.
The posterior of the difference between guards (G) and centres (C) in foul-drawing skill, and the same for $b$. Do guards draw more fouls than centres?
What happens to $\sigma_\theta$ and $\sigma_b$ compared with Task 2, and what does that tell you? Be careful about what a play is in this dataset before you interpret the position effects on $b$.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

A group-level predictor changes what a player is shrunk towards: not the league average, but the average of players like him. You now have a new location invariance (position means vs the intercept) - deal with it the same way as before.

When you interpret b by position, remember how a play gets into this dataset: the league reviews contact. What kind of contact involves centres on almost every possession?

2 Approach

Each player has exactly one position in the file, so build a player-level position index: fouls.groupby("disadvantaged").disadvantaged_position.first().reindex(d_names) and turn it into codes with pd.Categorical(..., categories=POSITIONS).codes.

pos_theta = pm.ZeroSumNormal("pos_theta", sigma, dims="position") gives seven position means that sum to zero, so mu stays identified. Then theta = pos_theta[d_pos_idx] + sigma_theta * z_theta (non-centred), and likewise for b. The guard-centre contrast is post["pos_theta"].sel(position="G") - post["pos_theta"].sel(position="C").

3 Skeleton

with pm.Model(coords={**coords, "position": POSITIONS}) as position_model:
    mu = pm.Normal("mu", -1, 1)
    pos_theta = pm.ZeroSumNormal("pos_theta", ___, dims="position")
    pos_b = pm.ZeroSumNormal("pos_b", ___, dims="position")
    z_theta = pm.Normal("z_theta", 0, 1, dims="disadvantaged")
    theta = pm.Deterministic("theta", pos_theta[___] + sigma_theta * z_theta, dims="disadvantaged")
    ...
    pm.Bernoulli("foul_called", logit_p=mu + theta[d_idx] - b[c_idx], observed=y)

Self-checks in the notebook (h.check("task3", …)):

guard_minus_centre_theta — posterior mean of pos_theta[G] - pos_theta[C] on the logit scale - the sign may surprise you; about -0.17
sigma_b — posterior mean of the within-position sd of the committing-player effect, about 0.28 (down from 0.41)

Task 4 · Check it where it could fail

The model is additive on the logit scale: a guard's skill is the same whether a guard or a centre is defending him. That is a strong claim with 49 position pairings to test it on.

Deliver a posterior predictive check of the foul-call rate for every pair of (disadvantaged position, committing position), for the models of Task 2 and Task 3, with a verdict: is either model adequate for the team's purpose? Where, if anywhere, does additivity fail?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

A PPC is only as sharp as the statistic. The overall call rate will be reproduced by almost any model with an intercept. Pick a statistic that the model's structure could get wrong: here, the call rate within each of the 49 position pairings.

2 Approach

You do not need pm.sample_posterior_predictive for 47,000 plays times 4,000 draws. Take a few hundred draws with az.extract(idata, var_names=["mu", "theta", "b"], num_samples=200), compute p = expit(mu + theta[d_idx] - b[c_idx]) as a (plays, draws) array in NumPy, simulate rng.random(p.shape) < p, and aggregate by pair with np.bincount(pair_idx, weights=...). Compare the observed rate in each pair with the 94% interval of the replicated rates; with 49 intervals expect about three misses by chance.

3 Skeleton

pair_idx, pair_names = pd.factorize(
    fouls.disadvantaged_position + " vs " + fouls.committing_position, sort=True)
pair_n = np.bincount(pair_idx)
pair_obs = np.bincount(pair_idx, weights=y) / pair_n

draws = az.extract(idata, var_names=["mu", "theta", "b"], num_samples=200)
theta = draws["theta"].transpose("disadvantaged", "sample").values
b = draws["b"].transpose("committing", "sample").values
p = expit(draws["mu"].values[None, :] + theta[d_idx] - ___)
y_rep = rng.random(p.shape) < p
rep = np.stack([np.bincount(pair_idx, weights=y_rep[:, k]) for k in range(200)], axis=1) / pair_n[:, None]
lo, hi = np.quantile(rep, [0.03, 0.97], axis=1)

Task 5 · What shrinkage buys you

Deliver

For every player, the model's estimate of his foul-drawing rate against an average defender, plotted against his raw rate, with the number of plays visible in the plot.
The top ten by raw rate next to the top ten by posterior mean. How many plays does the least-observed player in each list have?
Two sentences for the intern on what the model did to the 100% players, and why that is not "throwing data away".

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

To compare with a raw rate you need the model's answer on the probability scale, for a standardised opponent. What value of b represents "an average defender"?

2 Approach

expit(mu + theta_i) per posterior draw, averaged over draws, is player i's call rate against a defender with b = 0. Scatter it against the raw rate with marker size or colour showing the number of plays: low-n players collapse onto their position mean, high-n players stay near the diagonal (shifted by their opponents' quality).

3 Skeleton

theta_draws = az.extract(idata, var_names="theta").transpose("sample", "disadvantaged").values
mu_draws = az.extract(idata, var_names="mu").values
model_rate = expit(___).mean(0)        # against an average defender: which terms remain?
post_top = pd.Series(theta_draws.mean(0), index=d_names).sort_values(ascending=False).head(10)
drawn.plays[___].min()

Self-checks in the notebook (h.check("task5", …)):

min_plays_in_posterior_top10 — number of plays of the least-observed player among the top ten by posterior mean theta - dozens, not one or two

Task 6 · Rankings the front office can use

A posterior mean is still a point estimate, and the front office will read a ranked list as if rank 3 were better than rank 7.

Deliver

For the leading candidates: the posterior probability of being among the ten best foul-drawers in the league, and a 94% interval for each player's rank. Is anybody in the top ten with probability above one half?
The same treatment, briefly, for defenders who avoid the whistle.
A specific matchup: James Harden has the ball and Rudy Gobert is defending. The posterior of P(foul called) with an interval - and the same with an average centre in place of Gobert, to show what the defender is worth.
Three sentences for the front office on how to use (and not use) the list.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

A rank is a function of all the thetas at once, so compute it within each posterior draw and then summarise across draws. "Probability of being in the top ten" is then just a mean of an indicator.

2 Approach

With draws of shape (sample, player), (-draws).argsort(axis=1).argsort(axis=1) + 1 gives each player's rank in each draw. (ranks <= 10).mean(0) is P(top 10); np.quantile(ranks, [0.03, 0.97], axis=0) the rank interval.

For the matchup, combine the draws of mu, theta[James Harden] and b[Rudy Gobert] and push them through expit. For "an average centre" replace Gobert's b by the centre position mean.

3 Skeleton

ranks = (-theta_draws).argsort(axis=1).argsort(axis=1) + 1
p_top10 = pd.Series((ranks <= 10).mean(0), index=d_names).sort_values(ascending=False)

i = list(d_names).index("James Harden")
j = list(c_names).index("Rudy Gobert")
p = expit(mu_draws + theta_draws[:, i] - ___)
p.mean(), np.quantile(p, [0.03, 0.97])

Self-checks in the notebook (h.check("task6", …)):

max_p_top10 — the largest posterior probability any player has of being in the top ten foul-drawers - well short of certainty, roughly 0.3-0.6
p_harden_vs_gobert — posterior mean P(foul called) when Harden is the disadvantaged player and Gobert the committing player, about 0.25

Going further

One player, two skills. Most players appear in both roles. Give each player a bivariate effect $(\theta_i, b_i)$ with an LKJCholeskyCov prior. Are good foul-drawers also whistle-prone defenders?
Matchup interactions. Task 4 tests additivity only at the position level. Add a position-pair interaction (ZeroSumNormal over two dims) and compare with LOO. Is it worth 49 more parameters?
Referee error instead of whistles. Refit with the outcome "the call was wrong" (IC or INC). Do some players get the benefit of the doubt?
Does skill drift? The file has no dates, but the reports do. Get the season for each play from the original L2M data and let $\theta$ follow a random walk across seasons.

h.progress()

C10: 0/21 hints used, 11/11 checks passed.

Reference solution

A complete, executed solution with commentary on every task — including the ways the obvious model fails. Try the hint ladder first.

Open the worked solution →

#The brief

#How this notebook works

#The data

#Task 0 · Define the outcome, then look at raw rates

#Task 1 · A colleague's Rasch model

#Task 2 · Identify it, then worry about geometry

#Task 3 · Does position explain it?

#Task 4 · Check it where it could fail

#Task 5 · What shrinkage buys you

#Task 6 · Rankings the front office can use

#Going further

Reference solution

The brief

How this notebook works

The data

Task 0 · Define the outcome, then look at raw rates

Task 1 · A colleague's Rasch model

Task 2 · Identify it, then worry about geometry

Task 3 · Does position explain it?

Task 4 · Check it where it could fail

Task 5 · What shrinkage buys you

Task 6 · Rankings the front office can use

Going further