Sometimes it feels a bit silly when a simple statistical model has a fancy-sounding name. But it also feels good to drop the following in casual conversation: “Ah, then I recommend a Plackett-Luce model, a straightforward generalization of the Bradley–Terry model, you know”, when a friend wonders how they could model their, say, pinball championship dataset. Incidentally, in this post we’re going to model the result of the IFPA 18 World Pinball Championship using a Plackett-Luce model, implemented in Stan as a generalization of the Bradley–Terry model, you know.

I know neither who Bradley, Terry, Plackett, nor Luce were, but I know when their models could be useful:

• A Bradley-Terry model can be used to model data where you have pairwise comparisons between different items, and you are interested in the underlying “fitness” of the items. A concrete example is sports where each match is a “pairwise comparison” between two players or teams, and you assume each player or team has an underlying skill or ability.
• A Plackett-Luce model can be useful when you have several rankings between items and you’re, again, interested in the “fitness” of each item. This model could be used to assess the quality of different products when each participant has ranked the items from best to worst. Or, in a sports setting, it can be used to model the underlying skills of each player when the outcome isn’t wins or losses, but rankings. Just like you have in pinball championships.

So, both models can model the skills of a number of players/teams, the only difference being that a Bradley-Terry works with wins/losses and a Plackett-Luce model works with rankings (1st/2nd/3rd/etc).

Now we’re going to grab some data with rankings from the IFPA 18 World Pinball Championship, implement a Bayesian Plackett-Luce model in Stan, and then take it for a spin.

IFPA 18 World Pinball Championship dataset

Despite the name, the IFPA 18 World Pinball Championship took place in 2023. The IFPA Championship is generally considered the most prestigious pinball competition, but the main reason why we’re going to analyze it here is because the results of all the 480 pinball matches that went down between the 80 competing players are available in a single spreadsheet! It’s not a very tidy dataset, however, and it will need some tidying up. I won’t bore you with the details, unless you really want to know:

library(readxl)
library(tidyverse)

# This data is NOT in a tidy format, and so the code to tidy it up
# will also be fairly messy...
# The result of IFPA 18 World Pinball Championship was downloaded 
# from here: https://www.ifpapinball.com/ifpa18/
ifpa_xlsx_path <- "IFPA 18 World Pinball Championship live results.xlsx"

# Reading and tidying the sheet with the pinball machine names used in the tournament
machines_wide <- read_xlsx(ifpa_xlsx_path, sheet = "Machines", range = "B2:F22")
machines_long <- machines_wide |> 
  select(old = OLD, mid = MID, new = NEW) |> 
  pivot_longer(
    cols = everything(), names_to = "game_category", values_to = "game_name"
  ) |> 
  # The two machines were missing from the original data
  add_row(game_category = "old", game_name = "Dodge City") |>
  add_row(game_category = "old", game_name = "8 Ball") |>
  arrange(game_name)

# Reading in the match data which is spread over multiple Session sheets, 
# each containing several sub-tables with the results of the games.
ifpa_xlsx_info <- 
  expand_grid(
    session = paste("Session", 1:8),
    tibble(
      group = 1:20,
      group_pos = seq(1, 96, 5)
    )
  ) |> 
  rowwise() |> 
  mutate(session_df = list(
    read_xlsx(
      ifpa_xlsx_path, sheet = session, 
      range = paste0("A", group_pos, ":K", group_pos + 4),
      col_names = FALSE
    )
  )) |> 
  ungroup()

# A function to pivot the sub-tables into a tidy data frame
pivot_session_df <- function(s) {
  player_names <- s$...1[2:5]
  game1 <- s$...3[2:5]
  score1 <- s$...5[2:5]
  game2 <- s$...6[2:5]
  score2 <- s$...8[2:5]
  game3 <- s$...9[2:5]
  score3 <- s$...11[2:5]
  
  tibble(
    player_name = rep(player_names, 3),
    round = rep(1:3, each = 4),
    game_name = c(game1, game2, game3),
    score = c(score1, score2, score3)
  )  
}

