Fit football models using CmdStan

Description

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

Usage

stan_foot(
  data,
  model,
  predict = 0,
  ranking,
  dynamic_type,
  dynamic_weight = FALSE,
  dynamic_par = list(common_sd = FALSE, kl_variance = FALSE, periods_per_season = 2,
    spike = normal(9, 1.5), slab = normal(0, 3)),
  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. "dixon_coles": Dixon-Coles 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). The ranking points are matched to the teams of data by name, so the order of the rows does not matter. Every team of data must have exactly one value in each ranking period; teams of the ranking that are not in data are ignored.
dynamic_type	A character string specifying the type of dynamics in the model. Options are: "weekly": Weekly dynamic parameters. "seasonal": Seasonal dynamic parameters.
dynamic_weight	A logical value indicating whether to fit the weighted dynamic model of Macrì Demartino, Egidi and Torelli (2026), where the evolution of the attack and defence abilities is governed by team- and period-specific commensurate priors with spike-and-slab hyperpriors (default FALSE). It requires dynamic_type and a normal prior for the abilities (prior_par$ability = normal(...)), and it is not available for the "student_t" model. See Details.
dynamic_par	A list of options for the evolution variance of the dynamic models (ignored when dynamic_type is missing). Elements: common_sd: A logical value. If TRUE, attack and defence abilities share a single evolution standard deviation (sigma_common), as in Owen (2011). If FALSE (default), separate standard deviations sigma_att and sigma_def are estimated, as in Egidi et al. (2018). kl_variance: A logical value. If TRUE, the evolution variance is inflated in the periods that follow a summer break, as in Koopman and Lit (2015). Default FALSE. It requires dynamic_type = "seasonal" and more than periods_per_season training periods. periods_per_season: An integer (>= 2) giving the number of consecutive values of the periods column that form one season. Used only when kl_variance = TRUE. Default 2 (each season split into two halves). spike: Half-normal prior normal(location, scale) for the spike component of the commensurate precisions, used when dynamic_weight = TRUE. Default normal(9, 1.5). slab: Half-normal prior normal(location, scale) for the slab component of the commensurate precisions, used when dynamic_weight = TRUE. Default normal(0, 3). The options common_sd = TRUE, kl_variance = TRUE and dynamic_weight = TRUE are mutually exclusive.
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). In the dynamic models the prior is the distribution of the period-to-period increments of the abilities; with dynamic_weight = TRUE only normal is allowed. For student_t, cauchy and laplace the evolution standard deviations are scale parameters, not standard deviations of the increments. ability_sd: Prior distribution for the team-specific standard deviations. Possible distributions are normal, student_t, cauchy, laplace. 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 to mean 0 and standard deviation 0.5, i.e. the points minus their mean, divided by twice their standard deviation. "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}\).

The Dixon and Coles (1997) model keeps the two scoring rates \(\lambda_{1n}, \lambda_{2n}\) of the double Poisson model but multiplies the joint probability of low-scoring results by a dependence factor \(\tau\) to better capture the observed correlation in the \(0\text{-}0, 1\text{-}0, 0\text{-}1\) and \(1\text{-}1\) outcomes:

\[P(Y^H_n = y^H_n, Y^A_n = y^A_n) = \tau_{\lambda_{1n}, \lambda_{2n}}(y^H_n, y^A_n) \, \mathsf{Poisson}(y^H_n | \lambda_{1n}) \, \mathsf{Poisson}(y^A_n | \lambda_{2n}),\]

where

\[\tau_{\lambda_1, \lambda_2}(y^H, y^A) = \begin{cases} 1 - \lambda_1 \lambda_2 \rho & y^H = y^A = 0 \\ 1 + \lambda_1 \rho & y^H = 0, y^A = 1 \\ 1 + \lambda_2 \rho & y^H = 1, y^A = 0 \\ 1 - \rho & y^H = y^A = 1 \\ 1 & \text{otherwise,} \end{cases}\]

and \(\rho\) is the low-score dependence parameter, assigned a weakly-informative \(\mathrm{N}(0, 0.1)\) prior and constrained to keep the adjustment factor positive.

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, Macri 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 evolution standard deviations of the dynamic models can be specified in four alternative ways, through the arguments dynamic_par and dynamic_weight:

• Default (Egidi et al., 2018): two distinct standard deviations \(\sigma_{att}\) and \(\sigma_{def}\), shared by all the teams and periods (parameters sigma_att and sigma_def).

• dynamic_par = list(common_sd = TRUE) (Owen, 2011): a single evolution standard deviation \(\sigma_{att} = \sigma_{def} = \sigma\) (parameter sigma_common).

• dynamic_par = list(kl_variance = TRUE) (Koopman and Lit, 2015): the evolution variance is inflated in the periods that follow a summer break,

  \[\sigma^2_{k,\tau} = \sigma^2_{k} + \sigma^2_{break} \, I_{\tau}, \quad k \in \{att, def\},\]

  where \(I_{\tau} = 1\) if period \(\tau\) follows a summer break and \(0\) otherwise. The package assumes that each season is made of periods_per_season consecutive periods (default 2, i.e. two half-seasons), so that a summer break precedes the periods \(1 + m \, \mathrm{periods\_per\_season}\), \(m = 1, 2, \ldots\) (parameters sigma_att_kl, sigma_def_kl and sigma_break).

