Simulating the Bundesliga (Part 1)

Formulating the model

Working with football data is nice, because the number of goals often match a Poisson distribution pretty closely. I therefore fit a Poisson Generalized Linear Mixed Model (GLMM) with a log-link function, one binary predictor home := (side == "H"), and two grouping variables (team with index \(t\) and match_id with index \(m\)). The likelihood for our data (scored goals \(G\)) can be written like this:

\[ \begin{align} G_i &\sim \text{Poisson}(\mu_i) \newline \log \mu_i &= \alpha_{t[i]} + \beta_{t[i]} \text{home}_i + \delta_{m[i]} \end{align} \]

The coefficient \(\alpha\) is the log expected number of goals when \(\text{home}=0\), i.e. it measures the expected number of goals when team \(t\) plays not at home. When team \(t\) does play at home, i.e. \(\text{home}=1\), then \(\beta\) is the additional effect of playing at home, i.e. the home advantage. It is really only an advantage if \(\beta>0\), which is a very reasonable assumption – however, I don’t restrict \(\beta\) to be positive.

The \(\delta\) is allowed to vary for each match. This will capture unobserved factors which affect both teams (in the same way) in a match. Unobserved here means stuff that is not explicitly in the model, i.e. weather conditions, time of day etc., but I don’t expect \(\delta\) to be big. Not because the aren’t possibly any systematic effects, which are not part of the model, but rather because \(\delta\) is an amalgamation of many smaller (and likely counteracting) effects.

Since \(\alpha\) and \(\beta\) are allowed to vary by team and \(\delta\) varies by match, we need to add some structure to the model, otherwise we would be horribly overfitting. Here comes the hierarchical part of the model into play.

\[ \begin{align} \delta_m &\sim \text{Normal}(0, \sigma_{M}^2) \newline \begin{bmatrix} \alpha \newline \beta \end{bmatrix}_t &\sim \text{MultiNormal} \left( \begin{bmatrix} \mu_\alpha \newline \mu_\beta \end{bmatrix}, \Sigma \right) \end{align} \]

All \(\delta_m\) come from the same distribution with common standard deviation \(\sigma_M\). Here, \(\sigma_M\) is like a hyper-parameter, or tuning parameter which controls the amount of heterogeneity across matches. If \(\sigma_M\) is high, we’d expect a high number of very low scoring games, e.g. 0-0, and a high number of high scoring games, e.g 4-5. When \(\sigma_M\) is low, the match outcomes will follow the estimated expected goals per team more closely.

Expected goals per team are estimated by \(\alpha\) and \(\beta\) which come from a Multivariate Normal distribution. The mean parameters \(\mu_\alpha\) and \(\mu_\beta\) represent the average of team expected goals away from home (\(\mu_\alpha\)) and the average of team home advantages (\(\mu_\beta\)). These effects are allowed to be correlated through the covariance matrix \(\Sigma\).

Now we have to set priors for these hierarchical parameters.

\[ \begin{align} \mu_\alpha &\sim \text{Normal}(0.4, 1^2) \newline \mu_\beta &\sim \text{Normal}(0.3, 0.5^2) \newline \sigma_M &\sim \text{Exponential}(1) \newline \Sigma &= \begin{bmatrix} \sigma_\alpha & 0 \newline 0 & \sigma_\beta \end{bmatrix} \Omega \begin{bmatrix} \sigma_\alpha & 0 \newline 0 & \sigma_\beta \end{bmatrix} \newline \sigma_\alpha, \sigma_\beta & \sim \text{Exponential}(1) \newline \Omega &\sim \text{LKJ}(1) \end{align} \]

I chose a prior for \(\mu_\alpha\) which implies that the median of the team expected goals away from home is roughly \(\exp(0.4) \approx 1.5\). This means, I expect half of the teams score less than 1.5 goals on average away from home and half of the teams score more than 1.5 goals, when playing away. The prior on \(\mu_\beta\) follows a similar logic. It is set so that \(\exp(0.4 + 0.3) \approx 2\) roughly half of the teams playing at home score less than 2 goals on average and the other half scores more than 2 on average when at home.

All other priors are rstanarm defaults – I think in reality rstanarm does something a little more complex, but this should essentially be what’s going on: The standard deviations \(\sigma_M, \sigma_\alpha, \sigma_\beta\) control the variability of the corresponding hierarchical parameters \(\delta, \alpha, \beta\). Their Exponential(1) priors are most likely too wide, but it’s better to err on the side of flexibility here.

The correlation matrix \(\Omega\) is given a LKJ(1) prior which is uniform in \((-1,1)\). The correlation between \(\alpha\) and \(\beta\) could be positive, since I’d reckon that good teams (high \(\alpha\)) play even better at home (high \(\beta\)). On the other hand, bad teams (low \(\alpha\)) could follow the strategy to just play good at home (high \(\beta\)) and don’t care so much about winning away matches which would imply a negative correlation. I guess the uniform prior is fine.

Running the model