# Stitching the data together, and turning it into long format.
# Also, adding unique numerical identifiers for everything 
# as it's going to make things easier when writing the Stan model.
match_results <- 
  ifpa_xlsx_info |> 
  mutate(session_df = map(session_df, pivot_session_df)) |>
  unnest(session_df) |> 
  left_join( machines_long, by = "game_name") |>
  mutate(
    session = str_remove(session, "Session "),
    player_id = as.integer(factor(player_name)),
    game_id = as.integer(factor(game_name)),
    game_category_id = recode(game_category, `old` = 1, `mid` = 2, `new` = 3),
    rank = recode(score, `7` = 1, `5` = 2, `3` = 3, `1` = 4),
  ) |> 
  group_by(session, group, round) |> 
  mutate(round_id = cur_group_id()) |>
  select(
    session, group, round, round_id, player_name, player_id, score, rank,
    game_name, game_id, game_category, game_category_id
  )

write_csv(match_results, "ifpa-18-world-pinball-championship-match-results.csv")

The final tidy IFPA 18 Pinball Championship dataset includes the results from each of the eight qualifying sessions before the final tournament:

library(tidyverse)
match_results <- read_csv("ifpa-18-world-pinball-championship-match-results.csv")
match_results

# A tibble: 1,920 × 12
   session group round round_id player_name    player_id score  rank game_name  
     <dbl> <dbl> <dbl>    <dbl> <chr>              <dbl> <dbl> <dbl> <chr>      
 1       1     1     1        1 Michael Trepp         57     1     4 Little Joe 
 2       1     1     1        1 Bob Matthews          12     7     1 Little Joe 
 3       1     1     1        1 Mark Pearson          54     3     3 Little Joe 
 4       1     1     1        1 Escher Lefkoff        28     5     2 Little Joe 
 5       1     1     2        2 Michael Trepp         57     5     2 Jokerz     
 6       1     1     2        2 Bob Matthews          12     1     4 Jokerz     
 7       1     1     2        2 Mark Pearson          54     7     1 Jokerz     
 8       1     1     2        2 Escher Lefkoff        28     3     3 Jokerz     
 9       1     1     3        3 Michael Trepp         57     3     3 Indianapol…
10       1     1     3        3 Bob Matthews          12     1     4 Indianapol…
# ℹ 1,910 more rows
# ℹ 3 more variables: game_id <dbl>, game_category <chr>,
#   game_category_id <dbl>

In each round during these sessions, four players compete on the same pinball machine. The player that comes 1st gets 7 points, the second one gets 5 points, and so on. For example, looking at the first couple of rows above, we can see that when competing on the Little Joe (1972) machine, Bob Matthews won and Michael Trepp ended up last. During the pre-tournament sessions, players rotate to face off against most other players and compete on many different pinball machines. At the end of the sessions, each player’s score is tallied up, and the top 32 players proceed to the final tournament (not included in this dataset).

We’re now ready to model this data using a Plackett-Luce model!

The Plackett-Luce model

Despite the somewhat involved names, both the Bradley–Terry model and the Plackett-Luce model are fairly straightforward. In the simplest case, without covariates, both models assume that each player (or team/item/product) has an underlying skill. In themselves, these skill parameters don’t have any meaning. They only become meaningful when used to come up with the probability of each player winning. Let’s say n players compete, each with their own skill parameter skill_{1, 2, …, n}. Then the probability of each player winning is calculated as

$$p_1 = \frac{\exp({skill}_1)}{\sum(\exp({skill}_{1, 2, …, n}))}, \\[0.7em] p_2 = \frac{\exp({skill}_2)}{\sum(\exp({skill}_{1, 2, …, n}))}, \\[0.7em] … \\[0.5em] p_n = \frac{\exp({skill}_n)}{\sum(\exp({skill}_{1, 2, …, n}))}$$

