README.md
In opisthokonta/goalmodel: Prediction Models For Goals Scored
The goalmodel package let you build prediction models for the number of
goals scored in sport games. The models are primarily aimed at modelling
and predicting football (soccer) scores, but could also be applicable
for similar sports, such as hockey and handball.
If you find this package useful, please consider supporting the
development at https://ko-fi.com/opisthokonta
Installation
install.packages("devtools")
devtools::install_github("opisthokonta/goalmodel")
Whats new
Version 0.6.4 (development version)
- New function rho.ml to estimate the rho parameter of the Dixon-Coles
model.
- Fixed important bug that caused the defense parameters to be
slightly wrong in certain models. This affected mostly the defense
parameter for the first team in alphabetical order. In terms of
prediction error this affected the expected goals to be off by less
than 0.1 goals for games involving the first team. Models that were
affected were the default Poisson model and Gaussian models.
- Can now fit a hurdle model (See A. Owen (2017) The Application of
Hurdle Models to Accurately Model 0-0 Draws in Predictive Models of
Football Match Outcomes), but prediction methods are not yet
implemented.
- Small bugfixes in the
score_predictions function.
- README now has an example of how to use the
score_predictions
function.
Version 0.6
- New functions predict_btts() and pbtts() for computing
Both-Teams-To-Score probabilities.
- New function score_predictions(), for evaluating prediction accuracy
using Ranked Probability Scores, Brier scores, and Log scores.
- New function league_table() that computes a league table with points
and useful summary statistics.
Version 0.5
- New function expg_from_ou() - Finds the expected goals from
over/under probabilities.
- New function p1x2() - Compute win and draw probabilities from
expected goals.
- Bugfix: predict_goals() with return_df=TRUE did not work properly
when the team arguments were factors.
See NEWS.md for changes for complete version history.
Overview
The main function in this package is the goalmodel function, which is
used to fit a range of different models, including using different
probability distributions and predictors. It also allows for setting
parameters to a fixed value, which can be useful for more advanced
techniques such as two-step estimation.
The following probability distributions are available for fitting with
the goalmodel function via the model argument:
-
poisson, for the independent Poisson model, sometimes also called
the double Poisson model.
-
negbin, for the (independent) negative Binomial model, which
allows for more variance (overdisperison) than the Poisson model.
-
cmp, for the (independent) Conway-Maxwell-Poisson, which allows
for both underdispersion and overdispersion.
-
gaussian for the (independent) bivariate Gaussian (also known as
the normal distribution) model. This is useful for modelling scores
that might be decimal numbers.
In addition, the Dixon-Coles model (Dixon and Coles, 1997) is also
available. The model extend the Poisson model so that it modifies the
probabilities for low scores. The Dixon-Coles model can be fitted by
setting the dc argument to TRUE. Hence, the Dixon-Coles adjustment
can be used not only with Poisson models, but also with the Negative
Binomial and Conway-Maxwell-Poisson models.
With a fitted model, you would want to make predictions. This package
has functions for making different sorts of predictions:
-
predict_expg, for predicting expected goals.
-
predict_goals, for the scoreline probabilities.
-
predict_result, for the final result (home-draw-away or 1x2)
probabilities.
-
predict_ou, for over/under probabilities.
-
predict_btts, for both teams to score probabilities.
There are also functions that can be useful for making predictions
without relying on a fitted model, but instead uses expected goals
directly. You can use all of the probability distributions listed above
(except ‘gaussian’) as the underlying model, including the Dixon-Coles
adjustment.
-
p1x2
-
pbtts
There are a couple of functions for ‘reverse engineering’ results from
other models. These functions relies on the Poisson model.
-
expg_from_ou converts over/under probabilities to expected goals.
-
expg_from_probabilities converts 1x2 probabilities to expected
goals.
Misc
-
weights_dc, make weights to let results further back in time have
less influence on the estimates.
-
score_predictions computes ‘scoring rules’ to assess the quality
of predictions.
-
days_since_last_match calculates how long it was since the
previous game for each team.
-
matches_last_xdays calculates how many games each team have played
for a given period in the past.
There are also some functions to work directly with the
Conway-Maxwell-Poisson probability distribution.
-
dCMP, pCMP, and qCMP, are the probability, cumulative and
quantile functions for the CMP distribution.
-
eCMP and lambdaCMP converts between expected value and the
lambda parameter in the CMP distribution.
-
upsilon.ml estimates the dispersion parameter (upsilon) from data
and some given rate (lambda) or mean (mu) parameters.
The default model
The goalmodel package models the goal scoring intensity so that it is a
function of the two opposing teams attack and defense ratings. Let
λ1 and λ2) be the goal scoring intensities for the
two sides. The default model in the goalmodel package then models these
as
log(λ1) = μ + hfa + attack1 − defense2
log(λ2) = μ + attack2 − defense1
where μ is the intercept (approximately the average number of goals
scored) and hfa is the home field advantage given to team 1. By default
the number of goals scored, Y1 and Y2 are
distributed as
Y1 ∼ Poisson(λ1)
Y2 ∼ Poisson(λ2)
The default model can be modified in numerous ways as detailed in this
vignette.
Fitting the default model
Models are fitted with the goalmodel function. The minimum required
input are vectors containing the number of goals scored and the names of
the teams. Here I use data from the excellent
engsoccerdata
package.
library(goalmodel)
library(dplyr) # Useful for data manipulation.
library(engsoccerdata) # Some data.
library(Rcpp)
# Load data from English Premier League, 2011-12 season.
england %>%
filter(Season %in% c(2011),
tier==c(1)) %>%
mutate(Date = as.Date(Date)) -> england_2011
# Fit the default model, with the home team as team1.
gm_res <- goalmodel(goals1 = england_2011$hgoal, goals2 = england_2011$vgoal,
team1 = england_2011$home, team2=england_2011$visitor)
# Show the estimated attack and defense ratings and more.
summary(gm_res)
## Model sucsessfully fitted in 0.01 seconds
##
## Number of matches 380
## Number of teams 20
##
## Model Poisson
##
## Log Likelihood -1088.99
## AIC 2257.98
## R-squared 0.18
## Parameters (estimated) 40
## Parameters (fixed) 0
##
## Team Attack Defense
## Arsenal 0.36 0.03
## Aston Villa -0.33 -0.01
## Blackburn Rovers -0.04 -0.41
## Bolton Wanderers -0.09 -0.40
## Chelsea 0.23 0.10
## Everton -0.04 0.26
## Fulham -0.07 0.01
## Liverpool -0.10 0.26
## Manchester City 0.57 0.54
## Manchester United 0.53 0.41
## Newcastle United 0.08 0.01
## Norwich City 0.02 -0.25
## Queens Park Rangers -0.17 -0.24
## Stoke City -0.36 -0.01
## Sunderland -0.14 0.12
## Swansea City -0.16 0.02
## Tottenham Hotspur 0.24 0.21
## West Bromwich Albion -0.13 -0.00
## Wigan Athletic -0.19 -0.18
## Wolverhampton Wanderers -0.22 -0.45
## -------
## Intercept 0.13
## Home field advantage 0.27
The Dixon-Coles model
In a paper by Dixon and Coles (1997) they introduce a model that extend
the defaul model that modifies the probabilities for low scores, where
each team scores at most 1 goal each. The amount of adjustment is
determined by the parameter ρ (rho), which can be estimated from the
data. The Dixon-Coles model can be fitted by setting the dc argument to
TRUE.
# Fit the Dixon-Coles model.
gm_res_dc <- goalmodel(goals1 = england_2011$hgoal, goals2 = england_2011$vgoal,
team1 = england_2011$home, team2=england_2011$visitor,
dc=TRUE)
summary(gm_res_dc)
## Model sucsessfully fitted in 0.82 seconds
##
## Number of matches 380
## Number of teams 20
##
## Model Poisson
##
## Log Likelihood -1087.36
## AIC 2256.72
## R-squared 0.18
## Parameters (estimated) 41
## Parameters (fixed) 0
##
## Team Attack Defense
## Arsenal 0.37 0.03
## Aston Villa -0.31 -0.03
## Blackburn Rovers -0.06 -0.41
## Bolton Wanderers -0.08 -0.40
## Chelsea 0.23 0.09
## Everton -0.06 0.27
## Fulham -0.07 0.01
## Liverpool -0.11 0.25
## Manchester City 0.56 0.55
## Manchester United 0.52 0.43
## Newcastle United 0.10 0.00
## Norwich City 0.02 -0.25
## Queens Park Rangers -0.18 -0.23
## Stoke City -0.36 -0.01
## Sunderland -0.14 0.12
## Swansea City -0.15 0.01
## Tottenham Hotspur 0.24 0.21
## West Bromwich Albion -0.13 0.00
## Wigan Athletic -0.19 -0.17
## Wolverhampton Wanderers -0.21 -0.46
## -------
## Intercept 0.12
## Home field advantage 0.27
## Dixon-Coles adj. (rho) -0.13
Notice how the estimated ρ parameter is listed together with the other
parameters.
The Rue-Salvesen adjustment
In a paper by Rue and Salvesen (2001) they introduce an adjustment to
the goals scoring intesities λ1 and λ2:
log(λ1adj) = log(λ1) − γΔ1,2
log(λ2adj) = log(λ2) + γΔ1,2
where
Δ1,2 =
(attack1+defense1−attack2−defense2)/2
and γ is a parameter that modifies the amount of adjustment. This model
can be fitted by setting the rs argument to TRUE.
# Fit the model with the Rue-Salvesen adjustment.
gm_res_rs <- goalmodel(goals1 = england_2011$hgoal, goals2 = england_2011$vgoal,
team1 = england_2011$home, team2=england_2011$visitor,
rs=TRUE)
summary(gm_res_rs)
## Model sucsessfully fitted in 0.04 seconds
##
## Number of matches 380
## Number of teams 20
##
## Model Poisson
##
## Log Likelihood -1088.99
## AIC 2259.98
## R-squared 0.18
## Parameters (estimated) 41
## Parameters (fixed) 0
##
## Team Attack Defense
## Arsenal 0.36 0.03
## Aston Villa -0.33 -0.01
## Blackburn Rovers -0.04 -0.41
## Bolton Wanderers -0.09 -0.40
## Chelsea 0.23 0.10
## Everton -0.04 0.26
## Fulham -0.07 0.01
## Liverpool -0.10 0.26
## Manchester City 0.57 0.54
## Manchester United 0.53 0.41
## Newcastle United 0.08 0.01
## Norwich City 0.02 -0.25
## Queens Park Rangers -0.17 -0.24
## Stoke City -0.36 -0.01
## Sunderland -0.14 0.12
## Swansea City -0.16 0.02
## Tottenham Hotspur 0.24 0.21
## West Bromwich Albion -0.13 -0.00
## Wigan Athletic -0.19 -0.18
## Wolverhampton Wanderers -0.22 -0.45
## -------
## Intercept 0.13
## Home field advantage 0.27
## Rue-Salvesen adj. (gamma) 0.00
Making predictions
There is little use in fitting the model without making predictions.
Several functions are provided that make different types of predictions.
We can get the expected goals with the predict_expg function. Here the
return_df is set to TRUE, which returns a convenient data.frame.
to_predict1 <- c('Arsenal', 'Manchester United', 'Bolton Wanderers')
to_predict2 <- c('Fulham', 'Chelsea', 'Liverpool')
predict_expg(gm_res_dc, team1=to_predict1, team2=to_predict2, return_df = TRUE)
## team1 team2 expg1 expg2
## 1 Arsenal Fulham 2.117914 1.0257180
## 2 Manchester United Chelsea 2.291611 0.9287833
## 3 Bolton Wanderers Liverpool 1.069024 1.5088061
We can also get the probabilities of the final outcome (team1 win, draw,
team2 win) with the predict_result function.
predict_result(gm_res_dc, team1=to_predict1, team2=to_predict2, return_df = TRUE)
## team1 team2 p1 pd p2
## 1 Arsenal Fulham 0.6114897 0.2268230 0.1616873
## 2 Manchester United Chelsea 0.6685679 0.2051373 0.1262948
## 3 Bolton Wanderers Liverpool 0.2522617 0.2903403 0.4573981
Other functions are predict_goals which return matrices of the
probabilities of each scoreline, and predict_ou for getting the
probabilities for over/under total scores.
Weights
For predicton purposes it is a good idea to weight the influence of each
match on the estimates so that matches far back in time have less
influence than the more recent ones. The Dixon-Coles paper describe a
function that can help set the weights. This function is implemented in
the weights_dc function. The degree more recent matches are weighted
relative to earlier ones are determined by the parameter xi. Good values
of xi is usually somewhere between 0.001 and 0.003
Here we calculate the weights, plot them, and fit the default model with
the weights.
my_weights <- weights_dc(england_2011$Date, xi=0.0019)
length(my_weights)
## [1] 380
plot(england_2011$Date, my_weights)
gm_res_w <- goalmodel(goals1 = england_2011$hgoal, goals2 = england_2011$vgoal,
team1 = england_2011$home, team2=england_2011$visitor,
weights = my_weights)
Additional covariates
It is possible to use additional variables to predict the outcome. One
particular useful thing this can be used for is to have greater control
of the home field advantage. For example if you use data from several
leagues you might want to have a separate home field advantage for each
league, or be able to turn of the advantage for games played on neutral
ground.
Here is how to specify the default model with the home field advantage
specified as an additional covariate. Notice how the hfa argument is set
to FALSE, so that the default home field advantage factor is droped.
Also notice that only covariates for the first team is given, implicitly
setting all values of this covariate to 0 for team2.
# Create a n-by-1 matrix, containing only 1's, to indicate that all
# games have a home team.
xx1_hfa <- matrix(1, nrow=nrow(england_2011), ncol=1)
colnames(xx1_hfa) <- 'hfa' # must be named.
# Fit the default model, with the home team as team1.
gm_res_hfa <- goalmodel(goals1 = england_2011$hgoal, goals2 = england_2011$vgoal,
team1 = england_2011$home, team2=england_2011$visitor,
hfa = FALSE, # Turn of the built-in Home Field Advantage
x1 = xx1_hfa) # Provide the covariate matrix.
# the coefficients for the additional covariates are found in the 'beta' element in
# the parameter list. They will also show up in the summary.
gm_res_hfa$parameters$beta
## hfa
## 0.2680093
Alternatives to the Poisson model
By the default the Poisson distribution is used to model the goals. The
Poisson distribution has a lot of nice properties, but also some
limitations. For instance there is only one parameter that describes
both the expected value (the average) and the variance. It is also
limited to integer number of goals.
The Negative Binomial model
The Negative Binomial (NB) model is a more flexible model that allows
higher variance than the Poisson model. The NB model introduces another
parameter called the “dispersion” parameter. This parameter has to be
positive, and the greater the value of this parameter, the more variance
there is relative to the Poisson model with the same λ. When the
variance is greater than the expectation, it is referred to as
overdispersion. When the parameter is close to zero, it the result is
approximately the Poisson. The NB model can be fitted by setting the
“model” argument to ‘negbin’. In this example the dispersion parameter
is close to zero.
# Fit the model with overdispersion.
gm_res_nb <- goalmodel(goals1 = england_2011$hgoal, goals2 = england_2011$vgoal,
team1 = england_2011$home, team2=england_2011$visitor,
model='negbin')
gm_res_nb$parameters$dispersion
## [1] 4.539982e-05
The Conway-Maxwell-Poisson model
While the Negative Binomial model is more flexible than the Poisson
model in that it can model excess variance, the Conway-Maxwell-Poisson
(CMP) model is even more flexible. The CMP can be both over-dispersed
and under-dispersed. The dispersion parameter controlling the amount of
dispersion is called υ (upsilon), and behaves a bit different than the
dispersion parameter in the Negative Binomial. When υ = 1 the model
becomes the Poisson model. If υ is greater than 1 the variance is less
than the expectation, and we have under-dispersion. If υ is between 0
and 1, there is over-dipsersion.
Unfortunately, the CMP distribution is a bit difficult to work with for
several reasons. It is not built-in in R, so I have implemented some of
the relevant distribution functions (called dCMP and pCMP). It also has
a normalizing constant that is a bit difficult to compute, and must be
approximated. The distribution functions have an argument called ‘error’
that determines the amount of approximation, and has a default that is
precise enough for our needs.
The perhaps most confusing thing with the CMP model is that it has a
parameter called λ, and this parameter does NOT behave like the λ in
the Poisson or NB models. In those models the lambda is the same as the
expectation, but this is not the case for the CMP. To make matters
worse, there is no simple forumla for finding the expectation for a
given set of values of λ and υ, although an approximate formula
exists. Therefore I have included two functions to compute the
expectation (eCMP) and another to find the λ for a given set of values
for the expectation and dispersion (lambdaCMP).
To fit and use a model with the CMP distribution is easy, you can just
set model == ‘cmp’ in the goalmodel function. However, because of all
the complications with this model, estimating the parameters can be
really slow. Under the hood, I have tried to speed things up by using
approximations where I could, but it is still slow. I therefore recomend
using two-step estimation to fit this model. If you fit this model, and
see that there is overdispersion (dispersion \< 1), it can be a good
idea to use the NB model instead.
# Fit the CMP model, with the attack and defence parameters fixed to the values in the default model.
gm_res_cmp_2s <- goalmodel(goals1 = england_2011$hgoal, goals2 = england_2011$vgoal,
team1 = england_2011$home, team2=england_2011$visitor,
model='cmp', fixed_params = gm_res$parameters)
# Take a look at the dispersion. It indicates under-dispersion.
gm_res_cmp_2s$parameters$dispersion
## [1] 1.122462
The Gaussian model
It is also possible to use model = ‘gaussian’ to fit models using the
Gaussian (or normal) distribution. This option is intended to be used
when the scores are decimal numbers. This could be useful if instead of
actual goals scored, you wanted to fit a model to expected goals. You
can’t combine this with the Dixon-Coles adjustment, except perhaps using
a two step procedure. If you use a Gaussian model to make predictions,
the Poisson distribution will be used, and a warning will be given.
Fixed parameters - Two step estimation.
If you want to, you can set some of the parameters to a fixed value that
is kept constant during model fitting. In other words, you don’t
estimate the fixed parameters from the data, but set them manually
beforehand and all other parameters are estimated with this in mind.
This feature can be useful for doing two-step estimation. Sometimes the
model fitting can be unstable, especially when the model is complicated.
It can then be a good idea to first fit a simpler model, like the
default model, and then with the parameters from the simpler model kept
constant, estimate the rest of the parameters in a second step.
Here is an example of two-step fitting of the Dixon-Coles model.
# Fit the Dixon-Coles model, with most of the parameters fixed to the values in the default model.
gm_res_dc_2s <- goalmodel(goals1 = england_2011$hgoal, goals2 = england_2011$vgoal,
team1 = england_2011$home, team2=england_2011$visitor,
dc=TRUE, fixed_params = gm_res$parameters)
# Compare the two estimates of rho. They are pretty similar.
gm_res_dc_2s$parameters$rho
## [1] -0.1261445
gm_res_dc$parameters$rho
## [1] -0.1335741
Offset - Varying playing time
Sometimes the playing time can be extended to extra time. In football
this usually happens when there is a draw in a knock-out tournament,
where an additional 30 minutes is added. To handle this, we must add
what is called an offset to the model. You can read more about the
details on
wikipedia.
There is no offset option in the goalmodel package, but we can add it as
an extra covariate, and fix the parameter for that covariate to be 1. In
the example below we add a game from the 2011-12 FA Cup to the data set
which was extended to extra time. Then we create the offset variable and
fit the model with the fixed parameter. Note that the offset variable
must be log-transformed. In this example, the extra game that is added
is the only game in the data set where Middlesbrough is playing, so the
estimates will be a bit unstable. Also note that the number of minutes
played is divided by 90, the ordinary playing time, otherwise the model
fitting tend to get unstable.
# Get matches from the FA cup that was extended to extra time.
facup %>%
filter(Season == 2011,
round >= 4,
aet == 'yes') %>%
select(Date, home, visitor, hgoal, vgoal, aet) %>%
mutate(Date = as.Date(Date)) -> facup_2011
facup_2011
## Date home visitor hgoal vgoal aet
## 16005 2012-02-08 Middlesbrough Sunderland 1 2 yes
# Merge the additional game with the origianl data frame.
england_2011 %>%
bind_rows(facup_2011) %>%
mutate(aet = replace(aet, is.na(aet), 'no')) -> england_2011_2
# Create empty matrix for additional covariates
xx_offset <- matrix(nrow=nrow(england_2011_2))
colnames(xx_offset) <- 'Offset'
# Add data
xx_offset[,'Offset'] <- 90/90 # Per 90 minutes
# Extra time is 30 minutes added.
xx_offset[england_2011_2$aet == 'yes','Offset'] <- (90+30)/90
# Offset must be log-transformed.
xx_offset <- log(xx_offset)
# Take a look.
tail(xx_offset)
## Offset
## [376,] 0.0000000
## [377,] 0.0000000
## [378,] 0.0000000
## [379,] 0.0000000
## [380,] 0.0000000
## [381,] 0.2876821
# fit the model
gm_res_offset <- goalmodel(goals1 = england_2011_2$hgoal, goals2 = england_2011_2$vgoal,
team1 = england_2011_2$home, team2=england_2011_2$visitor,
# The offset must be added to both x1 and x2.
x1 = xx_offset, x2 = xx_offset,
# The offset parameter must be fixed to 1.
fixed_params = list(beta = c('Offset' = 1)))
summary(gm_res_offset)
## Model sucsessfully fitted in 0.54 seconds
##
## Number of matches 381
## Number of teams 21
##
## Model Poisson
##
## Log Likelihood -1091.30
## AIC 2266.60
## R-squared 0.18
## Parameters (estimated) 42
## Parameters (fixed) 1
##
## Team Attack Defense
## Arsenal 0.39 0.05
## Aston Villa -0.30 0.01
## Blackburn Rovers -0.02 -0.39
## Bolton Wanderers -0.06 -0.38
## Chelsea 0.25 0.12
## Everton -0.01 0.27
## Fulham -0.04 0.03
## Liverpool -0.08 0.28
## Manchester City 0.60 0.56
## Manchester United 0.56 0.43
## Middlesbrough -0.52 -0.39
## Newcastle United 0.11 0.03
## Norwich City 0.05 -0.23
## Queens Park Rangers -0.14 -0.22
## Stoke City -0.33 0.01
## Sunderland -0.11 0.14
## Swansea City -0.13 0.04
## Tottenham Hotspur 0.26 0.23
## West Bromwich Albion -0.11 0.02
## Wigan Athletic -0.17 -0.16
## Wolverhampton Wanderers -0.20 -0.43
## -------
## Intercept 0.12
## Home field advantage 0.27
## Offset 1.00
Modifying parameters
In addition to fixing parameters before the model is fit, you can also
manually set the parameters after the model is fitted. This can be
useful if you believe the attack or defence parameters given by the
model is inacurate because you have some additional information that is
not captured in the data.
Modifying the parameters can also be useful for the same reasons as you
might want to fit the model in two (or more) steps. For intance, if you
look at the log-likelihood of the default model, and the model with the
Rue-Salvesen adjustment, you will notice that they are almost identical.
This is probably due to the model being nearly unidentifiable. That does
not neccecarily mean that the Rue-Salvesen adjustment does not improve
prediction, it might mean that it is hard to estimate the amount of
adjustment based on the available data. In the paper by Rue & Salvesen
they find that setting γ to 0.1 is a good choice. Here is how this can
be done:
# Copy the default model.
gm_res_rs2 <- gm_res
# set the gamma parameter to 0.1
gm_res_rs2$parameters$gamma <- 0.1
# Make a few predictions to compare.
predict_result(gm_res, team1=to_predict1, team2=to_predict2, return_df = TRUE)
## team1 team2 p1 pd p2
## 1 Arsenal Fulham 0.6197118 0.2035688 0.1767195
## 2 Manchester United Chelsea 0.6736194 0.1843756 0.1420050
## 3 Bolton Wanderers Liverpool 0.2612112 0.2567706 0.4820183
predict_result(gm_res_rs2, team1=to_predict1, team2=to_predict2, return_df = TRUE)
## team1 team2 p1 pd p2
## 1 Arsenal Fulham 0.6047391 0.2083325 0.1869284
## 2 Manchester United Chelsea 0.6539250 0.1917532 0.1543218
## 3 Bolton Wanderers Liverpool 0.2780789 0.2603511 0.4615700
Evaluating predictions
To evaluate and compare the quality of predictions between different
models, you need a measure of how well the predictions match the
observed outcomes. The score_predictions function can help with this,
and implements three different scoring rules. Here is an example of how
it works. The example is done in-sample, but a proper evaluation should
compute scores out-of-sample.
my_predictions_default <- predict_result(gm_res, team1 = to_predict1, team2 = to_predict2, return_df = TRUE)
my_predictions_dc <- predict_result(gm_res_dc, team1 = to_predict1, team2 = to_predict2, return_df = TRUE)
england_2011 %>%
# Select the games to evaluate
mutate(teams_string = paste(home, '-', visitor)) %>%
filter(teams_string %in% paste(to_predict1, '-', to_predict2)) %>%
# Transform the result column (H, D, A) to reflect the column names of the prediction matrix.
mutate(result2 = c('H' = 'p1', 'D'='pd', 'A' = 'p2')[result]) %>%
select(home, visitor, result2) -> england_2011_to_score
predictions_scores_default <- score_predictions(predictions = my_predictions_default[,c('p1','pd','p2')],
observed = england_2011_to_score$result2,
score = c('log', 'brier', 'rps'))
predictions_scores_dc <- score_predictions(predictions = my_predictions_dc[,c('p1','pd','p2')],
observed = england_2011_to_score$result2,
score = c('log', 'brier', 'rps'))
# Total log score for default model, 4.27
sum(predictions_scores_default$log)
## [1] 4.273418
# Total log score for Dixon-Coles model, 4.33, which is better than for default.
sum(predictions_scores_dc$log)
## [1] 4.334922
Miscalaneous
Reverse engineering expected goals
The function predict_result that’s demonstrated above uses the
underlying statistical model and the expected goals from fitted
goalmodel to compute the probabilities for win-draw-lose. The function
expg_from_probabilities can be used to reverse this procedure, going
from probabilities to the underlying expected goals.
This can be used to reverse-engineer predictions from other models or
even bookmaker odds. The expected goals extracted from this method can
for example be used as input to the goalmodel to get attack and defense
ratings, and to make new predictions.
Here is some code demonstrating that it can recover the expected goals
which are made from a fitted goalmodel with the Poisson distribution.
# Win-Draw-Lose probabilities.
wdl_probs <- predict_result(gm_res, team1=to_predict1, team2=to_predict2,
return_df = TRUE)
wdl_probs
## team1 team2 p1 pd p2
## 1 Arsenal Fulham 0.6197118 0.2035688 0.1767195
## 2 Manchester United Chelsea 0.6736194 0.1843756 0.1420050
## 3 Bolton Wanderers Liverpool 0.2612112 0.2567706 0.4820183
expg_reveng <- expg_from_probabilities(wdl_probs[,c('p1', 'pd', 'p2')])
expg_reveng$expg
## [,1] [,2]
## [1,] 2.098592 1.0272203
## [2,] 2.278604 0.9458234
## [3,] 1.049001 1.5224753
# Compare with expected goals predictions from above.
predict_expg(gm_res_dc, team1=to_predict1, team2=to_predict2, return_df = TRUE)
## team1 team2 expg1 expg2
## 1 Arsenal Fulham 2.117914 1.0257180
## 2 Manchester United Chelsea 2.291611 0.9287833
## 3 Bolton Wanderers Liverpool 1.069024 1.5088061
You can also specify the Dixon-Coles parameter ρ (rho). Unfortunately it
is not possible to extract both the two expected goals and rho, since
there will typically be a large number of combinations of the three
parameters that yield the same probabilities. With rho given as a
constant it is however possible to extract the expected goals. Here is
how you can recover the expg from the DC model fitted above, using rho =
-0.13.
# Win-Draw-Lose probabilities from DC model
wdl_probs_dc <- predict_result(gm_res_dc, team1=to_predict1, team2=to_predict2,
return_df = TRUE)
wdl_probs_dc
## team1 team2 p1 pd p2
## 1 Arsenal Fulham 0.6114897 0.2268230 0.1616873
## 2 Manchester United Chelsea 0.6685679 0.2051373 0.1262948
## 3 Bolton Wanderers Liverpool 0.2522617 0.2903403 0.4573981
# Use the rho argument to specify rho.
expg_reveng_dc <- expg_from_probabilities(wdl_probs_dc[,c('p1', 'pd', 'p2')],
rho = -0.13)
expg_reveng_dc$expg
## [,1] [,2]
## [1,] 2.109461 1.0195470
## [2,] 2.282623 0.9229529
## [3,] 1.062882 1.5016514
# Compare with expected goals predictions from above.
predict_expg(gm_res_dc, team1=to_predict1, team2=to_predict2, return_df = TRUE)
## team1 team2 expg1 expg2
## 1 Arsenal Fulham 2.117914 1.0257180
## 2 Manchester United Chelsea 2.291611 0.9287833
## 3 Bolton Wanderers Liverpool 1.069024 1.5088061
Other packages
There are some other packages out there that have similar functionality
as this package. The regista
package implements a more general interface to fit the Dixon-Coles
model. The
fbRanks
package implements only the Dixon-Coles time weighting functionality,
but has a lot of other features. In
footBayes you can fit
models with Bayesian procedures. The package includes several models,
including bivariate Poisson models and models with time-varying attack
and defense parameters.
References
- Mark J. Dixon and Stuart G. Coles (1997) Modelling Association
Football Scores and Inefficiencies in the Football Betting Market.
- Håvard Rue and Øyvind Salvesen (2001) Prediction and Retrospective
Analysis of Soccer Matches in a League.
opisthokonta/goalmodel documentation built on April 3, 2024, 1:32 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
The goalmodel package let you build prediction models for the number of goals scored in sport games. The models are primarily aimed at modelling and predicting football (soccer) scores, but could also be applicable for similar sports, such as hockey and handball.
If you find this package useful, please consider supporting the development at https://ko-fi.com/opisthokonta
Installation
install.packages("devtools")
devtools::install_github("opisthokonta/goalmodel")
Whats new
Version 0.6.4 (development version)
- New function rho.ml to estimate the rho parameter of the Dixon-Coles model.
- Fixed important bug that caused the defense parameters to be slightly wrong in certain models. This affected mostly the defense parameter for the first team in alphabetical order. In terms of prediction error this affected the expected goals to be off by less than 0.1 goals for games involving the first team. Models that were affected were the default Poisson model and Gaussian models.
- Can now fit a hurdle model (See A. Owen (2017) The Application of Hurdle Models to Accurately Model 0-0 Draws in Predictive Models of Football Match Outcomes), but prediction methods are not yet implemented.
- Small bugfixes in the
score_predictionsfunction. - README now has an example of how to use the
score_predictionsfunction.
Version 0.6
- New functions predict_btts() and pbtts() for computing Both-Teams-To-Score probabilities.
- New function score_predictions(), for evaluating prediction accuracy using Ranked Probability Scores, Brier scores, and Log scores.
- New function league_table() that computes a league table with points and useful summary statistics.
Version 0.5
- New function expg_from_ou() - Finds the expected goals from over/under probabilities.
- New function p1x2() - Compute win and draw probabilities from expected goals.
- Bugfix: predict_goals() with return_df=TRUE did not work properly when the team arguments were factors.
See NEWS.md for changes for complete version history.
Overview
The main function in this package is the goalmodel function, which is
used to fit a range of different models, including using different
probability distributions and predictors. It also allows for setting
parameters to a fixed value, which can be useful for more advanced
techniques such as two-step estimation.
The following probability distributions are available for fitting with
the goalmodel function via the model argument:
-
poisson, for the independent Poisson model, sometimes also called the double Poisson model. -
negbin, for the (independent) negative Binomial model, which allows for more variance (overdisperison) than the Poisson model. -
cmp, for the (independent) Conway-Maxwell-Poisson, which allows for both underdispersion and overdispersion. -
gaussianfor the (independent) bivariate Gaussian (also known as the normal distribution) model. This is useful for modelling scores that might be decimal numbers.
In addition, the Dixon-Coles model (Dixon and Coles, 1997) is also
available. The model extend the Poisson model so that it modifies the
probabilities for low scores. The Dixon-Coles model can be fitted by
setting the dc argument to TRUE. Hence, the Dixon-Coles adjustment
can be used not only with Poisson models, but also with the Negative
Binomial and Conway-Maxwell-Poisson models.
With a fitted model, you would want to make predictions. This package has functions for making different sorts of predictions:
-
predict_expg, for predicting expected goals. -
predict_goals, for the scoreline probabilities. -
predict_result, for the final result (home-draw-away or 1x2) probabilities. -
predict_ou, for over/under probabilities. -
predict_btts, for both teams to score probabilities.
There are also functions that can be useful for making predictions without relying on a fitted model, but instead uses expected goals directly. You can use all of the probability distributions listed above (except ‘gaussian’) as the underlying model, including the Dixon-Coles adjustment.
-
p1x2 -
pbtts
There are a couple of functions for ‘reverse engineering’ results from other models. These functions relies on the Poisson model.
-
expg_from_ouconverts over/under probabilities to expected goals. -
expg_from_probabilitiesconverts 1x2 probabilities to expected goals.
Misc
-
weights_dc, make weights to let results further back in time have less influence on the estimates. -
score_predictionscomputes ‘scoring rules’ to assess the quality of predictions. -
days_since_last_matchcalculates how long it was since the previous game for each team. -
matches_last_xdayscalculates how many games each team have played for a given period in the past.
There are also some functions to work directly with the Conway-Maxwell-Poisson probability distribution.
-
dCMP,pCMP, andqCMP, are the probability, cumulative and quantile functions for the CMP distribution. -
eCMPandlambdaCMPconverts between expected value and the lambda parameter in the CMP distribution. -
upsilon.mlestimates the dispersion parameter (upsilon) from data and some given rate (lambda) or mean (mu) parameters.
The default model
The goalmodel package models the goal scoring intensity so that it is a function of the two opposing teams attack and defense ratings. Let λ1 and λ2) be the goal scoring intensities for the two sides. The default model in the goalmodel package then models these as
log(λ1) = μ + hfa + attack1 − defense2
log(λ2) = μ + attack2 − defense1
where μ is the intercept (approximately the average number of goals scored) and hfa is the home field advantage given to team 1. By default the number of goals scored, Y1 and Y2 are distributed as
Y1 ∼ Poisson(λ1)
Y2 ∼ Poisson(λ2)
The default model can be modified in numerous ways as detailed in this vignette.
Fitting the default model
Models are fitted with the goalmodel function. The minimum required input are vectors containing the number of goals scored and the names of the teams. Here I use data from the excellent engsoccerdata package.
library(goalmodel)
library(dplyr) # Useful for data manipulation.
library(engsoccerdata) # Some data.
library(Rcpp)
# Load data from English Premier League, 2011-12 season.
england %>%
filter(Season %in% c(2011),
tier==c(1)) %>%
mutate(Date = as.Date(Date)) -> england_2011
# Fit the default model, with the home team as team1.
gm_res <- goalmodel(goals1 = england_2011$hgoal, goals2 = england_2011$vgoal,
team1 = england_2011$home, team2=england_2011$visitor)
# Show the estimated attack and defense ratings and more.
summary(gm_res)
## Model sucsessfully fitted in 0.01 seconds
##
## Number of matches 380
## Number of teams 20
##
## Model Poisson
##
## Log Likelihood -1088.99
## AIC 2257.98
## R-squared 0.18
## Parameters (estimated) 40
## Parameters (fixed) 0
##
## Team Attack Defense
## Arsenal 0.36 0.03
## Aston Villa -0.33 -0.01
## Blackburn Rovers -0.04 -0.41
## Bolton Wanderers -0.09 -0.40
## Chelsea 0.23 0.10
## Everton -0.04 0.26
## Fulham -0.07 0.01
## Liverpool -0.10 0.26
## Manchester City 0.57 0.54
## Manchester United 0.53 0.41
## Newcastle United 0.08 0.01
## Norwich City 0.02 -0.25
## Queens Park Rangers -0.17 -0.24
## Stoke City -0.36 -0.01
## Sunderland -0.14 0.12
## Swansea City -0.16 0.02
## Tottenham Hotspur 0.24 0.21
## West Bromwich Albion -0.13 -0.00
## Wigan Athletic -0.19 -0.18
## Wolverhampton Wanderers -0.22 -0.45
## -------
## Intercept 0.13
## Home field advantage 0.27
The Dixon-Coles model
In a paper by Dixon and Coles (1997) they introduce a model that extend the defaul model that modifies the probabilities for low scores, where each team scores at most 1 goal each. The amount of adjustment is determined by the parameter ρ (rho), which can be estimated from the data. The Dixon-Coles model can be fitted by setting the dc argument to TRUE.
# Fit the Dixon-Coles model.
gm_res_dc <- goalmodel(goals1 = england_2011$hgoal, goals2 = england_2011$vgoal,
team1 = england_2011$home, team2=england_2011$visitor,
dc=TRUE)
summary(gm_res_dc)
## Model sucsessfully fitted in 0.82 seconds
##
## Number of matches 380
## Number of teams 20
##
## Model Poisson
##
## Log Likelihood -1087.36
## AIC 2256.72
## R-squared 0.18
## Parameters (estimated) 41
## Parameters (fixed) 0
##
## Team Attack Defense
## Arsenal 0.37 0.03
## Aston Villa -0.31 -0.03
## Blackburn Rovers -0.06 -0.41
## Bolton Wanderers -0.08 -0.40
## Chelsea 0.23 0.09
## Everton -0.06 0.27
## Fulham -0.07 0.01
## Liverpool -0.11 0.25
## Manchester City 0.56 0.55
## Manchester United 0.52 0.43
## Newcastle United 0.10 0.00
## Norwich City 0.02 -0.25
## Queens Park Rangers -0.18 -0.23
## Stoke City -0.36 -0.01
## Sunderland -0.14 0.12
## Swansea City -0.15 0.01
## Tottenham Hotspur 0.24 0.21
## West Bromwich Albion -0.13 0.00
## Wigan Athletic -0.19 -0.17
## Wolverhampton Wanderers -0.21 -0.46
## -------
## Intercept 0.12
## Home field advantage 0.27
## Dixon-Coles adj. (rho) -0.13
Notice how the estimated ρ parameter is listed together with the other parameters.
The Rue-Salvesen adjustment
In a paper by Rue and Salvesen (2001) they introduce an adjustment to the goals scoring intesities λ1 and λ2:
log(λ1adj) = log(λ1) − γΔ1,2
log(λ2adj) = log(λ2) + γΔ1,2
where
Δ1,2 = (attack1+defense1−attack2−defense2)/2
and γ is a parameter that modifies the amount of adjustment. This model can be fitted by setting the rs argument to TRUE.
# Fit the model with the Rue-Salvesen adjustment.
gm_res_rs <- goalmodel(goals1 = england_2011$hgoal, goals2 = england_2011$vgoal,
team1 = england_2011$home, team2=england_2011$visitor,
rs=TRUE)
summary(gm_res_rs)
## Model sucsessfully fitted in 0.04 seconds
##
## Number of matches 380
## Number of teams 20
##
## Model Poisson
##
## Log Likelihood -1088.99
## AIC 2259.98
## R-squared 0.18
## Parameters (estimated) 41
## Parameters (fixed) 0
##
## Team Attack Defense
## Arsenal 0.36 0.03
## Aston Villa -0.33 -0.01
## Blackburn Rovers -0.04 -0.41
## Bolton Wanderers -0.09 -0.40
## Chelsea 0.23 0.10
## Everton -0.04 0.26
## Fulham -0.07 0.01
## Liverpool -0.10 0.26
## Manchester City 0.57 0.54
## Manchester United 0.53 0.41
## Newcastle United 0.08 0.01
## Norwich City 0.02 -0.25
## Queens Park Rangers -0.17 -0.24
## Stoke City -0.36 -0.01
## Sunderland -0.14 0.12
## Swansea City -0.16 0.02
## Tottenham Hotspur 0.24 0.21
## West Bromwich Albion -0.13 -0.00
## Wigan Athletic -0.19 -0.18
## Wolverhampton Wanderers -0.22 -0.45
## -------
## Intercept 0.13
## Home field advantage 0.27
## Rue-Salvesen adj. (gamma) 0.00
Making predictions
There is little use in fitting the model without making predictions. Several functions are provided that make different types of predictions.
We can get the expected goals with the predict_expg function. Here the return_df is set to TRUE, which returns a convenient data.frame.
to_predict1 <- c('Arsenal', 'Manchester United', 'Bolton Wanderers')
to_predict2 <- c('Fulham', 'Chelsea', 'Liverpool')
predict_expg(gm_res_dc, team1=to_predict1, team2=to_predict2, return_df = TRUE)
## team1 team2 expg1 expg2
## 1 Arsenal Fulham 2.117914 1.0257180
## 2 Manchester United Chelsea 2.291611 0.9287833
## 3 Bolton Wanderers Liverpool 1.069024 1.5088061
We can also get the probabilities of the final outcome (team1 win, draw, team2 win) with the predict_result function.
predict_result(gm_res_dc, team1=to_predict1, team2=to_predict2, return_df = TRUE)
## team1 team2 p1 pd p2
## 1 Arsenal Fulham 0.6114897 0.2268230 0.1616873
## 2 Manchester United Chelsea 0.6685679 0.2051373 0.1262948
## 3 Bolton Wanderers Liverpool 0.2522617 0.2903403 0.4573981
Other functions are predict_goals which return matrices of the probabilities of each scoreline, and predict_ou for getting the probabilities for over/under total scores.
Weights
For predicton purposes it is a good idea to weight the influence of each match on the estimates so that matches far back in time have less influence than the more recent ones. The Dixon-Coles paper describe a function that can help set the weights. This function is implemented in the weights_dc function. The degree more recent matches are weighted relative to earlier ones are determined by the parameter xi. Good values of xi is usually somewhere between 0.001 and 0.003
Here we calculate the weights, plot them, and fit the default model with the weights.
my_weights <- weights_dc(england_2011$Date, xi=0.0019)
length(my_weights)
## [1] 380
plot(england_2011$Date, my_weights)
gm_res_w <- goalmodel(goals1 = england_2011$hgoal, goals2 = england_2011$vgoal,
team1 = england_2011$home, team2=england_2011$visitor,
weights = my_weights)
Additional covariates
It is possible to use additional variables to predict the outcome. One particular useful thing this can be used for is to have greater control of the home field advantage. For example if you use data from several leagues you might want to have a separate home field advantage for each league, or be able to turn of the advantage for games played on neutral ground.
Here is how to specify the default model with the home field advantage specified as an additional covariate. Notice how the hfa argument is set to FALSE, so that the default home field advantage factor is droped. Also notice that only covariates for the first team is given, implicitly setting all values of this covariate to 0 for team2.
# Create a n-by-1 matrix, containing only 1's, to indicate that all
# games have a home team.
xx1_hfa <- matrix(1, nrow=nrow(england_2011), ncol=1)
colnames(xx1_hfa) <- 'hfa' # must be named.
# Fit the default model, with the home team as team1.
gm_res_hfa <- goalmodel(goals1 = england_2011$hgoal, goals2 = england_2011$vgoal,
team1 = england_2011$home, team2=england_2011$visitor,
hfa = FALSE, # Turn of the built-in Home Field Advantage
x1 = xx1_hfa) # Provide the covariate matrix.
# the coefficients for the additional covariates are found in the 'beta' element in
# the parameter list. They will also show up in the summary.
gm_res_hfa$parameters$beta
## hfa
## 0.2680093
Alternatives to the Poisson model
By the default the Poisson distribution is used to model the goals. The Poisson distribution has a lot of nice properties, but also some limitations. For instance there is only one parameter that describes both the expected value (the average) and the variance. It is also limited to integer number of goals.
The Negative Binomial model
The Negative Binomial (NB) model is a more flexible model that allows higher variance than the Poisson model. The NB model introduces another parameter called the “dispersion” parameter. This parameter has to be positive, and the greater the value of this parameter, the more variance there is relative to the Poisson model with the same λ. When the variance is greater than the expectation, it is referred to as overdispersion. When the parameter is close to zero, it the result is approximately the Poisson. The NB model can be fitted by setting the “model” argument to ‘negbin’. In this example the dispersion parameter is close to zero.
# Fit the model with overdispersion.
gm_res_nb <- goalmodel(goals1 = england_2011$hgoal, goals2 = england_2011$vgoal,
team1 = england_2011$home, team2=england_2011$visitor,
model='negbin')
gm_res_nb$parameters$dispersion
## [1] 4.539982e-05
The Conway-Maxwell-Poisson model
While the Negative Binomial model is more flexible than the Poisson model in that it can model excess variance, the Conway-Maxwell-Poisson (CMP) model is even more flexible. The CMP can be both over-dispersed and under-dispersed. The dispersion parameter controlling the amount of dispersion is called υ (upsilon), and behaves a bit different than the dispersion parameter in the Negative Binomial. When υ = 1 the model becomes the Poisson model. If υ is greater than 1 the variance is less than the expectation, and we have under-dispersion. If υ is between 0 and 1, there is over-dipsersion.
Unfortunately, the CMP distribution is a bit difficult to work with for several reasons. It is not built-in in R, so I have implemented some of the relevant distribution functions (called dCMP and pCMP). It also has a normalizing constant that is a bit difficult to compute, and must be approximated. The distribution functions have an argument called ‘error’ that determines the amount of approximation, and has a default that is precise enough for our needs.
The perhaps most confusing thing with the CMP model is that it has a parameter called λ, and this parameter does NOT behave like the λ in the Poisson or NB models. In those models the lambda is the same as the expectation, but this is not the case for the CMP. To make matters worse, there is no simple forumla for finding the expectation for a given set of values of λ and υ, although an approximate formula exists. Therefore I have included two functions to compute the expectation (eCMP) and another to find the λ for a given set of values for the expectation and dispersion (lambdaCMP).
To fit and use a model with the CMP distribution is easy, you can just set model == ‘cmp’ in the goalmodel function. However, because of all the complications with this model, estimating the parameters can be really slow. Under the hood, I have tried to speed things up by using approximations where I could, but it is still slow. I therefore recomend using two-step estimation to fit this model. If you fit this model, and see that there is overdispersion (dispersion \< 1), it can be a good idea to use the NB model instead.
# Fit the CMP model, with the attack and defence parameters fixed to the values in the default model.
gm_res_cmp_2s <- goalmodel(goals1 = england_2011$hgoal, goals2 = england_2011$vgoal,
team1 = england_2011$home, team2=england_2011$visitor,
model='cmp', fixed_params = gm_res$parameters)
# Take a look at the dispersion. It indicates under-dispersion.
gm_res_cmp_2s$parameters$dispersion
## [1] 1.122462
The Gaussian model
It is also possible to use model = ‘gaussian’ to fit models using the Gaussian (or normal) distribution. This option is intended to be used when the scores are decimal numbers. This could be useful if instead of actual goals scored, you wanted to fit a model to expected goals. You can’t combine this with the Dixon-Coles adjustment, except perhaps using a two step procedure. If you use a Gaussian model to make predictions, the Poisson distribution will be used, and a warning will be given.
Fixed parameters - Two step estimation.
If you want to, you can set some of the parameters to a fixed value that is kept constant during model fitting. In other words, you don’t estimate the fixed parameters from the data, but set them manually beforehand and all other parameters are estimated with this in mind.
This feature can be useful for doing two-step estimation. Sometimes the model fitting can be unstable, especially when the model is complicated. It can then be a good idea to first fit a simpler model, like the default model, and then with the parameters from the simpler model kept constant, estimate the rest of the parameters in a second step.
Here is an example of two-step fitting of the Dixon-Coles model.
# Fit the Dixon-Coles model, with most of the parameters fixed to the values in the default model.
gm_res_dc_2s <- goalmodel(goals1 = england_2011$hgoal, goals2 = england_2011$vgoal,
team1 = england_2011$home, team2=england_2011$visitor,
dc=TRUE, fixed_params = gm_res$parameters)
# Compare the two estimates of rho. They are pretty similar.
gm_res_dc_2s$parameters$rho
## [1] -0.1261445
gm_res_dc$parameters$rho
## [1] -0.1335741
Offset - Varying playing time
Sometimes the playing time can be extended to extra time. In football this usually happens when there is a draw in a knock-out tournament, where an additional 30 minutes is added. To handle this, we must add what is called an offset to the model. You can read more about the details on wikipedia.
There is no offset option in the goalmodel package, but we can add it as an extra covariate, and fix the parameter for that covariate to be 1. In the example below we add a game from the 2011-12 FA Cup to the data set which was extended to extra time. Then we create the offset variable and fit the model with the fixed parameter. Note that the offset variable must be log-transformed. In this example, the extra game that is added is the only game in the data set where Middlesbrough is playing, so the estimates will be a bit unstable. Also note that the number of minutes played is divided by 90, the ordinary playing time, otherwise the model fitting tend to get unstable.
# Get matches from the FA cup that was extended to extra time.
facup %>%
filter(Season == 2011,
round >= 4,
aet == 'yes') %>%
select(Date, home, visitor, hgoal, vgoal, aet) %>%
mutate(Date = as.Date(Date)) -> facup_2011
facup_2011
## Date home visitor hgoal vgoal aet
## 16005 2012-02-08 Middlesbrough Sunderland 1 2 yes
# Merge the additional game with the origianl data frame.
england_2011 %>%
bind_rows(facup_2011) %>%
mutate(aet = replace(aet, is.na(aet), 'no')) -> england_2011_2
# Create empty matrix for additional covariates
xx_offset <- matrix(nrow=nrow(england_2011_2))
colnames(xx_offset) <- 'Offset'
# Add data
xx_offset[,'Offset'] <- 90/90 # Per 90 minutes
# Extra time is 30 minutes added.
xx_offset[england_2011_2$aet == 'yes','Offset'] <- (90+30)/90
# Offset must be log-transformed.
xx_offset <- log(xx_offset)
# Take a look.
tail(xx_offset)
## Offset
## [376,] 0.0000000
## [377,] 0.0000000
## [378,] 0.0000000
## [379,] 0.0000000
## [380,] 0.0000000
## [381,] 0.2876821
# fit the model
gm_res_offset <- goalmodel(goals1 = england_2011_2$hgoal, goals2 = england_2011_2$vgoal,
team1 = england_2011_2$home, team2=england_2011_2$visitor,
# The offset must be added to both x1 and x2.
x1 = xx_offset, x2 = xx_offset,
# The offset parameter must be fixed to 1.
fixed_params = list(beta = c('Offset' = 1)))
summary(gm_res_offset)
## Model sucsessfully fitted in 0.54 seconds
##
## Number of matches 381
## Number of teams 21
##
## Model Poisson
##
## Log Likelihood -1091.30
## AIC 2266.60
## R-squared 0.18
## Parameters (estimated) 42
## Parameters (fixed) 1
##
## Team Attack Defense
## Arsenal 0.39 0.05
## Aston Villa -0.30 0.01
## Blackburn Rovers -0.02 -0.39
## Bolton Wanderers -0.06 -0.38
## Chelsea 0.25 0.12
## Everton -0.01 0.27
## Fulham -0.04 0.03
## Liverpool -0.08 0.28
## Manchester City 0.60 0.56
## Manchester United 0.56 0.43
## Middlesbrough -0.52 -0.39
## Newcastle United 0.11 0.03
## Norwich City 0.05 -0.23
## Queens Park Rangers -0.14 -0.22
## Stoke City -0.33 0.01
## Sunderland -0.11 0.14
## Swansea City -0.13 0.04
## Tottenham Hotspur 0.26 0.23
## West Bromwich Albion -0.11 0.02
## Wigan Athletic -0.17 -0.16
## Wolverhampton Wanderers -0.20 -0.43
## -------
## Intercept 0.12
## Home field advantage 0.27
## Offset 1.00
Modifying parameters
In addition to fixing parameters before the model is fit, you can also manually set the parameters after the model is fitted. This can be useful if you believe the attack or defence parameters given by the model is inacurate because you have some additional information that is not captured in the data.
Modifying the parameters can also be useful for the same reasons as you might want to fit the model in two (or more) steps. For intance, if you look at the log-likelihood of the default model, and the model with the Rue-Salvesen adjustment, you will notice that they are almost identical. This is probably due to the model being nearly unidentifiable. That does not neccecarily mean that the Rue-Salvesen adjustment does not improve prediction, it might mean that it is hard to estimate the amount of adjustment based on the available data. In the paper by Rue & Salvesen they find that setting γ to 0.1 is a good choice. Here is how this can be done:
# Copy the default model.
gm_res_rs2 <- gm_res
# set the gamma parameter to 0.1
gm_res_rs2$parameters$gamma <- 0.1
# Make a few predictions to compare.
predict_result(gm_res, team1=to_predict1, team2=to_predict2, return_df = TRUE)
## team1 team2 p1 pd p2
## 1 Arsenal Fulham 0.6197118 0.2035688 0.1767195
## 2 Manchester United Chelsea 0.6736194 0.1843756 0.1420050
## 3 Bolton Wanderers Liverpool 0.2612112 0.2567706 0.4820183
predict_result(gm_res_rs2, team1=to_predict1, team2=to_predict2, return_df = TRUE)
## team1 team2 p1 pd p2
## 1 Arsenal Fulham 0.6047391 0.2083325 0.1869284
## 2 Manchester United Chelsea 0.6539250 0.1917532 0.1543218
## 3 Bolton Wanderers Liverpool 0.2780789 0.2603511 0.4615700
Evaluating predictions
To evaluate and compare the quality of predictions between different
models, you need a measure of how well the predictions match the
observed outcomes. The score_predictions function can help with this,
and implements three different scoring rules. Here is an example of how
it works. The example is done in-sample, but a proper evaluation should
compute scores out-of-sample.
my_predictions_default <- predict_result(gm_res, team1 = to_predict1, team2 = to_predict2, return_df = TRUE)
my_predictions_dc <- predict_result(gm_res_dc, team1 = to_predict1, team2 = to_predict2, return_df = TRUE)
england_2011 %>%
# Select the games to evaluate
mutate(teams_string = paste(home, '-', visitor)) %>%
filter(teams_string %in% paste(to_predict1, '-', to_predict2)) %>%
# Transform the result column (H, D, A) to reflect the column names of the prediction matrix.
mutate(result2 = c('H' = 'p1', 'D'='pd', 'A' = 'p2')[result]) %>%
select(home, visitor, result2) -> england_2011_to_score
predictions_scores_default <- score_predictions(predictions = my_predictions_default[,c('p1','pd','p2')],
observed = england_2011_to_score$result2,
score = c('log', 'brier', 'rps'))
predictions_scores_dc <- score_predictions(predictions = my_predictions_dc[,c('p1','pd','p2')],
observed = england_2011_to_score$result2,
score = c('log', 'brier', 'rps'))
# Total log score for default model, 4.27
sum(predictions_scores_default$log)
## [1] 4.273418
# Total log score for Dixon-Coles model, 4.33, which is better than for default.
sum(predictions_scores_dc$log)
## [1] 4.334922
Miscalaneous
Reverse engineering expected goals
The function predict_result that’s demonstrated above uses the underlying statistical model and the expected goals from fitted goalmodel to compute the probabilities for win-draw-lose. The function expg_from_probabilities can be used to reverse this procedure, going from probabilities to the underlying expected goals.
This can be used to reverse-engineer predictions from other models or even bookmaker odds. The expected goals extracted from this method can for example be used as input to the goalmodel to get attack and defense ratings, and to make new predictions.
Here is some code demonstrating that it can recover the expected goals which are made from a fitted goalmodel with the Poisson distribution.
# Win-Draw-Lose probabilities.
wdl_probs <- predict_result(gm_res, team1=to_predict1, team2=to_predict2,
return_df = TRUE)
wdl_probs
## team1 team2 p1 pd p2
## 1 Arsenal Fulham 0.6197118 0.2035688 0.1767195
## 2 Manchester United Chelsea 0.6736194 0.1843756 0.1420050
## 3 Bolton Wanderers Liverpool 0.2612112 0.2567706 0.4820183
expg_reveng <- expg_from_probabilities(wdl_probs[,c('p1', 'pd', 'p2')])
expg_reveng$expg
## [,1] [,2]
## [1,] 2.098592 1.0272203
## [2,] 2.278604 0.9458234
## [3,] 1.049001 1.5224753
# Compare with expected goals predictions from above.
predict_expg(gm_res_dc, team1=to_predict1, team2=to_predict2, return_df = TRUE)
## team1 team2 expg1 expg2
## 1 Arsenal Fulham 2.117914 1.0257180
## 2 Manchester United Chelsea 2.291611 0.9287833
## 3 Bolton Wanderers Liverpool 1.069024 1.5088061
You can also specify the Dixon-Coles parameter ρ (rho). Unfortunately it is not possible to extract both the two expected goals and rho, since there will typically be a large number of combinations of the three parameters that yield the same probabilities. With rho given as a constant it is however possible to extract the expected goals. Here is how you can recover the expg from the DC model fitted above, using rho = -0.13.
# Win-Draw-Lose probabilities from DC model
wdl_probs_dc <- predict_result(gm_res_dc, team1=to_predict1, team2=to_predict2,
return_df = TRUE)
wdl_probs_dc
## team1 team2 p1 pd p2
## 1 Arsenal Fulham 0.6114897 0.2268230 0.1616873
## 2 Manchester United Chelsea 0.6685679 0.2051373 0.1262948
## 3 Bolton Wanderers Liverpool 0.2522617 0.2903403 0.4573981
# Use the rho argument to specify rho.
expg_reveng_dc <- expg_from_probabilities(wdl_probs_dc[,c('p1', 'pd', 'p2')],
rho = -0.13)
expg_reveng_dc$expg
## [,1] [,2]
## [1,] 2.109461 1.0195470
## [2,] 2.282623 0.9229529
## [3,] 1.062882 1.5016514
# Compare with expected goals predictions from above.
predict_expg(gm_res_dc, team1=to_predict1, team2=to_predict2, return_df = TRUE)
## team1 team2 expg1 expg2
## 1 Arsenal Fulham 2.117914 1.0257180
## 2 Manchester United Chelsea 2.291611 0.9287833
## 3 Bolton Wanderers Liverpool 1.069024 1.5088061
Other packages
There are some other packages out there that have similar functionality as this package. The regista package implements a more general interface to fit the Dixon-Coles model. The fbRanks package implements only the Dixon-Coles time weighting functionality, but has a lot of other features. In footBayes you can fit models with Bayesian procedures. The package includes several models, including bivariate Poisson models and models with time-varying attack and defense parameters.
References
- Mark J. Dixon and Stuart G. Coles (1997) Modelling Association Football Scores and Inefficiencies in the Football Betting Market.
- Håvard Rue and Øyvind Salvesen (2001) Prediction and Retrospective Analysis of Soccer Matches in a League.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.