We specified a probabilistic model, with a likelihood and priors. In Bayesian statistics we combine the likelihood of the data and our prior beliefs and compute the posterior distribution via Bayes rule.

\[p(\theta|D) = \frac{p(D|\theta) \times p(\theta)}{p(D)},\] or

\[\text{posterior} = \frac{\text{likelihood} \times \text{prior}}{\text{evidence}}\]

The posterior is thus a summary of what we know and belief about the world. If you are not familiar with Bayes theorem, you can check out this awesome video by Grant Sanderson.

Computing the posterior of non-trivial models is often extremely hard, if not impossible. But luckily we can use computers and well crafted programs to help us.

Now that we have formulated the model, let’s fit it.

Bayesian GLMMs can easily be run using the rstanarm::stan_glmer() function. The rstanarm package aims at offering Bayesian alternatives to frequently used frequentest methods – no pun intended. In this case it mimics the glmer function from the lme4 package, so we can use the lme4 formula syntax. rstanarm is build on Stan which is a probabilistic programming language for very efficient Markov Chain Monte Carlo (MCMC) – a way to produce samples from the posterior distribution of our model. Before we run the model we set the options to tell R that we have multiple cores available for MCMC.

options(mc.cores = parallel::detectCores())

poisson_glmm <- stan_glmer(
  G ~ side + (side|team) + (1|match_id),
  data = buli1920_long_obs, 
  family = poisson(link = "log"),
  prior_intercept = normal(0.4, 1, autoscale = FALSE),
  prior = normal(0.3, 0.5, autoscale = FALSE),
  iter = 2000,
  refresh = 500)

The model ran without problems in about 20 seconds on my laptop and produced 4 Markov Chains with 1,000 draws from the posterior each. So for all parameters in our model we have 4,000 plausible values. Let’s see the results.

Results

poisson_glmm

stan_glmer
 family:       poisson [log]
 formula:      G ~ side + (side | team) + (1 | match_id)
 observations: 448
------
            Median MAD_SD
(Intercept) 0.4    0.1   
sideH       0.1    0.1   

Error terms:
 Groups   Name        Std.Dev. Corr
 match_id (Intercept) 0.076        
 team     (Intercept) 0.247        
          sideH       0.199    0.19
Num. levels: match_id 224, team 18 

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

This summary should be looking familiar to anyone, who used the lme4 package before. It shows summary statistics of the posterior distribution for all parameters.

The posterior median for \(\alpha\) is 0.4. That’s really close to my prior belief! I’m a bit surprised – although, I probably shouldn’t be, but that’s a different and more philosophical debate. The posterior median for \(\beta\) is quite a bit lower than my prior belief. \(\beta\) is mostly positive though, so there is a measurable home advantage.

Match heterogeneity denoted by \(\delta\) in our model has a very low spread with the posterior median of \(\sigma_M\) being roughly 0.08. Variability across teams in \(\sigma_\alpha\) (median = 0.25) and \(\sigma_\beta\) (median = 0.20) is more prominent. However, I think it’s not surprising to find team differences.

Finally, we see that the correlation between scoring goals (away) and home advantage is positive (0.16): Better teams are expected to have an even bigger home advantage. This can also be seen when we plot the median expected goals for home and way matches of each team implied by the model.

If a team lies on the dashed diagonal line, it’s expected goals will be the same in home and away matches. Teams that lie above the dashed line are expected to score more goals at home than away.

The home advantage is seen for nearly every team (except for Werder Bremen, who really struggled at home this season). We can also see that the cloud of points is tilted slightly above the 45 degree line. This is due to the positive correlation between \(\alpha_t\) and \(\beta_t\) in the model above: Dortmund, Bayern Munich, and RB Leipzig are expected to score more away goals than other teams and even more home goals relative to other teams. Low scoring teams such as Fortuna Düsseldorf, Werder Bremen, and Paderborn on the other hand do not have a strong home advantage – Werder Bremen even has a slight home disadvantage.

Since we estimated a full Bayesian model, we have complete quantification of uncertainty for all the parameters in the model. Instead of looking just at the medians, we can inspect the whole posterior expected goal distribution for each team.

This second plot draws contour lines for each team separately to avoid cluttering. If you are not familiar with contour lines, think of them as the height indications on a map, where the higher you go, the more probable is the location you are standing in. Location is just a specific combination of expected away goals and expected home goals and the “height” indicates how plausible this location is.

Here we can see differences between teams in a bit more detail.

For example, Bayern Munich and Borrusia Dortmund both score a lot of goals. Dortmund has the stronger home advantage, while Bayern is better at scoring away from home.

Comparing Gladbach and Leipzig we see that Leipzig tends to score more goals, but Gladbach has the clearer home advantage and scores a bit more consistently, which you can see in the slightly tighter contour lines.

Taking a look at the teams at the bottom of the table, Paderborn and Werder Bremen have comparable away goal expectations, while Bremen has a considerably lower range of plausible expected home goals, indicated by the contour lines which fall below the dashed diagonal line.