stan_foot: Fit football models using CmdStan

In LeoEgidi/footBayes: Fitting Bayesian and MLE Football Models

Fit football models using CmdStan

Description

Fits football goal-based models using Stan via the CmdStan backend. Supported models include: double Poisson, bivariate Poisson, Skellam, Student's t, diagonal-inflated bivariate Poisson, zero-inflated Skellam, and negative Binomial.

Usage

stan_foot(
  data,
  model,
  predict = 0,
  ranking,
  dynamic_type,
  prior_par = list(ability = normal(0, NULL), ability_sd = cauchy(0, 5), home = normal(0,
    5)),
  home_effect = TRUE,
  norm_method = "none",
  ranking_map = NULL,
  method = "MCMC",
  ...
)

Arguments

data	A data frame containing match data with columns: periods: Time point of each observation (integer >= 1). home_team: Home team's name (character string). away_team: Away team's name (character string). home_goals: Goals scored by the home team (integer >= 0). away_goals: Goals scored by the away team (integer >= 0).
model	A character string specifying the Stan model to fit. Options are: "double_pois": Double Poisson model. "biv_pois": Bivariate Poisson model. "neg_bin": Negative Binomial model. "skellam": Skellam model. "student_t": Student's t model. "diag_infl_biv_pois": Diagonal-inflated bivariate Poisson model. "zero_infl_skellam": Zero-inflated Skellam model.
predict	An integer specifying the number of out-of-sample matches for prediction. If missing, the function fits the model to the entire dataset without making predictions.
ranking	An optional "btdFoot" class element or a data frame containing ranking points for teams with the following columns: periods: Time periods corresponding to the rankings (integer >= 1). team: Team names matching those in data (character string). rank_points: Ranking points for each team (numeric).
dynamic_type	A character string specifying the type of dynamics in the model. Options are: "weekly": Weekly dynamic parameters. "seasonal": Seasonal dynamic parameters.
prior_par	A list specifying the prior distributions for the parameters of interest: ability: Prior distribution for team-specific abilities. Possible distributions are normal, student_t, cauchy, laplace. Default is normal(0, NULL). ability_sd: Prior distribution for the team-specific standard deviations. See the prior argument for more details. Default is cauchy(0, 5). home: Prior distribution for the home effect (home). Applicable only if home_effect = TRUE. Only normal priors are allowed. Default is normal(0, 5). See the rstanarm package for more details on specifying priors.
home_effect	A logical value indicating the inclusion of a home effect in the model. (default is TRUE).
norm_method	A character string specifying the method used to normalize team-specific ranking points. Options are: "none": No normalization (default). "standard": Standardization (mean 0, standard deviation 1). "mad": Median Absolute Deviation normalization. "min_max": Min-max scaling to [0,1].
ranking_map	An optional vector mapping ranking periods to data periods. If not provided and the number of ranking periods matches the number of data periods, a direct mapping is assumed.
method	A character string specifying the method used to obtain the Bayesian estimates. Options are: "MCMC": Markov chain Monte Carlo algorithm (default). "VI": Automatic differentiation variational inference algorithms. "pathfinder": Pathfinder variational inference algorithm. "laplace": Laplace algorithm.
...	Additional arguments passed to cmdstanr (e.g., iter_sampling, chains, parallel_chains).

Details

Let (y^{H}_{n}, y^{A}_{n}) denote the observed number of goals scored by the home and the away team in the n-th game, respectively. A general bivariate Poisson model allowing for goals' correlation (Karlis & Ntzoufras, 2003) is the following:

Y^H_n, Y^A_n| \lambda_{1n}, \lambda_{2n}, \lambda_{3n} \sim \mathsf{BivPoisson}(\lambda_{1n}, \lambda_{2n}, \lambda_{3n})

\log(\lambda_{1n}) = \mu+att_{h_n} + def_{a_n}

\log(\lambda_{2n}) = att_{a_n} + def_{h_n}

\log(\lambda_{3n}) =\beta_0,