The exp () makes sure that the result is always positive, even for negative skills, and by dividing by ∑(exp(skill_{1, 2, …, n})) we make the transformed skills sum to one like good probabilities should (this transformation is known as the softmax function). For example, if we have four players with skills c(0.8, 1.0, -1.0, 0.0) the probability distribution p over each player winning would be:

skills <- c(0.8, 1.0, -1.0, 0.0)
p <- exp(skills) / sum(exp(skills))

barplot(p, col = "salmon", ylab = "Probability of winning", names.arg = 
  paste(
    "Player", 1:4, "\nskill:", skills, "\nexp(skill):", round(exp(skills), 2),
    "\np: ", round(p, 2)
  )
)

The sole difference between the Bradley–Terry model and the Plackett-Luce model is in what data they model: Bradley-Terry models the winner of a competition between two players, while Plackett-Luce models the rankings in a competition with several players. It does this by assuming that the performance of each player isn’t influenced by the other players, which allows for modeling rankings as a series of competitions. Perhaps it’s easiest to explain by showing the “generative model” in R, for a competition with our four players from above:

players <- 1:4 # Our four competing players

# Sampling the winner of the bunch
p <- exp(skills[players]) / sum(exp(skills[players]))
first <- sample(players, prob = p, size = 1)

# Excluding the player who came first, who will win second place?
players <- players[players != first] 
p <- exp(skills[players]) / sum(exp(skills[players]))
second <- sample(players, prob = p, size = 1)

# Excluding the players who came first and second, who will win third place?
players <- players[players != second]
p <- exp(skills[players]) / sum(exp(skills[players]))
third <- sample(players, prob = p, size = 1)

# The player who's left gets fourth place
fourth <- players[players != third]

ranking <- c(first, second, third, fourth)
ranking

[1] 4 1 2 3

That’s basically the simple version of the Plackett-Luce model. In our case we, of course, know the rankings of all the competitions in the IFPA 18 Pinball Championship so we don’t want to simulate the data using fixed skill parameters. Instead, we’re now ready to do the Bayesian trick where we assume the player skills are unknown parameters and use the data to figure them out.

A Plackett-Luce model in Stan

While we could work with the tidied match_results in Stan, the model becomes easier to implement if we extract only the parts of the data that we need. Here, that’s the number of players (n_players), the number of rounds (n_round), and a matrix with the 1st, 2nd, 3rd, and 4th player_id for each round (player_ranks).

stan_data <- list(
  n_players = max(match_results$player_id),
  n_rounds = max(match_results$round_id),
  player_ranks = match_results |> 
    select(round_id, rank, player_id) |> 
    pivot_wider(names_from = rank, values_from = player_id) |> 
    arrange(round_id) |>
    select(`1`, `2`, `3`, `4`) |> 
    as.matrix()
)
head(stan_data$player_ranks)

      1  2  3  4
[1,] 12 28 54 57
[2,] 54 57 28 12
[3,] 28 54 57 12
[4,] 64 37  2 62
[5,] 64 62 37  2
[6,] 62  2 37 64

For example, looking at the first couple of rows in stan_data$player_ranks above, we can (again) see that in round one, on the first row, player_id: 12 (that is, Bob Matthews) won and player_id: 57 (Michael Trepp) came last.

Also, while we could implement the likelihood part of this model directly using sum()s and exp()s, a shortcut is to use the categorical_logit() distribution. The parameter to this distribution is a vector which will be softmax transformed (just what we want in the Plackett-Luce model) into the probability that each of the categories represented by this vector will be selected. For example, this sampling statement defines the likelihood that the first player, with skill 0.8 would win: 1 ~ categorical_logit([0.8, 1.0, -1.0, 0.0]).

Given that we’ve got our data nicely packaged up in stan_data and that we can use the categorical_logit() distribution, the Stan model definition becomes fairly straightforward:

