Rugby — hierarchical Poisson with crossed team random effects¶

From-scratch Metropolis-within-Gibbs; the PyMC Six-Nations example¶

Module: rugby_gibbs.py. The PyMC "Hierarchical model for Rugby prediction" (Baio & Blangiardo 2010), 2014 Six Nations. This is the richest of the three hierarchical-Poisson projects: every team carries two random effects — an attack and a defence — that appear, crossed, across many matches.

project random effect(s) conjugacy sampler
Pumps one Gamma rate per pump conjugate (almost) pure Gibbs
Epil one Normal intercept per patient non-conjugate Metropolis-within-Gibbs
Rugby (here) two crossed Normal effects per team (attack, defence) non-conjugate Metropolis-within-Gibbs

Model¶

Each match (home team $h$, away team $a$) produces two scores, each Poisson: $$\text{home score}\sim\text{Poisson}(\theta_h),\quad \log\theta_h=\mu+\text{home}+\text{att}_h-\text{def}_a,$$ $$\text{away score}\sim\text{Poisson}(\theta_a),\quad \log\theta_a=\mu\ \ \ \ \ \ \ +\text{att}_a-\text{def}_h.$$ A team scores more with a stronger attack $\text{att}$ and against a weaker defence $\text{def}$ (higher $\text{def}_t$ ⇒ fewer points conceded). The home side gets a fixed bonus $\text{home}$. Hierarchical priors tie the teams together: $$\text{att}_t\sim N(0,\sigma_{\text{att}}^2),\quad \text{def}_t\sim N(0,\sigma_{\text{def}}^2),\quad \mu,\text{home}\sim N(0,10^2),\quad \sigma^2_{\bullet}\sim\text{IG}(a_0,b_0).$$

Algorithm (crossed Metropolis-within-Gibbs)¶

We stack the 15 matches into 30 scoring events, each with linear predictor $\eta=\mu+\text{home}\cdot\mathbb 1_{\text{home}}+\text{att}_{\text{scorer}}-\text{def}_{\text{opp}}$. One sweep:

  1. $(\mu,\text{home})$ — block random-walk Metropolis.
  2. $\text{att}_t$ — per-team RW-Metropolis, vectorised over teams (step $\propto1/\sqrt{\text{curvature}_t}$, accept/reject all teams at once via np.bincount over the scorer index). Enters with $+$ sign.
  3. $\text{def}_t$ — same, indexed by the opponent; enters with $-$ sign.
  4. Centre att and def (subtract their means, absorb into $\mu$) — the sum-to-zero translation move, since att/def levels are confounded with $\mu$ (att adds, def subtracts in every event). Exact Gibbs step; without it the level mixes terribly.
  5. $\sigma^2_{\text{att}},\sigma^2_{\text{def}}$ — Inverse-Gamma conjugate (sum-to-zero ⇒ $T-1$ effective d.f.).

This is the same machinery as Epil's GLMM (hpois_gibbs.py), generalised to two crossed effect vectors instead of one nested intercept.

In [1]:
import numpy as np, matplotlib.pyplot as plt
from rugby_gibbs import rugby_data, rugby_gibbs, summary

y, scorer, opp, is_home, teams = rugby_data('rugby.csv', 2014); T = len(teams)
out = rugby_gibbs(y, scorer, opp, is_home, T, R=30000, burn=5000, seed=1)
print(f"acceptance  (mu,home) {out['accept_fix']:.2f}   att {out['accept_att']:.2f}   def {out['accept_def']:.2f}")
print(f"mu   = {out['mu'].mean():.3f}   (baseline score exp(mu) = {np.exp(out['mu'].mean()):.1f})")
print(f"home = {out['home'].mean():.3f}  95% CI [{np.percentile(out['home'],2.5):.2f}, {np.percentile(out['home'],97.5):.2f}]  -> home advantage x{np.exp(out['home'].mean()):.2f}")
print(f"sd_att = {out['sd_att'].mean():.3f}   sd_def = {out['sd_def'].mean():.3f}\n")
att = out['att'].mean(0); deff = out['deff'].mean(0)
print(f"{'team':<10}{'attack':>9}{'defence':>9}   (sorted by attack)")
for t in np.argsort(-att):
    print(f"{teams[t]:<10}{att[t]:>9.3f}{deff[t]:>9.3f}")

acceptance  (mu,home) 0.48   att 0.63   def 0.63
mu   = 2.661   (baseline score exp(mu) = 14.3)
home = 0.366  95% CI [0.19, 0.54]  -> home advantage x1.44
sd_att = 0.594   sd_def = 0.617

team         attack  defence   (sorted by attack)
England       0.416    0.247
Ireland       0.187    0.569
Wales         0.135    0.142
France        0.064   -0.144
Italy        -0.243   -0.523
Scotland     -0.559   -0.292
In [2]:
fig, ax = plt.subplots(1, 2, figsize=(12, 4.6))
# team-strength map: attack (x) vs defence (y); top-right = strong both ways
ax[0].axhline(0, color='grey', lw=.7); ax[0].axvline(0, color='grey', lw=.7)
ax[0].scatter(att, deff, s=70, color='steelblue', zorder=3, edgecolor='k', lw=.4)
for t in range(T):
    ax[0].annotate(teams[t], (att[t], deff[t]), fontsize=9, xytext=(5, 4), textcoords='offset points')
ax[0].set_xlabel('attack  (higher = scores more)'); ax[0].set_ylabel('defence  (higher = concedes less)')
ax[0].set_title('Team-strength map (2014 Six Nations)')
# home-advantage posterior
ax[1].hist(out['home'], bins=50, density=True, color='darkorange', alpha=.8)
ax[1].axvline(out['home'].mean(), color='k', ls='--', lw=1, label=f"mean {out['home'].mean():.2f}")
ax[1].axvline(0, color='firebrick', ls=':', lw=1, label='no home effect')
ax[1].set_xlabel('home advantage (log scale)'); ax[1].set_yticks([])
ax[1].set_title(f"Home advantage  (P>0 = {(out['home']>0).mean():.3f})"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.savefig('rugby_strength.png', dpi=120, bbox_inches='tight'); plt.show()
No description has been provided for this image

Results¶

  • Home advantage is real: $\text{home}\approx0.37$, i.e. a host scores about 1.4× as much, with essentially all posterior mass above 0.
  • Attack: England top (they put 52 on Italy, 29 on Wales), then Ireland and Wales; Scotland weakest (scored 0, 3, 6 in three games). Defence: Ireland best (highest def — they conceded the fewest points and won the 2014 title on points difference), Italy worst (conceded 46, 30, 52).
  • The team-strength map separates the table: Ireland and England sit top-right (strong attack and defence); Italy and Scotland bottom-left. This is exactly the 2014 final standings, recovered from 15 matches by pooling through $\sigma_{\text{att}},\sigma_{\text{def}}$.

Takeaways¶

  • Crossed random effects, same engine. Two effect vectors (attack/defence), each entering every match through a different team index, are handled by the same vectorised Metropolis-within-Gibbs as Epil's single nested intercept — np.bincount over the scorer index for attack, over the opponent index for defence.
  • Sum-to-zero is essential. att and def levels are confounded with $\mu$; the centring/translation move makes the level mix in one step (cf. the identical fix in hpois_gibbs.py).
  • Caveat: rugby points arrive in lumps of 3/5/7 (penalties/tries/conversions), so the Poisson is a deliberate simplification — but it captures the relative team strengths well, which is the point of the example. Cross-checks (PyMC NUTS, R glmer) follow.