where the case \lambda_{3n}=0 reduces to the double Poisson model (Baio & Blangiardo, 2010). \lambda_{1n}, \lambda_{2n} represent the scoring rates for the home and the away team, respectively, where: \mu is the home effect; the parameters att_T and def_T represent the attack and the defence abilities, respectively, for each team T, T=1,\ldots,N_T; the nested indexes h_{n}, a_{n}=1,\ldots,N_T denote the home and the away team playing in the n-th game, respectively. Attack/defence parameters are imposed a sum-to-zero constraint to achieve identifiability and assigned some weakly-informative prior distributions:

att_T \sim \mathrm{N}(\mu_{att}, \sigma_{att})

def_T \sim \mathrm{N}(\mu_{def}, \sigma_{def}),

with hyperparameters \mu_{att}, \sigma_{att}, \mu_{def}, \sigma_{def}.

Instead of using the marginal number of goals, another alternative is to modelling directly the score difference (y^{H}_{n}- y^{A}_{n}). We can use the Poisson-difference distribution (or Skellam distribution) to model goal difference in the n-th match (Karlis & Ntzoufras, 2009):

y^{H}_{n}- y^{A}_{n}| \lambda_{1n}, \lambda_{2n} \sim PD(\lambda_{1n}, \lambda_{2n}),

and the scoring rates \lambda_{1n}, \lambda_{2n} are unchanged with respect to the bivariate/double Poisson model. If we want to use a continue distribution, we can use a student t distribution with 7 degrees of freedom (Gelman, 2014):

y^{H}_{n}- y^{A}_{n} \sim t(7, ab_{h_{n}}-ab_{a(n)}, \sigma_y)

ab_t \sim \mathrm{N}(\mu + b \times {prior\_score}_t, sigma_{ab}),

where ab_t is the overall ability for the t-th team, whereas prior\_score_t is a prior measure of team's strength (for instance a ranking).

These model rely on the assumption of static parameters. However, we could assume dynamics in the attach/defence abilities (Owen, 2011; Egidi et al., 2018, Macrì Demartino et al., 2024) in terms of weeks or seasons through the argument dynamic_type. In such a framework, for a given number of times 1, \ldots, \mathcal{T}, the models above would be unchanged, but the priors for the abilities parameters at each time \tau, \tau=2,\ldots, \mathcal{T}, would be:

att_{T, \tau} \sim \mathrm{N}({att}_{T, \tau-1}, \sigma_{att})

def_{T, \tau} \sim \mathrm{N}({def}_{T, \tau-1}, \sigma_{def}),

whereas for \tau=1 we have:

att_{T, 1} \sim \mathrm{N}(\mu_{att}, \sigma_{att})

def_{T, 1} \sim \mathrm{N}(\mu_{def}, \sigma_{def}).

Of course, the identifiability constraint must be imposed for each time \tau.

The current version of the package allows for the fit of a diagonal-inflated bivariate Poisson and a zero-inflated Skellam model in the spirit of (Karlis & Ntzoufras, 2003) to better capture draw occurrences. See the vignette for further details.

Value

An object of class "stanFoot", which is a list containing:

• fit: The CmdStanFit object returned by cmdstanr.

• data: The input data.

• stan_data: The data list passed to Stan.

• stan_code: The Stan code of the underline model.

• stan_args: The optional cmdstanr parameters passed to (...).

• alg_method: The inference algorithm used to obtain the Bayesian estimates.

Author(s)

Leonardo Egidi [e-mail], Roberto Macrì Demartino [e-mail], and Vasilis Palaskas [e-mail].

References

Baio, G. and Blangiardo, M. (2010). Bayesian hierarchical model for the prediction of football results. Journal of Applied Statistics 37(2), 253-264.

Egidi, L., Pauli, F., and Torelli, N. (2018). Combining historical data and bookmakers' odds in modelling football scores. Statistical Modelling, 18(5-6), 436-459.

Gelman, A. (2014). Stan goes to the World Cup. From "Statistical Modeling, Causal Inference, and Social Science" blog.

Macrì Demartino, R., Egidi, L. and Torelli, N. Alternative ranking measures to predict international football results. Computational Statistics (2024), 1-19.

Karlis, D. and Ntzoufras, I. (2003). Analysis of sports data by using bivariate poisson models. Journal of the Royal Statistical Society: Series D (The Statistician) 52(3), 381-393.

