Estimating home court advantage in Lithuanian Basketball League with Gaussian Processes

It was a bit more than a year ago that I built my first real Bayesian model. As part of Georgia Tech’s Bayesian Statistics class, I attempted to quantify home court advantage in Lithuanian Basketball League. I got that it is ~2.3-2.6 points, somewhat below what has been found in NBA. Since then, I learned quite a bit more about Bayesian modelling, and I thought I’ll revisit the problem. In case you’re interesting in the previous work or want to grab code that scrapes all the data I am using in this post, I suggest to check out a dedicated repo I have for that project.

As usual, my blog posts include all the code – but if you want to follow through a full notebook with everything, you can find it @ https://github.com/kamicollo/blog-posts.

Previously, I approached the problem by modelling the outcome of the game as a function of pair-wise team strength differences, plus the home court advantage of one of the teams. The simplest variation of the models I explored was:

$$ O^{s,k}_{i,j} \sim \text{Normal}(\alpha^{s}_{i,j} + β_i H^{k}_i – β_j H^{k}_j, \epsilon) \\ \alpha^{s}_{i,j} \sim \text{Normal}(0, 10) \\ \beta_{t} \sim \text{Normal}(4, 10) \\ \epsilon^{s} \sim \text{HalfCauchy}(5) \\ \\ $$

where $O^{k}_{i,j}$ is the score differential of game $k$ in season $s$ between teams $i$ and $j$, $\alpha^{s}_{i,j}$ is the pair-wise strength differential between teams $i$ and $j$ in season $s$, and $H^{k}_t$ is an indicator variable whether the team $t$ is playing at home. The main parameters of interest are $\beta_{t}$.

Looking back at the modelling approach, a couple of things immediatelly caught my eye:

• It’s interesting I chose to model pair-wise strengths instead of latent team strengths themselves. I wrote some justification for it in the accompanying report, but reading it today it doesn’t make so much sense. Especially given that I modelled them indepedently. Surely they should be correlated!
• The prior for $\beta$ is way too flat. It implies that the home court advantage is likely to take values as extreme as -15 to +20 points. I chose mean of $4$ based on prior research, but the combination with such a high variance wasn’t the smartest choice🤦‍♂️.
• I fit the above model on one season of data and also on multiple seasons, where I represented each team/season combination as a separate entity. I chose that because team strengths don’t stay constant throughout seasons – but looking back, treating them as completely independent was’t the best approach either.

In any way – let’s fit this model one more time (but with better priors) so that it’s available as a baseline. I’ll be using data across all seasons since 1995, excluding KMT cup games (the scraper I wrote last year still worked, so I downloaded a bit more data, too). Here’s how it looks like.

Let’s see if the model converged based on max $\hat{R}s$ across the parameter space. Let’s also visualize the home court advantage estimates (limited to teams that still exist as of 2022-2023 season).

Everything looks good from convergence perspective, but the home court advantage estimates vary a lot among teams. If you are familiar with Lithuanian basketball, you will also notice that they correlate with the strength of the teams (Žalgiris and Rytas are the two teams that have been battling for the champion title pretty much every season with a few exceptions). I suspect this model specification leads to some “spillover” effect from the team strength itself into the home court advantage estimates – although it could be that those teams have strongest fan support, too. What would other approaches yield?

Latent team strengths¶

I’ll start with something simple:

• Instead of modelling pairwise strengths (where only a few observations are available each season), I’ll simply model latent team strengths, and assume that the score differential is a function of team’s strengths’ differences.
• I’ll introduce a hierarchical structure to the home court advantages. I think the model may benefit from partial pooling, especially given the wide range of estimates seen in the previous specification.

The model could be written as (in practice, I use a non-centered parameterization for the hierarchical effects): $$ O^{s,k}_{i,j} \sim \text{Normal}(\alpha^{s}_{i} – \alpha^{s}_{j} + β_i H^{k}_i – β_j H^{k}_j, \epsilon) \\ \alpha^{s}_{i} \sim \text{Normal}(0, 5) \\ \bar{\beta} \sim \text{Normal}(4, 2) \\ γ \sim \text{HalfNormal}(2) \\ \beta_{t} \sim \text{Normal}(\bar{\beta}, γ) \\ \epsilon^{s} \sim \text{HalfCauchy}(5) \\ $$

This model also converges nicely. Looking at home court advantage estimates, we observe some regularization taking place, but it still attributes the largest effects to the strongest teams. Perhaps that’s indeed true – they have the largest arenas, after all.

Let’s compare it to the previous model using leave-one-out cross validation. Turns out, it’s actually better.

We also get team strengths – which is an interesting variable on its own. Let’s visualize them over time. Anyone familiar with Lithuanian basketball will recognize some real trends here. Because of how the model is specified, we also get estimates of team strengths in seasons they did not exist – that’s why some of the credible intervals are so wide.

Let’s specifically compare the strengths of two archrivals – Kauno Žalgiris and Vilniaus Rytas. This model isn’t very certain on who is stronger this year! Žalgiris is stronger, but the model is gives around 50% that the differential is less that ~3 points.

Improving the model with Gaussian Processes¶

Could we make this model even better? One obvious area of improvement is the fact that so far the model treats all seasons as independent observations. What if, instead, team’s strength was modelled as a function of season, but allowing for correlation between seasons? That’s a good use case for Gaussian Processes.

The full model specification is thus as follows:

$$ O^{s,k}_{i,j} \sim \text{Normal}(\alpha^{s}_{i} – \alpha^{s}_{j} + β_i H^{k}_i – β_j H^{k}_j, \epsilon) \\ \alpha^{s}_{t} \sim \text{GP}_t(0, f(s, l_t, \eta_t^{2})) \\ l_t \sim \frac{1}{2}Exp(1) \\ \eta_t^{2} \sim Exp(1) \\ \bar{\beta} \sim \text{Normal}(4, 2) \\ γ \sim \text{HalfNormal}(2) \\ \beta_{t} \sim \text{Normal}(\bar{\beta}, γ) \\ \epsilon^{s} \sim \text{HalfCauchy}(5) \\ $$

The one complication with this specification is that it means fitting one Gaussian Process per team. A naive for loop to specify them is not the way we want to go (I tried – compilation of the model alone took 18 minutes). An alternative would be to have just one GP, where teams are encoded in a one-of-k fashion. That would result in a (seasons x teams) x (teams + 1) dataset, though (the chosen specification is effectively seasons x teams). I tried that specification, too, but after running the sampling for over an hour, I gave up. The specifications used in this blog post, on the other hand, take just a few minutes to sample.

It should be possible to define such a model using pyMC’s GP module and leveraging pytensor.scan which is a “loop-like” functionality for such situations. But I could not get past some initial errors, and I did not want to scan everything – I wanted to see if I could solve at least some of the steps via broadcasting. Plus, it was educational in getting to understand pyMC’s GP internals.

It took a while and I had some really strong flashbacks of linaear algebra monsters from my Master’s, but I made it work!

First, we need an Exponential kernel that can handle 3D data:

Just to make sure it works correctly, we can compare results of it vs. the built-in pyMC class. Note that in the actual model fitting, I use a slightly different parameterization (controlled via pymc_compat parameter). For context on how they are different, see a forum post on pyMC Discourse.

And now we can use that in a pyMC model, replicating most of what the GP model does under the hood:

This model also samples pretty well.

Let’s see what latent team strengths we got this time. For comparison, I’ll overlay the strengths learned by the previous model, too.

I don’t know about you, but I find it really cool how GP’s work. With just 2 parameters per team ($ls$ and $\eta$), the model learned these sort of intricate shapes. I also find it fascinating what it learned. When it comes to teams that are constant title contenders (Žalgiris and Rytas), the functions learned are very smooth. The same applies to other teams that are often battling for 3-4th position – e.g. Neptūnas and Šiauliai. At the same time, it’s not that GP’s simply smoothened out what the previous model learned via independent parameters. For some teams – e.g. Lietkabelis or Labas Gas, the team strengths learned are nearly as wiggly as the ones learned previously.

What about home court advantages themselves? They are quite different, too. The strongest teams no longer have such dominant advantages; in fact, excluding Wolves (which is a new team established just this season), the highest mean homecourt advantage now belongs to Šiauliai.

Is this model actually better, from cross-validation perspective? It’s just a bit worse than the previous model, it turns out.

Gaussian Process with a latent team baseline strength¶

Could we make better? Let’s do one modification to the GP model. So far, the GP’s were zero-mean, i.e. the entire team strength function was fully dependent on the season. Let’s change that by modelling a latent “baseline” team strength which would be constant over time, while the GP will then capture variations around that baseline.

You could argue whether this assumption is a really good one – to me, it represents the fact that some teams are more established than others. It is also an opportunity to provide some informative priors to the model. Specifically, I set the priors as 10 for Žalgiris, 8 for Rytas, 5 for Lietkabelis, Šiauliai, Neptūnas and Juventus, while all the other teams’ strength is set to 0.

Based on ELPD, this model specification is the best so far!

And this is how team strengths look like over time (note that because of how the specification is constructed, the baseline team strength gets propagated backwarrds to seasons where the team did not exist).

What about home court advantages? In this specification, most of the teams have very similar levels. It got me concerned for a moment that the model was simply passing through the prior information and not learning anything, but the results below are robust to home court advantage prior – setting it it $\text{Normal}(0,2)$ leads to nearly identical results.

Here’s a different visualization – if we were trust the final model, Šiauliai commands the largest home court advantage, at nearly 4 points, followed by Rytas, Wolves and Lietkabelis, while the overall home court advantage across the league is somewhere between 3 and 4 points.

Finally, what do all the models think about Žalgiris vs. Rytas in this season? Well, zero falls into the 94% credible interval no matter which model you look at. Better grab that popcorn!

Addendum – strength variations within a season¶

Originally I collected this data by season. But as I kept on working through different model specifications, I realized it would be even more fun to model team strength variations over time. Players get injured, teams develop better understanding of each other – those are just a few factors that affect team strength within the season. Could we use Gaussian Processes to model that, too?

That’s the last model I tried fitting:

• Reduced only to data from 2019 and onwards (otherwise it takes a bit too long!)
• Each team no longer gets a baseline team strength, but instead – a baseline season strength (drawn from Normal, I didn’t want to go GP-in-GP route :D)
• Gaussian Process is fit for each team to capture variations across months

Here’s how the team strengths look like (this time, without credible intervals):

But is this model any good? Unfortunately – does not seem so. To understand if the model learned any true in-season variations, it’s useful to look at the lengthscale parameters. They are pretty much equal to priors – so the model did not really learn much in that regard. Perhaps a better way would be to model the strength as a function of two-dimensional input (season indicator + month indicator). Oh well – that’s for next time.