plackett_luce_model_code <- "
data {
  int n_players; 
  int n_rounds;
  int player_ranks[n_rounds, 4];
}

parameters {
  vector[n_players] skills;
  real<lower=0> skills_sigma;
}

model {
  // A vector to hold the skills for each round
  // Just to limit the amunt of indexing we'll need to do
  vector[4] round_skills; 
  
  // It's important that the distribution over skills is anchored at a 
  // fixed value, here 0.0. Otherwise, as skills are relative to each other,
  // they wouldn't be identifiable.
  skills ~ normal(0, skills_sigma);
  skills_sigma ~ cauchy(0, 1);
  
  for (round_i in 1:n_rounds) {
    round_skills = skills[player_ranks[round_i]];
    
    // The likelihood of the winner winning out of the 4 players
    1 ~ categorical_logit(round_skills[1:4]);
    // The likelihood of the 2nd place winning out of the 3 remaining players
    1 ~ categorical_logit(round_skills[2:4]);
    // The likelihood of the 3rd place winning out of the 2 remaining players
    1 ~ categorical_logit(round_skills[3:4]);
    // And the remaining 4th place is guaranteed to 'win' against themselves...
  }
}
"

In the above model, I’ve placed a hierarchical distribution on the skills parameters. This works well for this dataset, where there are a lot of players and plenty of data. With less data, one might want to just put a fixed prior here (say, skills ~ normal(0, 1.0)). The model above is somewhat inflexible in that it only works for the case where there are exactly four players in each round, but I hope you can see that it’s straightforward to tweak it to allow for other number of players, as well.

Now we can finally go ahead and fit this model!

library(rstan)

model_fit <- stan(
  model_code = plackett_luce_model_code, data = stan_data, 
  chains = 4, iter = 20000, cores = 4
)

Skimping a bit on checking model convergence, we can, at least, see that both the trace plots and the number of effective samples look reasonable.

traceplot(model_fit, pars = c("skills_sigma", "skills[1]", "skills[2]", "skills[3]"))

model_fit |> 
  extract(pars = c("skills_sigma", "skills[1]", "skills[2]", "skills[3]"), permuted = FALSE) |> 
  monitor()

Inference for the input samples (4 chains: each with iter = 10000; warmup = 5000):

               Q5  Q50 Q95 Mean  SD  Rhat Bulk_ESS Tail_ESS
skills_sigma  0.2  0.3 0.4  0.3 0.1     1     2252     2456
skills[1]    -0.4  0.0 0.3  0.0 0.2     1    37109    14113
skills[2]    -0.5 -0.1 0.2 -0.1 0.2     1    29867    13139
skills[3]    -0.8 -0.4 0.0 -0.4 0.2     1     8267    13065

For each parameter, Bulk_ESS and Tail_ESS are crude measures of 
effective sample size for bulk and tail quantities respectively (an ESS > 100 
per chain is considered good), and Rhat is the potential scale reduction 
factor on rank normalized split chains (at convergence, Rhat <= 1.05).

That was a simple Plackett-Luce model in Stan. But finally, let’s have a look at the fitted Plackett-Luce model of the IFPA 18 World Pinball Championship.

The fitted IFPA 18 World Pinball Championship model

First, let’s extract the player skill parameters. For simplicity I’ll use the median point estimates and 95% probability intervals here, rather than working with the full distributions.

skills_summary <- summary(model_fit, pars = "skills")$summary

player_skill <- match_results |> 
  summarise(score = sum(score), .by = c(player_id, player_name)) |>
  arrange(player_id) |>
  mutate(
    median_skill = skills_summary[, "50%"],
    lower_skill = skills_summary[, "2.5%"],
    upper_skill = skills_summary[, "97.5%"]
  ) |> 
  arrange(desc(median_skill))
player_skill

