Rugby — R cross-check via lme4::glmer (crossed-effects Poisson GLMM)¶

Companion to rugby_python.ipynb / rugby_pymc.ipynb. The crossed attack/defence model is exactly a Poisson GLMM with two crossed random effects, which lme4::glmer fits by maximum likelihood (Laplace approximation) — the same frequentist engine used for Epil in hpois_R.ipynb, here with two grouping factors.

Key reparameterisation. In the model, defence enters as $-\text{def}_{\text{opp}}$. glmer only adds random effects, so we fit $\log\theta=\mu+\text{home}\cdot\mathbb1_{\text{home}}+\underbrace{u_{\text{scorer}}}_{\text{attack}}+\underbrace{v_{\text{opp}}}_{=-\text{def}}$ and recover $\text{def}_t=-v_t$. A team that is easy to score against has a large positive $v$ (poor defence), i.e. a negative $\text{def}$.

In [1]:
suppressMessages(library(lme4))
d <- read.csv('rugby.csv'); d <- d[d$year == 2014, ]
# stack to 30 scoring events: home-scoring rows then away-scoring rows
long <- data.frame(
  y       = c(d$home_score, d$away_score),
  is_home = c(rep(1, nrow(d)), rep(0, nrow(d))),
  scorer  = c(d$home_team, d$away_team),
  opp     = c(d$away_team, d$home_team)
)
fit <- glmer(y ~ is_home + (1 | scorer) + (1 | opp), data = long, family = poisson)
fe <- fixef(fit); vc <- as.data.frame(VarCorr(fit))
cat(sprintf("mu (intercept) = %.3f\n", fe[["(Intercept)"]]))
cat(sprintf("home (is_home) = %.3f   -> home advantage x%.2f\n", fe[["is_home"]], exp(fe[["is_home"]])))
cat(sprintf("sd_att (scorer) = %.3f   sd_def (opp) = %.3f\n",
            vc$sdcor[vc$grp=="scorer"], vc$sdcor[vc$grp=="opp"]))

mu (intercept) = 2.670

home (is_home) = 0.368   -> home advantage x1.44

sd_att (scorer) = 0.324   sd_def (opp) = 0.375
In [2]:
re  <- ranef(fit)
att <- re$scorer[, 1]; names(att) <- rownames(re$scorer)
defn <- -re$opp[, 1];  names(defn) <- rownames(re$opp)          # negate: opp effect = -defence
teams <- sort(unique(long$scorer))
tab <- data.frame(team = teams,
                  attack  = round(att[teams], 3),
                  defence = round(defn[teams], 3))
tab <- tab[order(-tab$attack), ]
print(tab, row.names = FALSE)

png('rugby_R.png', width=560, height=480, res=110)
plot(att[teams], defn[teams], pch=19, col='steelblue', cex=1.4,
     xlab='attack (BLUP)', ylab='defence  (= -opp effect)',
     main='glmer team-strength map (2014 Six Nations)')
abline(h=0, col='grey'); abline(v=0, col='grey')
text(att[teams], defn[teams], labels=teams, pos=4, cex=0.8)
invisible(dev.off()); cat('\nwrote rugby_R.png\n')

     team attack defence
  England  0.395   0.241
  Ireland  0.180   0.529
    Wales  0.127   0.134
   France  0.052  -0.135
    Italy -0.228  -0.513
 Scotland -0.502  -0.288

wrote rugby_R.png

Results¶

  • glmer recovers the same fixed effects as the two Bayesian engines: μ ≈ 2.67, home ≈ 0.37 (×1.44 home advantage). The variance components come out smaller — σ_att ≈ 0.32, σ_def ≈ 0.38 — than the Bayesian posterior means (0.47–0.62): with only 6 teams per grouping factor, maximum-likelihood variance estimates are biased downward (and have no prior to hold them up), the classic "few groups" GLMM issue. This shrinks the BLUPs a little more than the Bayesian posteriors.
  • Even so, the attack/defence BLUPs match the from-scratch and PyMC posterior means in both ordering and (closely) magnitude: England top attack, Ireland best defence, Scotland/Italy weak. (BLUPs are shrunk point estimates, so they sit a touch closer to zero, but the ranking is identical.)
  • Three engines, one crossed-effects model: from-scratch Metropolis-within-Gibbs ≈ PyMC NUTS ≈ R glmer (Laplace ML). The only modelling subtlety is the sign flip on the opponent effect — once defence is read as $-v_{\text{opp}}$, the GLMM is the rugby model.
  • Same caveat as the other notebooks: Poisson is a simplification of lumpy rugby scoring, used here because the example is about relative team strength, which it captures well.