You might be reading other sites and see lots of posts on ELO models this post isn’t about building an ELO model yourself that’s another post. This post is about how to apply Bradley Terry models for AFL

# Why Bother?

Each of the 18 footy teams plays 22 games a year, we can see who has played who and who won each game. We could just use the AFL ladder but each team doesn’t play each other equally some teams have have easier schedules than others. What strength of shedule means is that we can’t use the ladder as the best measure of teams as some times might have a worse win/loss record but simply have had a tougher draw than opposing teams who might have had it a lot easier.

A model based approach can help address this problem.

A simple Bradley Terry model treats outcomes of games as an independent Bernoulli random variable with a Bernoulli distribution $$p_{ij}$$

The log odds corresponding to the probability $$p_{ij}$$ that team i beats team j is

$$log\frac{p_{ij}}{1-p_{ij}} = \beta{i} - \beta{j}$$

The problem is that this is over-parametized which means that its exactly the same if we were to add a fixed constant to all the values of $$\beta{i}$$

The model as it stands now doesn’t have home ground advantage. We can incorporate that by including an intercept term $$\alpha$$

$$log\frac{p_{ij}}{1-p_{ij}} = \alpha + \beta{i} - \beta{j}$$

By rearranging the equation we can see how it increases the log-odds of the home team winning by a constant $$\alpha$$

library(fitzRoy)
library(tidyverse)
## -- Attaching packages -------------------------------- tidyverse 1.2.1 --
## v ggplot2 2.2.1     v purrr   0.2.5
## v tibble  1.4.2     v dplyr   0.7.5
## v tidyr   0.8.1     v stringr 1.3.1
## v readr   1.1.1     v forcats 0.3.0
## -- Conflicts ----------------------------------- tidyverse_conflicts() --
## x dplyr::lag()    masks stats::lag()
df<-fitzRoy::get_match_results()
num_teams=18
df1<-df%>%filter(Season==2017)%>%
mutate(Y=if_else(Home.Points>Away.Points, 1, 0))

teams=unique(df1$Home.Team) df1$Home.Team[df1$Home.Team == "Richmond"] <- 1 df1$Home.Team[df1$Home.Team == "Essendon"] <- 2 df1$Home.Team[df1$Home.Team == "Port Adelaide"] <- 3 df1$Home.Team[df1$Home.Team == "Hawthorn"] <- 4 df1$Home.Team[df1$Home.Team == "Gold Coast"] <- 5 df1$Home.Team[df1$Home.Team == "GWS"] <- 6 df1$Home.Team[df1$Home.Team == "Melbourne"] <- 7 df1$Home.Team[df1$Home.Team == "West Coast"] <- 8 df1$Home.Team[df1$Home.Team == "Adelaide"] <- 9 df1$Home.Team[df1$Home.Team == "North Melbourne"] <- 10 df1$Home.Team[df1$Home.Team == "Carlton"] <- 11 df1$Home.Team[df1$Home.Team == "Collingwood"] <- 12 df1$Home.Team[df1$Home.Team == "Brisbane Lions"] <- 13 df1$Home.Team[df1$Home.Team == "Fremantle"] <- 14 df1$Home.Team[df1$Home.Team == "Footscray"] <- 15 df1$Home.Team[df1$Home.Team == "Sydney"] <- 16 df1$Home.Team[df1$Home.Team == "Geelong"] <- 17 df1$Home.Team[df1$Home.Team == "St Kilda"] <- 18 df1$Home.Team<-as.integer(df1$Home.Team) df1$Away.Team[df1$Away.Team == "Richmond"] <- 1 df1$Away.Team[df1$Away.Team == "Essendon"] <- 2 df1$Away.Team[df1$Away.Team == "Port Adelaide"] <- 3 df1$Away.Team[df1$Away.Team == "Hawthorn"] <- 4 df1$Away.Team[df1$Away.Team == "Gold Coast"] <- 5 df1$Away.Team[df1$Away.Team == "GWS"] <- 6 df1$Away.Team[df1$Away.Team == "Melbourne"] <- 7 df1$Away.Team[df1$Away.Team == "West Coast"] <- 8 df1$Away.Team[df1$Away.Team == "Adelaide"] <- 9 df1$Away.Team[df1$Away.Team == "North Melbourne"] <- 10 df1$Away.Team[df1$Away.Team == "Carlton"] <- 11 df1$Away.Team[df1$Away.Team == "Collingwood"] <- 12 df1$Away.Team[df1$Away.Team == "Brisbane Lions"] <- 13 df1$Away.Team[df1$Away.Team == "Fremantle"] <- 14 df1$Away.Team[df1$Away.Team == "Footscray"] <- 15 df1$Away.Team[df1$Away.Team == "Sydney"] <- 16 df1$Away.Team[df1$Away.Team == "Geelong"] <- 17 df1$Away.Team[df1$Away.Team == "St Kilda"] <- 18 df1$Away.Team<-as.integer(df1$Away.Team) loglik = function(theta, Home.Team, Away.Team, Y) { alpha = theta[1] beta = c(0, theta[-1]) params = alpha + beta[Home.Team] - beta[Away.Team] return(sum(Y * params - log(1 + exp(params)))) } theta0 = rep(0, num_teams) result = optim(theta0, loglik, Home=df1$Home.Team, Away=df1$Away.Team, Y=df1$Y,
method='BFGS', control=list('fnscale'=-1))

coefs = c(0, result$par[-1]) ranking = order(coefs, decreasing=TRUE) ranking ## [1] 9 1 17 6 16 3 8 7 18 2 15 4 12 14 11 10 5 13 # So what is our HGA using a Bradley Terry Model? result$par[1]
## [1] 0.4792202
exp(-result$par[1]) ## [1] 0.6192661 # Is HGA significant? loglik_noalpha = function(theta, Home, Away, Y) { beta = c(0, theta) params = beta[Home] - beta[Away] return(sum(Y * params - log(1 + exp(params)))) } theta0 = rep(0, num_teams - 1) result_noalpha = optim(theta0, loglik_noalpha, Home=df1$Home.Team, Away=df1$Away.Team, Y=df1$Y,
method='BFGS', control=list('fnscale'=-1))

print(result$value) ## [1] -116.2975 print(result_noalpha$value)
## [1] -120.8496
statistic = -2 * (result_noalpha$value - result$value)
p_value = 1 - pchisq(statistic, df=1)
print(statistic)
## [1] 9.104269
print(p_value)
## [1] 0.002550136

$$P(i beats j)=2/3, P(j beats k)= 2/3 then P(i beats k)=4/5$$