• dynamic_weight = TRUE (Macrì Demartino et al., 2026): weighted dynamic model with commensurate priors. For each team \(T\), ability \(k \in \{att, def\}\) and period \(\tau \ge 2\),

  \[k_{T, \tau} \sim \mathrm{N}(k_{T, \tau-1}, 1/\phi_{k, T, \tau}),\]

  where the commensurate precision \(\phi_{k, T, \tau}\) is assigned a two-component mixture of half-normal distributions (spike-and-slab),

  \[\phi_{k, T, \tau} \sim p_T \, \mathrm{N}^{+}(\mu_s, \psi_s) + (1 - p_T) \, \mathrm{N}^{+}(\mu_l, \psi_l), \quad p_T \sim \mathrm{Beta}(1, 1).\]

  The spike \(\mathrm{N}^{+}(\mu_s, \psi_s)\) is concentrated on large precisions (strong borrowing of information from the previous period), whereas the slab \(\mathrm{N}^{+}(\mu_l, \psi_l)\) is diffuse on small precisions (weak borrowing). The hyperparameters are set through dynamic_par$spike and dynamic_par$slab (parameters prob_spike, comm_prec_att, comm_prec_def, comm_sd_att and comm_sd_def).

To better capture draw occurrences, the package also allows for the fit of a diagonal-inflated bivariate Poisson and a zero-inflated Skellam model. The diagonal-inflated bivariate Poisson model (Karlis & Ntzoufras, 2003) is a mixture of a bivariate Poisson and a distribution over the draws:

\[P(Y^H_n = x, Y^A_n = y) = (1 - p) \, \mathsf{BivPoisson}(x, y | \lambda_{1n}, \lambda_{2n}, \lambda_{3n}) + p \, \theta_x \, I(x = y),\]

where \(p\) is the probability of the inflation component (parameter prob_of_draws, with a \(\mathrm{U}(0, 1)\) prior) and \(\theta_j\), \(j = 0, \ldots, 3\), is a discrete distribution over the draws 0-0, 1-1, 2-2 and 3-3 (parameter draw_dist, with a uniform Dirichlet prior); draws with more than three goals per team are not inflated. 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 legidi@units.it, Roberto Macrì Demartino roberto.macridemartino@deams.units.it, and Vasilis Palaskas vasilis.palaskas94@gmail.com.

References

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

Dixon, M. J. and Coles, S. G. (1997). Modelling association football scores and inefficiencies in the football betting market. Journal of the Royal Statistical Society: Series C (Applied Statistics), 46(2), 265-280.

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.

Hobbs, B. P., Carlin, B. P., Mandrekar, S. J. and Sargent, D. J. (2011). Hierarchical commensurate and power prior models for adaptive incorporation of historical information in clinical trials. Biometrics, 67(3), 1047-1056.

Koopman, S. J. and Lit, R. (2015). A dynamic bivariate Poisson model for analysing and forecasting match results in the English Premier League. Journal of the Royal Statistical Society: Series A (Statistics in Society), 178(1), 167-186.

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

Macrì Demartino, R., Egidi, L. and Torelli, N. (2026). Bayesian weighted discrete-time dynamic models for association football prediction. Journal of the Royal Statistical Society Series C: Applied Statistics, qlag032. doi:10.1093/jrsssc/qlag032

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)

  data("italy")
  italy <- as_tibble(italy)

  # Four seasons of the Italian Serie A, each split into two halves
  # (periods 1-8): the summer break precedes periods 3, 5 and 7
  italy_2018_2021 <- italy %>%
    select(Season, home, visitor, hgoal, vgoal) %>%
    filter(Season %in% c("2018", "2019", "2020", "2021")) %>%
    group_by(Season) %>%
    mutate(half = if_else(row_number() <= n() / 2, 1, 2)) %>%
    ungroup() %>%
    mutate(periods = 2 * (as.numeric(Season) - 2018) + half) %>%
    select(periods,
      home_team = home, away_team = visitor,
      home_goals = hgoal, away_goals = vgoal
    )

  # Koopman & Lit (2015): variance inflation at the summer breaks
  fit_kl <- stan_foot(
    data = italy_2018_2021,
    model = "biv_pois",
    dynamic_type = "seasonal",
    dynamic_par = list(kl_variance = TRUE, periods_per_season = 2),
    iter_sampling = 1000,
    chains = 4,
    parallel_chains = 4
  )
  print(fit_kl, pars = c("sigma_att_kl", "sigma_def_kl", "sigma_break"))

  # Weighted dynamic model with commensurate priors
  # (Macrì Demartino et al., 2026)
  fit_wdm <- stan_foot(
    data = italy_2018_2021,
    model = "double_pois",
    dynamic_type = "seasonal",
    dynamic_weight = TRUE,
    dynamic_par = list(spike = normal(9, 1.5), slab = normal(0, 3)),
    iter_sampling = 1000,
    chains = 4,
    parallel_chains = 4
  )
  print(fit_wdm, pars = "prob_spike")
}

## End(Not run)

\(\,\)