STAT 479 Lecture 13

Bradley-Terry Models II

Recap of Lecture 12

  • Bradley-Terry model: team \(j\) has latent strength \(\lambda_{j}\)
  • \(\textrm{Log-odds}(i \textrm{ beats } j) = \lambda_{i} - \lambda_{j}\) \[ \mathbb{P}(i \textrm{ beats } j ) = \frac{1}{1 + e^{-(\lambda_{i} - \lambda_{j})}} = \frac{e^{\lambda_{i}}}{e^{\lambda_{i}} + e^{\lambda_{j}}}. \]
  • 2024 NCAA D1 Women’s Hockey results: \(\hat{\lambda}_{\textrm{WIS}} \approx 4.68\) & \(\hat{\lambda}_{\textrm{OSU}} \approx 3.18\)
1/(1 + exp(-1 * (4.68 - 3.18)))
[1] 0.8175745

BT with Home Advantage I

  • Basic BT model: \(\mathbb{P}(i \textrm{ beats } j)\) does not depend on location
  • Elaborated model: introduce a latent intercept \(\lambda_{0}\) so that \[ \begin{align*} \mathbb{P}(\textrm{team i beats team j at i's home}) &= \frac{e^{\lambda_{i} + \lambda_{0}}}{e^{\lambda_{i} + \lambda_{0}} + e^{\lambda_{j}}} \\ & \\ &= \frac{1}{1 + e^{-1 \times (\lambda_{i} + \lambda_{0} - \lambda_{j})}} \end{align*} \]

BT w/ Home Advantage II

  • Probability that \(i\) beats \(j\) now depends on location
  • At \(i\)’s home: \(\frac{e^{\lambda_{0} + \lambda_{i}}}{e^{\lambda_{0} + \lambda_{i}} + e^{\lambda_{j}}}\)
  • At \(j\)’s home: \(\frac{e^{\lambda_{i}}}{e^{\lambda_{i}} + e^{\lambda_{j} + \lambda_{0}}}\)
  • At a neutral site: \(\frac{e^{\lambda_{i}}}{e^{\lambda_{i}} + e^{\lambda_{j}}}\)

Data Preparation I

  • Data scraped from USCHO does not record game location
  • Step 1: is game at neutral site
  • Manually impute using some heuristics
    • Most regular season games at listed home team’s home
    • Tournaments on case-by-case basis
  • See course notes for details

Data Preparation II

  • Add neutral site indicator
load("wd1hockey_regseason_2024_2025.RData")

no_ties <-
  no_ties |>
  dplyr::mutate(
    neutral = dplyr::case_when(
      grepl("IceBreaker", Notes) ~ 1,
      grepl("East/West", Notes) & Home != "Minnesota" ~ 1,
      grepl("Nutmeg", Notes) & Home != "Sacred Heart" ~ 1,
      grepl("Mayor", Notes) & Home == "Rensselaer" ~ 1,
      grepl("Smashville", Notes) ~ 1,
      grepl("Beanpot", Notes) & Home != "Northeastern" ~ 1, 
      Date == "3/7/2025" & Home == "St. Lawrence" & Opponent == "Colgate" ~ 1,
      Date == "3/7/2025" & Home == "Minnesota" & Opponent == "Ohio State" ~ 1,
      Date == "3/8/2025" & Home == "Minnesota" & Opponent == "Wisconsin" ~ 1,
      .default = 0))

Data Preparation III

  • BradleyTerry2 expects data in following format
    • One row per combination of home team, away team, and location
    • Columns counting home team and away team wins
  • New columns home.team and away.team
    • These variables are themselves data frames
    • Columns recording team identity & whether they are at home
unik_teams <- sort(unique(c(no_ties$Home, no_ties$Opponent)))

results <-
  no_ties |>
  dplyr::rename(home.team = Home, away.team = Opponent) |>
  dplyr::group_by(home.team, away.team, neutral) |> 
  dplyr::summarise(
    home.win = sum(Home_Winner), 
    away.win = sum(Opp_Winner), .groups = 'drop') |> 
  dplyr::mutate(
    home.team = factor(home.team, levels = unik_teams), 
    away.team = factor(away.team,levels = unik_teams),
    home.athome = ifelse(neutral == 1, 0, 1),
    away.athome = 0)

Model Fitting

tmp_df <- data.frame(home.win = results$home.win, away.win = results$away.win)
tmp_df$home.team <- data.frame(team = results$home.team, at.home = results$home.athome)
tmp_df$away.team <- data.frame(team = results$away.team, at.home = results$away.athome)
fit <-
  BradleyTerry2::BTm(
    outcome = cbind(home.win, away.win),
    player1 = home.team, player2 = away.team,
    formula = ~ team + at.home,
    refcat = "New Hampshire",
    id = "team",
    data = tmp_df) 
summary(fit)

Call:
BradleyTerry2::BTm(outcome = cbind(home.win, away.win), player1 = home.team, 
    player2 = away.team, formula = ~team + at.home, id = "team", 
    refcat = "New Hampshire", data = tmp_df)

Coefficients:
                      Estimate Std. Error z value Pr(>|z|)    
teamAssumption        -4.74585    1.24679  -3.806 0.000141 ***
teamBemidji State     -0.75271    0.80275  -0.938 0.348415    
teamBoston College     0.88528    0.52440   1.688 0.091375 .  
teamBoston University  1.13510    0.52380   2.167 0.030229 *  
teamBrown              0.14970    0.65910   0.227 0.820326    
teamClarkson           1.22543    0.62242   1.969 0.048975 *  
teamColgate            2.07470    0.65875   3.149 0.001636 ** 
teamConnecticut        1.11779    0.53080   2.106 0.035217 *  
teamCornell            2.63595    0.71812   3.671 0.000242 ***
teamDartmouth         -0.79392    0.69582  -1.141 0.253875    
teamFranklin Pierce   -3.94938    1.23673  -3.193 0.001406 ** 
teamHarvard           -2.17071    0.87974  -2.467 0.013608 *  
teamHoly Cross        -0.50866    0.52815  -0.963 0.335499    
teamLindenwood        -1.47362    0.71219  -2.069 0.038532 *  
teamLIU               -3.07687    1.21642  -2.529 0.011425 *  
teamMaine             -0.38395    0.53066  -0.724 0.469350    
teamMercyhurst         0.51099    0.64323   0.794 0.426955    
teamMerrimack         -0.78126    0.52339  -1.493 0.135519    
teamMinnesota          2.76577    0.74846   3.695 0.000220 ***
teamMinnesota Duluth   2.17957    0.76679   2.842 0.004477 ** 
teamMinnesota State    0.78237    0.78409   0.998 0.318375    
teamNortheastern       1.02645    0.51130   2.008 0.044695 *  
teamOhio State         3.27350    0.79365   4.125 3.71e-05 ***
teamPenn State         2.02649    0.71347   2.840 0.004507 ** 
teamPost              -4.42233    1.23567  -3.579 0.000345 ***
teamPrinceton          0.90657    0.66577   1.362 0.173298    
teamProvidence         0.66011    0.52149   1.266 0.205585    
teamQuinnipiac         1.07109    0.59894   1.788 0.073726 .  
teamRensselaer        -0.29970    0.61486  -0.487 0.625954    
teamRIT               -0.21963    0.66882  -0.328 0.742616    
teamRobert Morris     -1.61446    0.71727  -2.251 0.024397 *  
teamSacred Heart      -3.33864    1.22272  -2.731 0.006324 ** 
teamSt. Anselm        -4.02309    1.20138  -3.349 0.000812 ***
teamSt. Cloud State    1.28838    0.75363   1.710 0.087349 .  
teamSt. Lawrence       1.31564    0.62136   2.117 0.034229 *  
teamSt. Michael's     -5.91267    1.29972  -4.549 5.39e-06 ***
teamSt. Thomas         0.08436    0.80906   0.104 0.916954    
teamStonehill         -4.17369    1.23819  -3.371 0.000749 ***
teamSyracuse          -0.39362    0.65720  -0.599 0.549213    
teamUnion             -0.05669    0.59147  -0.096 0.923646    
teamVermont           -0.56482    0.52659  -1.073 0.283454    
teamWisconsin          4.77395    0.99956   4.776 1.79e-06 ***
teamYale               0.82238    0.63286   1.299 0.193779    
at.home                0.31079    0.09542   3.257 0.001126 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 817.92  on 476  degrees of freedom
Residual deviance: 455.43  on 432  degrees of freedom
AIC: 661.38

Number of Fisher Scoring iterations: 6

Fitted Probabilities

lambda0_hat <- coef(fit)["at.home"] 
lambda_hat <- BradleyTerry2::BTabilities(fit)
cat("Intercept: ", lambda0_hat, "\n")
Intercept:  0.3107851 
cat("WIS: ", lambda_hat["Wisconsin", "ability"], "\n")
WIS:  4.773946 
cat("OSU: ", lambda_hat["Ohio State", "ability"], "\n")
OSU:  3.273501
  • Log-odds that WIS beats OSU:
    • At WIS: \(\hat{\lambda}_{\textrm{WIS}} + \hat{\lambda}_{0} - \hat{\lambda}_{\textrm{OSU}}.\)
    • At OSU: \(\hat{\lambda}_{\textrm{WIS}} - (\hat{\lambda_{\textrm{OSU}}} + \hat{\lambda}_{0}).\)
    • Neutral site: \(\hat{\lambda}_{\textrm{WIS}} - \hat{\lambda}_{\textrm{OSU}}.\)
  • diff: quantifies different in team strengths (\(\textrm{WIS} - \textrm{OSU}\))
diff <- lambda_hat["Wisconsin","ability"] - lambda_hat["Ohio State", "ability"]
cat("At WIS:", round(100 * 1/(1 + exp(-1 * (diff + lambda0_hat))), digits = 2), "%\n")
At WIS: 85.95 %
cat("AT OSU:", round(100 * 1/(1 + exp(-1 * (diff - lambda0_hat))), digits = 2), "%\n")
AT OSU: 76.67 %
cat("AT neutral site:", round(100 * 1/(1 + exp(-1 * (diff))), digits = 2), "%\n")
AT neutral site: 81.76 %

Frozen 4 Simulation

Background

  • Single-elimination tournament w/ last four teams
  • WIS, OSU, COR, and MIN
  • First semi-final (SF1): WIS vs MIN
  • Second semi-final (SF2): OSU vs COR
  • Finals: winner SF1 vs winner SF2
  • All games played at Ridder Arena (MIN’s home ice)

Computing Matchup Probabilities

  • diff: difference in team strengths (log-odds scale)
  • prob: \(\mathbb{P}(\textrm{higher seed beats lower seed})\)
seeds <- data.frame(
  Team = c("Wisconsin", "Ohio State", "Cornell", "Minnesota"),
  Seed = 1:4)

possible_matchups <-
  expand.grid(Hi = seeds$Team, Lo = seeds$Team) |>
  as.data.frame() |>
  dplyr::inner_join(y = seeds |> dplyr::rename(Hi = Team, Hi.Seed=Seed), by = "Hi") |>
  dplyr::inner_join(y = seeds |> dplyr::rename(Lo = Team, Lo.Seed=Seed), by = "Lo") |>
  dplyr::filter(Hi.Seed < Lo.Seed) |>
  dplyr::mutate(neutral = ifelse(Hi == "Minnesota" | Lo == "Minnesota", 0, 1)) |>
  dplyr::mutate(
    diff = lambda_hat[Hi, "ability"] - lambda_hat[Lo, "ability"],
    prob = dplyr::case_when(
      neutral == 1 ~ 1/(1 + exp(-1 * diff)),
      neutral == 0 & Hi == "Minnesota" ~ 1/(1 + exp(-1 * (diff + lambda0_hat))),
      neutral == 0 & Lo == "Minnesota" ~ 1/(1 + exp(-1 * (diff - lambda0_hat)))))

Matchup Probabilities

          Hi         Lo Hi.Seed Lo.Seed neutral       diff      prob
1  Wisconsin Ohio State       1       2       1  1.5004449 0.8176408
2  Wisconsin    Cornell       1       3       1  2.1379945 0.8945416
3 Ohio State    Cornell       2       3       1  0.6375496 0.6541993
4  Wisconsin  Minnesota       1       4       0  2.0081754 0.8451936
5 Ohio State  Minnesota       2       4       0  0.5077305 0.5490778
6    Cornell  Minnesota       3       4       0 -0.1298191 0.3915970
semis <- 
  data.frame(Hi = c("Wisconsin", "Ohio State"), Lo = c("Minnesota", "Cornell")) |>
  dplyr::left_join(possible_matchups |> dplyr::select(Hi, Lo, prob), by = c("Hi", "Lo"))
semis
          Hi        Lo      prob
1  Wisconsin Minnesota 0.8451936
2 Ohio State   Cornell 0.6541993

Semi-Final Simulation

semis <- 
  data.frame(Hi = c("Wisconsin", "Ohio State"), Lo = c("Minnesota", "Cornell")) |>
  dplyr::left_join(possible_matchups |> dplyr::select(Hi, Lo, prob), by = c("Hi", "Lo"))

set.seed(479)
sf_winners <- c(NA, NA)

sf_outcomes <- rbinom(n = nrow(semis), size = 1, prob = semis$prob) 
for(i in 1:nrow(semis)){
  if(sf_outcomes[i] == 1) sf_winners[i] <- semis$Hi[i] 
  else sf_winners[i] <- semis$Lo[i] 
}
cat("Semi-final outcomes:", sf_outcomes, "\n")
Semi-final outcomes: 0 1 
cat("Semi-final winners:", sf_winners, "\n")
Semi-final winners: Minnesota Ohio State

Finals Simulation

  • Extract semi-final winners
  • Get matchup probability for finals
finals <- 
  data.frame(Team1 = sf_winners[1], Team2 = sf_winners[2]) |>
  dplyr::left_join(y = seeds |> dplyr::rename(Team1 = Team, Team1.Seed = Seed), by = "Team1") |>
  dplyr::left_join(y = seeds |> dplyr::rename(Team2 = Team, Team2.Seed = Seed), by = "Team2") |>
  dplyr::mutate(
    Hi = ifelse(Team1.Seed < Team2.Seed, Team1, Team2),
    Lo = ifelse(Team1.Seed < Team2.Seed, Team2, Team1)) |>
  dplyr::select(Hi, Lo) |> 
  dplyr::left_join(
    y = possible_matchups |> dplyr::select(Hi, Lo, prob), by = c("Hi", "Lo"))

finals
          Hi        Lo      prob
1 Ohio State Minnesota 0.5490778
final_outcome <- rbinom(n = 1, size = 1, prob = finals$prob)
if(final_outcome == 1){
  final_winner <- finals$Hi[1]
} else{
  final_winner <- finals$Lo[1]
}

cat("Finals outcome:", final_outcome, "\n")
Finals outcome: 1 
cat("Finals winner:", final_winner, "\n")
Finals winner: Ohio State

Repeated Simulations

  • Embed earlier code into large for loop
  • In each iteration, we will save
    • Semi-finals winners:SF1_winner and SF2_winner
    • Teams in finals: Finals_Hi and Finals_Lo
    • Eventual Champion Champion
  • Using simulations, we will estimate
    • \(\mathbb{P}(\textrm{WIS makes it to and wins the finals})\)
    • \(\mathbb{P}(\textrm{WIS wins the finals} \vert \textrm{Wisconsin makes it to the finals}).\)

Simulation Snapshot

  • Results of the first 10 simulations
   SF1_Winner SF2_Winner  Finals_Hi  Finals_Lo   Champion
1   Wisconsin Ohio State  Wisconsin Ohio State  Wisconsin
2   Wisconsin Ohio State  Wisconsin Ohio State Ohio State
3   Wisconsin    Cornell  Wisconsin    Cornell  Wisconsin
4   Wisconsin Ohio State  Wisconsin Ohio State  Wisconsin
5   Wisconsin    Cornell  Wisconsin    Cornell  Wisconsin
6   Wisconsin Ohio State  Wisconsin Ohio State  Wisconsin
7   Wisconsin Ohio State  Wisconsin Ohio State  Wisconsin
8   Wisconsin Ohio State  Wisconsin Ohio State  Wisconsin
9   Wisconsin    Cornell  Wisconsin    Cornell  Wisconsin
10  Minnesota Ohio State Ohio State  Minnesota  Minnesota
table(results$Champion)

   Cornell  Minnesota Ohio State  Wisconsin 
       502        765       1544       7189

Estimating Complex Probabilities

  • Across the 10,000 simulations:
    • WIS made the finals 8496 times
    • WIS made and won the finals 7189 times
  • \(\mathbb{P}(\textrm{WIS wins the finals} \vert \textrm{WIS makes it to the finals}) \approx 84.6\%.\)
sum(results$SF1_Winner == "Wisconsin" & 
      results$Champion == "Wisconsin")/sum(results$SF1_Winner == "Wisconsin")
[1] 0.8461629

Alternative Rules

  • What if higher seed always gets to play at home?
  • See lecture notes for simulation
  • Very similar code but uses different matchup probabilities
  • Under alternative, WIS wins championship about 77% of the time (as opposed to 71% at neutral site)

Introducing Covariates

  • What if we had team-level covariates?
    • Offensive/defensive rating, record to date, other measures of strength
  • Possible to elaborate model: \(\lambda_{j} = U_{j} + \boldsymbol{\mathbf{x}}_{j}^{\top}\beta\)
  • Under elaborated model \[ \log\left(\frac{\mathbb{P}(i \textrm{ beats } j)}{\mathbb{P}(j \textrm{ beats } i)}\right) = U_{i} - U_{j} + (\boldsymbol{\mathbf{x}}_{i} - \boldsymbol{\mathbf{x}}_{j})^{\top}\beta \]
    • Depends on difference in observed covariates & latent team-specific random effect