# A tibble: 80 × 6
   player_id player_name              score median_skill lower_skill upper_skill
       <dbl> <chr>                    <dbl>        <dbl>       <dbl>       <dbl>
 1        39 Johannes Ostermeier        138        0.669      0.134        1.28 
 2        28 Escher Lefkoff             130        0.570      0.0811       1.15 
 3        54 Mark Pearson               114        0.362     -0.0724       0.864
 4        47 Keith Elwin                116        0.338     -0.0973       0.851
 5        21 Daniele Celestino Accia…   114        0.327     -0.111        0.835
 6        78 Viggo Löwgren              118        0.307     -0.119        0.804
 7        80 Zach Sharpe                116        0.294     -0.137        0.786
 8        36 Jason Zahler               106        0.240     -0.183        0.722
 9        48 Keri Wing                  108        0.230     -0.184        0.692
10        45 Josh Sharpe                110        0.226     -0.191        0.701
# ℹ 70 more rows

Not surprisingly, we find that Johannes Ostermeier, who ended up winning the whole tournament, also got the highest skill estimate. We can now also plot all the skill estimates:

player_skill |> 
  mutate(player_name = fct_reorder(player_name, median_skill)) |>
  ggplot(aes(x = median_skill, y = player_name)) +
  geom_point() +
  geom_errorbarh(aes(xmin = lower_skill, xmax = upper_skill), height = 0) +
  labs(
    title = "Player skill estimates for IFPA 18 World Pinball Championship",
    subtitle = "Posterior medians with 95% probablilty intervals",
    x = "Skill", y = ""
  ) +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 6))

It’s somewhat hard to interpret these skill estimates on their own, but two things to note here:

• Despite the large number of rounds played (480), the uncertainty in the skill parameters is large.
• Overall, players have very similar skill. This is not so surprising, as these are all top pinball players. If I had competed here, my skill estimate would be far far to the left.

To verify that these skill parameters make some sense, we could also plot them against each player’s final score from the qualifying sessions:

player_skill |> 
  ggplot(aes(x = median_skill, y = score)) +
  geom_point() +
  labs(
    title = "Player skill vs final score for IFPA 18 World Pinball Championship",
    x = "Median player skill", y = "Final Score"
  ) +
  theme_minimal()

The final score and the estimated skill are so strongly correlated that one might wonder if there was any point at all in going through the trouble of fitting a Plackett-Luce model. However, we can do much more interesting stuff with the skill estimates than we can with just the scores. For example, we can calculate the probability of players winning in different matchups. Say the top player played against the three players with the lowest skill:

highest_and_lowest_players <- union(
 slice_max(player_skill, median_skill, n = 1),
 slice_min(player_skill, median_skill, n = 3)
)
highest_and_lowest_players

# A tibble: 4 × 6
  player_id player_name         score median_skill lower_skill upper_skill
      <dbl> <chr>               <dbl>        <dbl>       <dbl>       <dbl>
1        39 Johannes Ostermeier   138        0.669       0.134      1.28  
2        77 Vid Kuklec             58       -0.574      -1.12      -0.107 
3        24 Didier Dujardin        56       -0.565      -1.11      -0.109 
4         7 Artur Natorski         62       -0.466      -0.979     -0.0334

skills <- highest_and_lowest_players$median_skill
p <- exp(skills) / sum(exp(skills))

barplot(p, col = "aquamarine", ylab = "Probability of winning", names.arg = 
  paste(highest_and_lowest_players$player_name, "\np: ", round(p, 2))
)

Here, according to this model, Johannes would have a 53% probability of winning when facing Vid, Didier, and Artur. Again, this shows that in a single game of pinball, even when the top player meets the lowest scoring players in the World Pinball Championship, it’s still far from guaranteed that the top player would win. But at the IFPA 18 World Pinball Championship Johannes did go all the way and won the final game as shown in this live stream.

All code for this post can be found in this Quarto markdown file.