Karlis, D. and Ntzoufras,I. (2009). Bayesian modelling of football outcomes: Using the Skellam's distribution for the goal difference. IMA Journal of Management Mathematics 20(2), 133-145.

Owen, A. (2011). Dynamic Bayesian forecasting models of football match outcomes with estimation of the evolution variance parameter. IMA Journal of Management Mathematics, 22(2), 99-113.

Examples

## Not run: 
if (instantiate::stan_cmdstan_exists()) {
  library(dplyr)

  # Example usage with ranking
  data("italy")
  italy <- as_tibble(italy)
  italy_2021 <- italy %>%
    select(Season, home, visitor, hgoal, vgoal) %>%
    filter(Season == "2021")


  teams <- unique(italy_2021$home)
  n_rows <- 20

  # Create fake ranking
  ranking <- data.frame(
    periods = rep(1, n_rows),
    team = sample(teams, n_rows, replace = FALSE),
    rank_points = sample(0:60, n_rows, replace = FALSE)
  )

  ranking <- ranking %>%
    arrange(periods, desc(rank_points))


  colnames(italy_2021) <- c("periods", "home_team", "away_team", "home_goals", "away_goals")

  fit_with_ranking <- stan_foot(
    data = italy_2021,
    model = "diag_infl_biv_pois",
    ranking = ranking,
    home_effect = TRUE,
    prior_par = list(
      ability = student_t(4, 0, NULL),
      ability_sd = cauchy(0, 3),
      home = normal(1, 10)
    ),
    norm_method = "mad",
    iter_sampling = 1000,
    chains = 2,
    parallel_chains = 2,
    adapt_delta = 0.95,
    max_treedepth = 15
  )

  # Print a summary of the model fit
  print(fit_with_ranking, pars = c("att", "def"))



  ### Use Italian Serie A from 2000 to 2002

  data("italy")
  italy <- as_tibble(italy)
  italy_2000_2002 <- italy %>%
    dplyr::select(Season, home, visitor, hgoal, vgoal) %>%
    dplyr::filter(Season == "2000" | Season == "2001" | Season == "2002")

  colnames(italy_2000_2002) <- c("periods", "home_team", "away_team", "home_goals", "away_goals")

  ### Fit Stan models
  ## no dynamics, no predictions

  fit_1 <- stan_foot(
    data = italy_2000_2002,
    model = "double_pois"
  ) # double poisson
  print(fit_1, pars = c(
    "home", "sigma_att",
    "sigma_def"
  ))

  fit_2 <- stan_foot(
    data = italy_2000_2002,
    model = "biv_pois"
  ) # bivariate poisson
  print(fit_2, pars = c(
    "home", "rho",
    "sigma_att", "sigma_def"
  ))

  fit_3 <- stan_foot(
    data = italy_2000_2002,
    mode = "skellam"
  ) # skellam
  print(fit_3, pars = c(
    "home", "sigma_att",
    "sigma_def"
  ))

  fit_4 <- stan_foot(
    data = italy_2000_2002,
    model = "student_t"
  ) # student_t
  print(fit_4, pars = c("beta"))

  ## seasonal dynamics, no prediction

  fit_5 <- stan_foot(
    data = italy_2000_2002,
    model = "double_pois",
    dynamic_type = "seasonal"
  ) # double poisson
  print(fit_5, pars = c(
    "home", "sigma_att",
    "sigma_def"
  ))

  ## seasonal dynamics, prediction for the last season

  fit_6 <- stan_foot(
    data = italy_2000_2002,
    model = "double_pois",
    dynamic_type = "seasonal",
    predict = 170
  ) # double poisson
  print(fit_6, pars = c(
    "home", "sigma_att",
    "sigma_def"
  ))

  ## other priors' options
  # double poisson with
  # student_t priors for teams abilities
  # and laplace prior for the hyper sds

  fit_p <- stan_foot(
    data = italy_2000_2002,
    model = "double_pois",
    prior_par = list(
      ability = student_t(4, 0, NULL),
      ability_sd = laplace(0, 1),
      home = normal(1, 10)
    )
  )

  print(fit_p, pars = c(
    "home", "sigma_att",
    "sigma_def"
  ))
}

## End(Not run)

LeoEgidi/footBayes documentation built on June 2, 2025, 11:32 a.m.