#> 1 logTie -0.114 -0.114 0.09 0.09 -0.256 0.039 1.00 405. 458.
#> 2 home 0.205 0.203 0.102 0.102 0.034 0.37 1.00 335. 366.
```
The `btd_foot` function offers flexibility in how the posterior distributions are obtained by allowing users to choose the estimation method via the `method` argument. This argument enables the user to balance computational efficiency with the precision of uncertainty estimates. The available options include:
- ```MCMC```: Uses Hamiltonian Monte Carlo (HMC) (Default).
- ```VI```: Uses the Automatic differentiation variational inference (ADVI) algorithm.
- ```pathfinder```: Uses the Pathfinder variational inference algorithm.
- ```laplace```: Uses the Laplace algorithm.
Below is an example that demonstrates how to use Variational Inference as the estimation method:
``` r
# Variational Inference Example
fit_result_vi <- btd_foot(
data = italy_2020_2021,
dynamic_rank = TRUE,
rank_measure = "mean",
method = "VI"
)
```
## Visualization tools
After fitting the model, the ```btd_foot``` function outputs an object of class ```btdFoot```. This object can be used by the ```plot_btdPosterior``` function, which is designed to visualize the posterior distributions for the log-strenghts, log-tie and home effect parameters, providing insights into the uncertainty and variability of estimates. Depending on the specifications, it can generate either boxplots or density plots of the posterior distributions. In the case of a dynamic Bayesian BTD model, the function produces a sequence of boxplots for each time period and for each team specified in the ```teams``` argument, or for all analyzed teams if this argument is omitted.
By default, ```plot_btdPosterior``` generates boxplots of the posterior distributions of teams' log-strengths. This is particularly useful for comparing the central tendency and variability of team strengths across different periods or between teams.
``` r
# Dynamic Ranking
plot_btdPosterior(fit_result_dyn)
```
``` r
# Static Ranking
plot_btdPosterior(fit_result_stat)
```
It is also possible to select specific teams to plot by using the ```teams``` argument
``` r
# Dynamic Ranking
plot_btdPosterior(fit_result_dyn,
teams = c("AC Milan", "AS Roma", "Juventus", "Inter"),
ncol = 2
)
```
``` r
# Static Ranking
plot_btdPosterior(fit_result_stat,
teams = c("AC Milan", "AS Roma", "Juventus", "Inter"),
ncol = 2
)
```
Additionally, it is possible to plot the posterior density by setting ```plot_type = "density"```
``` r
# Dynamic Ranking
plot_btdPosterior(fit_result_dyn,
teams = c("AC Milan", "AS Roma", "Juventus", "Inter"),
plot_type = "density",
scales = "free_y"
)
```
``` r
# Static Ranking
plot_btdPosterior(fit_result_stat,
teams = c("AC Milan", "AS Roma", "Juventus", "Inter"),
plot_type = "density",
scales = "free_y"
)
```
Furthermore, the ```plot_logStrength``` function allows for plotting team rankings based on different summary measures of the posterior log-strengths, including the median, mean, and maximum a posteriori (MAP) estimates used in the ```btd_foot``` function. Specifically, this function generates a ggplot object that represents the log-strength values of teams, facilitating the comparison of team performance across different periods or among various teams. Similar to the ```plot_btdPosterior``` function, specific teams can be selected using the ```teams``` argument.
``` r
# Dynamic Ranking
plot_logStrength(fit_result_dyn,
teams = c("AC Milan", "AS Roma", "Juventus", "Inter")
)
```
# Goal-based models fit
## Static fit
To start with some analysis, let's now ```italy``` data about the Italian Serie A, specifically season 2000/2001: the season consists of $T=18$ teams, we start fitting a static bivariate Poisson model using:
- The likelihood approach: the ```mle_foot``` function returns the MLE estimates along with 95\% profile-likelihood deviance confidence intervals (by default) and Wald-type confidence intervals. The user can specify the desired confidence interval with the optional argument ```interval = c("profile", "Wald")```.
- The Bayesian approach: the ```stan_foot``` function produces an Hamiltonian Monte Carlo posterior sampling by using the underlying ```rstan``` ecosystem. The user can choose the number of iterations (```iter```), the number of Markov chains (```chains```), and other optional arguments values. The output is an object of class ```stanFoot```.
At this stage, we are currently ignoring any time-dependence in our parameters, considering them to be **static** across distinct match-times.
``` r
### Use Italian Serie A 2000/2001
## with 'dplyr' environment
#
# library(dplyr)
# italy <- as_tibble(italy)
# italy_2000<- italy %>%
# dplyr::select(Season, home, visitor, hgoal,vgoal) #%>%
# dplyr::filter(Season=="2000")
# italy_2000
## alternatively, you can use the basic 'subsetting' code,
## not using the 'dplyr' environment:
data("italy")
italy <- as.data.frame(italy)
italy_2000 <- subset(
italy[, c(2, 3, 4, 6, 7)],
Season == "2000"
)
colnames(italy_2000) <- c("periods", "home_team", "away_team", "home_goals", "away_goals")
### Fit Stan models
## no dynamics, no predictions
## 4 Markov chains, 'n_iter' iterations each
n_iter <- 200 # number of MCMC iterations after the burn-in
fit1_stan <- stan_foot(
data = italy_2000,
model = "biv_pois",
chains = 4,
# parallel_chains = 4,
iter_sampling = n_iter
) # biv poisson
```
Similarly to the ```btd_foot``` function, the custom ```print``` function for an object of class ```stanFoot``` provides the usual Bayesian model summaries , including posterior means, medians, standard deviations, percentiles at 2.5\%, 25\%, 75\%, 97.5\% level, effective sample size (```n_eff```) and Gelman-Rubin statistic (```Rhat```). It accepts the same arguments as described previously.
``` r
## Print of model summary for parameters:
print(fit1_stan,
pars = c(
"home", "rho", "sigma_att",
"sigma_def", "att", "def"
),
teams = c("AC Milan", "AS Roma")
)
#> Summary of Stan football model
#> ------------------------------
#>
#> Posterior summaries for model parameters:
#> # A tibble: 8 × 10
#> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
#>
#> 1 home 0.283 0.287 0.067 0.064 0.164 0.386 1.00 809. 579.
#> 2 rho -1.75 -1.71 0.314 0.301 -2.34 -1.32 1.01 683. 380.
#> 3 sigma_att 0.207 0.201 0.065 0.061 0.108 0.318 1.02 240. 296.
#> 4 sigma_def 0.223 0.218 0.072 0.067 0.116 0.348 1.03 167. 72.9
#> 5 att[AS Roma] 0.262 0.266 0.134 0.141 0.03 0.478 1.01 541. 558.
#> 6 att[AC Milan] 0.131 0.129 0.116 0.114 -0.042 0.327 1.01 844. 422.
#> 7 def[AS Roma] -0.2 -0.187 0.147 0.14 -0.454 0.018 1.01 693. 303.
#> 8 def[AC Milan] 0 0.004 0.139 0.137 -0.239 0.227 1.00 1486. 633.
```
The **Gelman-Rubin statistic** $\hat{R}$ (```Rhat```) is below the threshold 1.1 for all the parameters, whereas the **effective sample size** (```n_eff```), measuring the approximate number of iid replications from the Markov chains, does not appear to be problematic. Thus, HMC sampling reached the convergence.
As we could expect, there is a positive effect from the home-effect (posterior mean about 0.3), and this implies that if two teams are equally good (meaning that their attack and defence abilities almost coincide), assuming that the constant intercept $\mu \approx 0$, we get that the average number of goals for the home-team will be $\lambda_{1} = \exp \{0.3 \} \approx 1.35$, against $\lambda_{2} = \exp \{0 \} = 1$.
In the model above, we are assuming that the covariance $\lambda_{3n}$ is constant and not depending on the match and/or on teams characteristics/further covariates:
\begin{align}
\lambda_{3n} =&\ \exp\{\rho\}\\
\rho \sim & \ \mathrm{N}^+(0,1),\\
\end{align}
where $\rho$ is assigned an half-Gaussian distribution with standard deviation equal to 1. According to the fit above, this means that in the model above we get an estimate of $\lambda_{3n}= \exp\{-4.25\} \approx 0.014$, suggesting a low, despite non-null, amount of goals-correlation existing in the 2000/2001 Italian Serie A. Of course, in the next package's version, the user will be allowed to specify a more general linear predictor for $\log(\lambda_{3n})$, as outlined in the BP presentation above, along with some prior distributions for the parameters involved in the covariance formulation.
We can also depict the marginal posterior distributions for $\rho$ (and eventually for the other fixed-effects parameters) using the ```bayesplot``` package for Bayesian visualizations:
``` r
## Marginal posterior with bayesplot
posterior1 <- fit1_stan$fit$draws(format = "matrix")
mcmc_areas(posterior1, pars = c(
"home", "rho",
"sigma_att", "sigma_def"
)) +
theme_bw()
```
We can also access the original BP Stan code by typing:
``` r
### Model's code extraction
fit1_stan$stan_code
#> [1] "functions{"
#> [2] ""
#> [3] " real bipois_lpmf(array[] int r , real mu1,real mu2,real mu3) {"
#> [4] " real ss;"
#> [5] " real log_s;"
#> [6] " real mus;"
#> [7] " int miny;"
#> [8] ""
#> [9] " miny = min(r[1], r[2]);"
#> [10] ""
#> [11] " ss = poisson_lpmf(r[1] | mu1) + poisson_lpmf(r[2] | mu2) -"
#> [12] " exp(mu3);"
#> [13] " if(miny > 0) {"
#> [14] " mus = -mu1-mu2+mu3;"
#> [15] " log_s = ss;"
#> [16] ""
#> [17] " for(k in 1:miny) {"
#> [18] " log_s = log_s + log(r[1] - k + 1) + mus"
#> [19] " + log(r[2] - k + 1)"
#> [20] " - log(k);"
#> [21] " ss = log_sum_exp(ss, log_s);"
#> [22] " }"
#> [23] " }"
#> [24] " return(ss);"
#> [25] " }"
#> [26] ""
#> [27] " }"
#> [28] " data{"
#> [29] " int N; // number of games"
#> [30] " int N_prev;"
#> [31] " array[N,2] int y;"
#> [32] " int nteams;"
#> [33] " array[N] int instants_rank;"
#> [34] " int ntimes_rank; // dynamic periods for ranking"
#> [35] " array[N] int team1;"
#> [36] " array[N] int team2;"
#> [37] " array[N_prev]int team1_prev;"
#> [38] " array[N_prev] int team2_prev;"
#> [39] " matrix[ntimes_rank,nteams] ranking;"
#> [40] " int ind_home;"
#> [41] " real mean_home; // Mean for home effect"
#> [42] " real sd_home; // Standard deviation for home effect"
#> [43] ""
#> [44] " // priors part"
#> [45] " int prior_dist_num; // 1 gaussian, 2 t, 3 cauchy, 4 laplace"
#> [46] " int prior_dist_sd_num; // 1 gaussian, 2 t, 3 cauchy, 4 laplace"
#> [47] ""
#> [48] " real hyper_df;"
#> [49] " real hyper_location;"
#> [50] ""
#> [51] " real hyper_sd_df;"
#> [52] " real hyper_sd_location;"
#> [53] " real hyper_sd_scale;"
#> [54] " }"
#> [55] " parameters{"
#> [56] " vector[nteams] att_raw;"
#> [57] " vector[nteams] def_raw;"
#> [58] " real sigma_att;"
#> [59] " real sigma_def;"
#> [60] " real home;"
#> [61] " real rho;"
#> [62] " real gamma;"
#> [63] " }"
#> [64] " transformed parameters{"
#> [65] " real adj_h_eff; // Adjusted home effect"
#> [66] " vector[nteams] att;"
#> [67] " vector[nteams] def;"
#> [68] " array[N] vector[3] theta;"
#> [69] ""
#> [70] " for (t in 1:nteams){"
#> [71] " att[t] = att_raw[t]-mean(att_raw);"
#> [72] " def[t] = def_raw[t]-mean(def_raw);"
#> [73] " }"
#> [74] ""
#> [75] " adj_h_eff = home * ind_home;"
#> [76] ""
#> [77] " for (n in 1:N){"
#> [78] " theta[n,1] = exp(adj_h_eff+att[team1[n]]+def[team2[n]]+"
#> [79] " (gamma/2)*(ranking[instants_rank[n], team1[n]]-ranking[instants_rank[n], team2[n]]));"
#> [80] " theta[n,2] = exp(att[team2[n]]+def[team1[n]]-"
#> [81] " (gamma/2)*(ranking[instants_rank[n], team1[n]]-ranking[instants_rank[n], team2[n]]));"
#> [82] " theta[n,3] = exp(rho);"
#> [83] " }"
#> [84] " }"
#> [85] " model{"
#> [86] " // log-priors for team-specific abilities"
#> [87] " for (t in 1:(nteams)){"
#> [88] " if (prior_dist_num == 1){"
#> [89] " target+= normal_lpdf(att_raw[t]|hyper_location, sigma_att);"
#> [90] " target+= normal_lpdf(def_raw[t]|hyper_location, sigma_def);"
#> [91] " }"
#> [92] " else if (prior_dist_num == 2){"
#> [93] " target+= student_t_lpdf(att_raw[t]|hyper_df, hyper_location, sigma_att);"
#> [94] " target+= student_t_lpdf(def_raw[t]|hyper_df, hyper_location, sigma_def);"
#> [95] " }"
#> [96] " else if (prior_dist_num == 3){"
#> [97] " target+= cauchy_lpdf(att_raw[t]|hyper_location, sigma_att);"
#> [98] " target+= cauchy_lpdf(def_raw[t]|hyper_location, sigma_def);"
#> [99] " }"
#> [100] " else if (prior_dist_num == 4){"
#> [101] " target+= double_exponential_lpdf(att_raw[t]|hyper_location, sigma_att);"
#> [102] " target+= double_exponential_lpdf(def_raw[t]|hyper_location, sigma_def);"
#> [103] " }"
#> [104] " }"
#> [105] ""
#> [106] ""
#> [107] " // log-hyperpriors for sd parameters"
#> [108] " if (prior_dist_sd_num == 1 ){"
#> [109] " target+=normal_lpdf(sigma_att|hyper_sd_location, hyper_sd_scale);"
#> [110] " target+=normal_lpdf(sigma_def|hyper_sd_location, hyper_sd_scale);"
#> [111] " }"
#> [112] " else if (prior_dist_sd_num == 2){"
#> [113] " target+=student_t_lpdf(sigma_att|hyper_sd_df, hyper_sd_location, hyper_sd_scale);"
#> [114] " target+=student_t_lpdf(sigma_def|hyper_sd_df, hyper_sd_location, hyper_sd_scale);"
#> [115] " }"
#> [116] " else if (prior_dist_sd_num == 3){"
#> [117] " target+=cauchy_lpdf(sigma_att|hyper_sd_location, hyper_sd_scale);"
#> [118] " target+=cauchy_lpdf(sigma_def|hyper_sd_location, hyper_sd_scale);"
#> [119] " }"
#> [120] " else if (prior_dist_sd_num == 4){"
#> [121] " target+=double_exponential_lpdf(sigma_att|hyper_sd_location, hyper_sd_scale);"
#> [122] " target+=double_exponential_lpdf(sigma_def|hyper_sd_location, hyper_sd_scale);"
#> [123] " }"
#> [124] ""
#> [125] " // log-priors fixed effects"
#> [126] " target+=normal_lpdf(home|mean_home,sd_home);"
#> [127] " target+=normal_lpdf(rho|0,1);"
#> [128] " target+=normal_lpdf(gamma|0,1);"
#> [129] ""
#> [130] " // likelihood"
#> [131] " for (n in 1:N){"
#> [132] " //target+=bipois_lpmf(y[n,]| theta[n,1],"
#> [133] " // theta[n,2], theta[n,3]);"
#> [134] " target+=poisson_lpmf(y[n,1]| theta[n,1]+theta[n,3]);"
#> [135] " target+=poisson_lpmf(y[n,2]| theta[n,2]+theta[n,3]);"
#> [136] " }"
#> [137] " }"
#> [138] " generated quantities{"
#> [139] " array[N,2]int y_rep;"
#> [140] " array[N_prev,2] int y_prev;"
#> [141] " array[N_prev] vector[3] theta_prev;"
#> [142] " vector[N] log_lik;"
#> [143] " array[N] int diff_y_rep;"
#> [144] ""
#> [145] " //in-sample replications"
#> [146] " for (n in 1:N){"
#> [147] " y_rep[n,1] = poisson_rng(theta[n,1]+theta[n,3]);"
#> [148] " y_rep[n,2] = poisson_rng(theta[n,2]+theta[n,3]);"
#> [149] " diff_y_rep[n] = y_rep[n,1] - y_rep[n,2];"
#> [150] " log_lik[n] = poisson_lpmf(y[n,1]| theta[n,1]+theta[n,3])+"
#> [151] " poisson_lpmf(y[n,2]| theta[n,2]+theta[n,3]);"
#> [152] " //bipois_lpmf(y[n,]| theta[n,1],"
#> [153] " // theta[n,2], theta[n,3]);"
#> [154] " }"
#> [155] ""
#> [156] " //out-of-sample predictions"
#> [157] " if (N_prev > 0) {"
#> [158] " for (n in 1:N_prev){"
#> [159] " theta_prev[n,1] = exp(adj_h_eff+att[team1_prev[n]]+"
#> [160] " def[team2_prev[n]]+"
#> [161] " (gamma/2)*(ranking[instants_rank[N],team1_prev[n]]-ranking[instants_rank[N],team2_prev[n]]));"
#> [162] " theta_prev[n,2] = exp(att[team2_prev[n]]+"
#> [163] " def[team1_prev[n]]-"
#> [164] " (gamma/2)*(ranking[instants_rank[N],team1_prev[n]]-ranking[instants_rank[N],team2_prev[n]]));"
#> [165] " theta_prev[n,3] = exp(rho);"
#> [166] " y_prev[n,1] = poisson_rng(theta_prev[n,1]+theta_prev[n,3]);"
#> [167] " y_prev[n,2] = poisson_rng(theta_prev[n,2]+theta_prev[n,3]);"
#> [168] " }"
#> [169] " }"
#> [170] " }"
```
In addition the `stan_foot` function provides flexibility in obtaining posterior distributions by allowing users to choose the estimation method via the `method` argument. The available options include:
- ```MCMC```: Uses Hamiltonian Monte Carlo (HMC) (Default).
- ```VI```: Uses the Automatic differentiation variational inference (ADVI) algorithm.
- ```pathfinder```: Uses the Pathfinder variational inference algorithm.
- ```laplace```: Uses the Laplace algorithm.
Below is an example that demonstrates how to use Pathfinder algorithm as the estimation method:
``` r
# Pathfinder algorithm example
fit1_stan_path <- stan_foot(
data = italy_2000,
model = "biv_pois",
method = "pathfinder"
) # biv poisson
```
We fit now the same model under the MLE approach with Wald-type confidence intervals. We can then print the MLE estimates, e.g for the parameters $\rho$ and $\text{home}$:
``` r
### Fit MLE models
## no dynamics, no predictions
## Wald intervals
fit1_mle <- mle_foot(
data = italy_2000,
model = "biv_pois",
interval = "Wald"
) # mle biv poisson
fit1_mle$home_effect
#> 2.5% mle 97.5%
#> [1,] 0.2 0.3 0.39
```
We got a very similar estimate to the Bayesian model for the home-effect.
### Changing default priors
One of the common practices in Bayesian statistics is to change the priors and perform some **sensitivity tests**. The default priors for the team-specific abilities and their related team-level standard deviations are:
\begin{align}
\text{att}_t &\sim \mathrm{N}(\mu_{\text{att}}, \sigma_{\text{att}}),\\
\text{def}_t &\sim \mathrm{N}(\mu_{\text{def}}, \sigma_{\text{def}}),\\
\sigma_{\text{att}}, \sigma_{\text{def}} &\sim \mathsf{Cauchy}^+(0,5),
\end{align}
where $\mathsf{Cauchy}^+$ denotes the half-Cauchy distribution with support $[0, +\infty)$. However, the user is free to elicit some different priors, possibly choosing one among the following distributions: Gaussian (```normal```), student-$t$ (```student_t```), Cauchy (```cauchy```) and Laplace (```laplace```). The ```ability``` optional argument allows to specify the priors for the team-specific parameters $\text{att}$ and $\text{def}$, whereas the optional argument ```ability_sd``` allows to assign a prior to the group-level standard deviations $\sigma_{\text{att}}, \sigma_{\text{def}}$. For instance, for each team $t, t=1,\ldots,T$, we could consider:
\begin{align}
\text{att}_t &\sim t(4, \mu_{\text{att}}, \sigma_{\text{att}}),\\
\text{def}_t &\sim t(4, \mu_{\text{def}}, \sigma_{\text{def}}),\\
\sigma_{\text{att}}, \sigma_{\text{def}} &\sim \mathsf{Laplace}^+(0,1),
\end{align}
where $t(\text{df}, \mu, \sigma)$ denotes a student-$t$ distribution with df degrees of freedom, location $\mu$ and scale $\sigma$, whereas $\mathsf{Laplace}^+$ denotes a half-Laplace distribution.
``` r
### Fit Stan models
## changing priors
## student-t for team-specific abilities, laplace for sds
fit1_stan_t <- stan_foot(
data = italy_2000,
model = "biv_pois",
chains = 4,
prior_par = list(
ability = student_t(4, 0, NULL),
ability_sd = laplace(0, 1),
home = normal(0, 10)
),
# parallel_chains = 4,
iter_sampling = n_iter
) # biv poisson
```
Then, we can compare the marginal posteriors from the two models, the one with Gaussian team-specific abilities and the default $\mathsf{Cauchy}(0,5)$ for the team-level sds, and the other one specified above, with student-$t$ distributed team-specific abilities and the $\mathsf{Laplace}^+(0,1)$ for the team-level sds. We depict here the comparison for the attack team-level sds only:
``` r
## comparing posteriors
posterior1_t <- fit1_stan_t$fit$draws(format = "matrix")
model_names <- c("Default", "Stud+Laplace")
color_scheme_set(scheme = "gray")
gl_posterior <- cbind(
posterior1[, "sigma_att"],
posterior1_t[, "sigma_att"]
)
colnames(gl_posterior) <- c("sigma_att", "sigma_att_t")
mcmc_areas(gl_posterior, pars = c("sigma_att", "sigma_att_t")) +
xaxis_text(on = TRUE, size = ggplot2::rel(2.9)) +
yaxis_text(on = TRUE, size = ggplot2::rel(2.9)) +
scale_y_discrete(labels = ((parse(text = model_names)))) +
ggtitle("Att/def sds") +
theme(plot.title = element_text(hjust = 0.5, size = rel(2.6))) +
theme_bw()
```
The student+laplace prior induces a lower amount of group-variability in the $\sigma_{\text{att}}$ marginal posterior distribution (then, a larger shrinkage towards the grand mean $\mu_{\text{att}}$).
When specifying the prior for the team-specific parameters through the argument ```ability```, you are not allowed to fix the group-level standard deviations $\sigma_{\text{att}}, \sigma_{\text{def}}$ to some numerical values. Rather, they need to be assigned a reasonable prior distribution.
For such reason, the most appropriate specification for the ```ability``` argument is ```ability = 'dist'(0, NULL)```, where the scale argument is set to ```NULL``` (otherwise, a warning message is occurring).
## Dynamic fit
A structural limitation in the previous models is the assumption of static team-specific parameters, namely teams are assumed to have a constant performance across time, as determined by the attack and defence abilities (att, def). However, teams' performance tend to be *dynamic* and change across different years, if not different weeks. Many factors contribute to this football aspect:
(i) teams act during the summer/winter players' transfermarket, by dramatically changing their rosters;
(ii) some teams' players could be injured in some periods, by affecting the global quality of the team in some matches;
(iii) coaches could be dismissed from their teams due to some non satisfactory results;
(iv) some teams could improve/worsen their attitudes due to the so-called turnover;
and many others. Again, we could assume **dynamics** in the attach/defence
abilities as in @owen2011dynamic , @egidi2018combining and @demartino2024alternative in terms of weeks or seasons. In such 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 **auto-regressive** of order 1:
\begin{align}
\text{att}_{t, \tau} & \sim \mathrm{N}({\text{att}}_{t, \tau-1}, \sigma_{\text{att}})\\
\text{def}_{t, \tau} &\sim \mathrm{N}({\text{def}}_{t, \tau-1}, \sigma_{\text{def}}),
\end{align}
whereas for $\tau=1$ we have:
\begin{align}
\text{att}_{t, 1} & \sim \mathrm{N}(\mu_{\text{att}}, \sigma_{\text{att}})\\
\text{def}_{t, 1} &\sim \mathrm{N}(\mu_{\text{def}}, \sigma_{\text{def}}),
\end{align}
with hyperparameters $\mu_{\text{att}}, \mu_{\text{def}}, \sigma_{\text{att}}, \sigma_{\text{def}}$ and with the standard deviations capturing the time's/evolution's variation of both the teams' skills (here assumed to be constant across time and teams).
Of course, the identifiability constraint must be imposed for each time $\tau$.
We can use the ```dynamic_type``` argument in the ```stan_foot``` function, with possible options ```'seasonal'``` or ```'weekly'``` in order to consider more seasons (no examples are given in this course) or more week-times within a single season, respectively. Let's fit a weekly-dynamic parameters model on the Serie A 2000/2001 season:
``` r
### Fit Stan models
## seasonal dynamics, no predictions
## 2 Markov chains, 'n_iter' iterations each
fit2_stan <- stan_foot(
data = italy_2000,
model = "biv_pois",
dynamic_type = "weekly",
# parallel_chains = 2,
chains = 2,
iter_sampling = n_iter
) # biv poisson
```
``` r
print(fit2_stan, pars = c(
"home", "rho", "sigma_att",
"sigma_def"
))
#> Summary of Stan football model
#> ------------------------------
#>
#> Posterior summaries for model parameters:
#> # A tibble: 4 × 10
#> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
#>
#> 1 home 0.28 0.286 0.069 0.072 0.168 0.387 1.00 521. 285.
#> 2 rho -1.76 -1.71 0.367 0.301 -2.39 -1.30 1.00 460. 251.
#> 3 sigma_att 0.048 0.047 0.009 0.007 0.035 0.066 1.22 7.20 11.4
#> 4 sigma_def 0.058 0.056 0.012 0.012 0.042 0.08 1.30 5.56 24.1
```
From the printed summary, we may note that the degree of goals' correlation seems to be again very small here. Moreover, the Gelman-Rubin statistic for $\sigma_{\text{att}}$ is relatively high, whereas the effective sample sizes for $\sigma_{\text{att}}$ and $\sigma_{\text{def}}$ are quite low. This is suggesting possible inefficiencies during the HMC sampling and that a *model-reparametrization* could be suited and effective at this stage. Another option is to play a bit with the prior specification for $\sigma_{\text{att}}$ and $\sigma_{\text{def}}$, for instance by specifying a prior inducing less shrinkage in the team-specific abilities.
To deal with these issues and broaden the set of candidate models, let's fit also a dynamic double-Poisson model with the ```double_pois``` option for the argument ```model```:
``` r
### Fit Stan models
## weekly dynamics, no predictions
## 2 chains, 'n_iter' iterations each
fit3_stan <- stan_foot(
data = italy_2000,
model = "double_pois",
dynamic_type = "weekly",
# parallel_chains = 2,
chains = 2,
iter_sampling = n_iter
) # double poisson
```
``` r
print(fit3_stan, pars = c(
"home", "sigma_att",
"sigma_def"
))
#> Summary of Stan football model
#> ------------------------------
#>
#> Posterior summaries for model parameters:
#> # A tibble: 3 × 10
#> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
#>
#> 1 home 0.407 0.408 0.052 0.056 0.32 0.494 1.01 412. 231.
#> 2 sigma_att 0.049 0.046 0.014 0.015 0.032 0.075 1.85 3.04 36.9
#> 3 sigma_def 0.053 0.05 0.011 0.011 0.039 0.074 1.31 7.03 19.3
```
The fitting problems mentioned above remain also for the double Poisson model...Thus, it's time to play a little bit with the prior distributions. Also in the dynamic approach we can change the default priors for the team-specific abilities and their standard deviations, respectively, through the optional arguments ```ability``` and ```ability_sd```. The specification follows almost analogously the static case: with the first argument we may specify the prior's family for the team-specific abilities and the specific priors for $\text{att}_{t,1}, \text{def}_{t,1}$ along with the hyper-prior location $\mu_{\text{att}}, \mu_{\text{def}}$, whereas $\sigma_{\text{att}}$ and $\sigma_{\text{def}}$ need to be assigned some proper prior distribution. Assume to fit the same double Poisson model, but here we suppose student-$t$ distributed team-specific abilities with 4 degrees of freedom to eventually capture more extreme team-specific abilities (more variability, i.e. less shrinkage), along with a $\text{Cauchy}^+(0,25)$ for their standard deviations (to better capture a possible larger evolution variability):
\begin{align}
\text{att}_{t, \tau} &\ \sim t(4, {\text{att}}_{t, \tau-1}, \sigma_{\text{att}})\\
\text{def}_{t, \tau} &\ \sim t(4, {\text{def}}_{t, \tau-1}, \sigma_{\text{def}})\\
\sigma_{\text{att}}, \sigma_{\text{def}} &\ \sim \mathsf{Cauchy}^+(0,25),
\end{align}
``` r
### Fit Stan models
## weekly dynamics, no predictions
## 2 chains, 'n_iter' iterations each
fit3_stan_t <- stan_foot(
data = italy_2000,
model = "double_pois",
prior_par = list(
ability = student_t(4, 0, NULL),
ability_sd = cauchy(0, 25),
home = normal(0, 5)
),
dynamic_type = "weekly",
# parallel_chains = 2,
chains = 2,
iter_sampling = n_iter
) # double poisson
```
``` r
print(fit3_stan_t, pars = c(
"home", "sigma_att",
"sigma_def"
))
#> Summary of Stan football model
#> ------------------------------
#>
#> Posterior summaries for model parameters:
#> # A tibble: 3 × 10
#> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
#>
#> 1 home 0.41 0.411 0.043 0.041 0.333 0.48 1.03 131. 208.
#> 2 sigma_att 0.033 0.032 0.008 0.009 0.022 0.046 1.47 3.98 53.3
#> 3 sigma_def 0.037 0.04 0.015 0.019 0.014 0.06 2.10 2.77 22.4
```
As we may conclude, the situation has been only slightly improved.
# Using ```btd_foot``` and ```stan_foot``` sequentially
One of the main features of the ```footBayes``` package is the ability to first compute the log-strengths using the Bayesian Bradley-Terry-Davidson model through the ```btd_foot``` function and then fit a Bayesian goal-based model using the ```stan_foot``` by passing the previously estimated historical team strengths as an additional covariate via the ```ranking``` argument. Specifically, this argument accepts an element of class ```btdFoot```or alternatively a data frame 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 (relative strengths) for each team (numeric).
The ability to use either ```btdFoot``` objects or custom ranking data frames provides flexibility in how the user can model team strengths, allowing for the integration of various ranking systems, such as the [FIFA ranking](https://inside.fifa.com/fifa-world-ranking).
``` r
# Dynamic Bradley-Terry-Davidson model
data("italy")
italy_2020_2021_rank <- italy %>%
dplyr::select(Season, home, visitor, hgoal, vgoal) %>%
dplyr::filter(Season == "2020" | Season == "2021") %>%
dplyr::mutate(match_outcome = dplyr::case_when(
hgoal > vgoal ~ 1, # Home team wins
hgoal == vgoal ~ 2, # Draw
hgoal < vgoal ~ 3 # Away team wins
)) %>%
dplyr::filter(dplyr::row_number() <= 570) %>%
dplyr::mutate(periods = dplyr::case_when(
dplyr::row_number() <= 190 ~ 1,
dplyr::row_number() <= 380 ~ 2,
dplyr::row_number() <= 570 ~ 3
)) %>%
dplyr::select(periods,
home_team = home,
away_team = visitor, match_outcome
)
fit_btd_dyn <- btd_foot(
data = italy_2020_2021_rank,
dynamic_rank = TRUE,
rank_measure = "median",
iter_sampling = 1000,
# parallel_chains = 2,
chains = 2,
adapt_delta = 0.9,
max_treedepth = 12
)
# Dynamic Bivariate Poisson Model
italy_2020_2021_fit <- italy %>%
dplyr::select(Season, home, visitor, hgoal, vgoal) %>%
dplyr::filter(Season == "2020" | Season == "2021") %>%
dplyr::filter(dplyr::row_number() <= 570) %>%
dplyr::mutate(periods = dplyr::case_when(
dplyr::row_number() <= 190 ~ 1,
dplyr::row_number() <= 380 ~ 2,
dplyr::row_number() <= 570 ~ 3
)) %>% # Assign periods based on match number
dplyr::select(periods,
home_team = home,
away_team = visitor, home_goals = hgoal, away_goals = vgoal
)
fit_stan_rank <- stan_foot(
data = italy_2020_2021_fit,
model = "biv_pois",
ranking = fit_btd_dyn,
predict = 180,
prior_par = list(
ability = student_t(4, 0, NULL),
ability_sd = cauchy(0, 25),
home = normal(0, 5)
),
dynamic_type = "season",
chains = 2,
# parallel_chains = 2,
iter_sampling = 1000
)
```
``` r
print(fit_stan_rank,
pars = c("home", "rho", "sigma_att", "sigma_def")
)
#> Summary of Stan football model
#> ------------------------------
#>
#> Posterior summaries for model parameters:
#> # A tibble: 4 × 10
#> variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
#>
#> 1 home 0.159 0.159 0.068 0.064 0.047 0.272 1 2342. 1669.
#> 2 rho -1.10 -1.09 0.162 0.163 -1.38 -0.849 1 2388. 1576.
#> 3 sigma_att 0.025 0.02 0.02 0.015 0.004 0.064 1.06 16.1 42.5
#> 4 sigma_def 0.03 0.026 0.02 0.022 0.005 0.066 1.18 12.9 44.1
```
# Model estimates and visualization tools
## Plotting and interpreting team-specific abilities
Once the model has been fitted, there is a large amount of interesting summaries to explore. The function ```foot_abilities``` allows to depict posterior/confidence intervals for global attack and defense abilities on the considered data (attack abilities are plotted in red, whereas defense abilities in blue colors). The *higher the attack and the lower the defence* for a given team, and *the better is the overall team's strength*.
We can produce the team-specific abilities for the two static fits above, ```fit1_mle``` (MLE) and ```fit1_stan``` (Bayes), with red bars for the attack and blue bars for the defence, respectively:
``` r
## Plotting abilities: credible and confidence 95% intervals
foot_abilities(object = fit1_stan, data = italy_2000)
foot_abilities(object = fit1_mle, data = italy_2000)
```
AS Roma, the team winning the Serie A 2000/2001, is associated with the highest attack ability and the lowest defence ability according to both the models. In general, the models seem to well capture the static abilities: AS Roma, Lazio Roma and Juventus (1st, 3rd and 2nd at the end of that season, respectively) are rated as the best teams in terms of their abilities, whereas AS Bari, SSC Napoli and Vicenza Calcio (all relegated at the end of the season) have the worst abilities.
We can also depict the team-specific dynamic plots for the dynamic models:
``` r
## Plotting abilities: credible and confidence 95% intervals
foot_abilities(fit2_stan, italy_2000)
```
As we can see, dynamic abilities naturally evolve over the time: better teams (AS Roma, Lazio Roma, Juventus, Parma) are associated with increasing attack abilities and decreasing defence abilities, whereas the worst ones (AS Bari, SSC NApoli, and Hellas Verona) exhibit decreasing attacking skills and increasing defensive skills. The reason for these increasing/decreasing behaviours is straightforward: at the beginning, all the attack/defence parameters have been initialized to have location equal to 0. The user is free to change the location, and in the final package's version he will also have the possibility to elicit different team-specific hyper-prior locations.
# Model checking
Checking the model fit is a relevant and vital statistical task. To this purpose, we can evaluate hypothetical replications $\mathcal{D}^{\text{rep}}$ under the **posterior predictive distribution**
$$p(\mathcal{D}^{\text{rep}}| \mathcal{D}) = \int p(\mathcal{D}^{\text{rep}}| \boldsymbol{\theta}) \pi(\boldsymbol{\theta}| \mathcal{D}) d\boldsymbol{\theta},$$
and check whether these replicated values are somehow close to the observed data $\mathcal{D}$. These methods comparing hypothetical replications with the observed data are named **posterior predictive checks** and have great theoretical and applied appeal in Bayesian inference. See @gelman2014bayesian for an overview.
The function ```pp_foot``` allows to obtain:
- an aggregated plot depicting the observed frequencies of the observed goal differences $Z_n=X_n-Y_n, \ n=1,\ldots,N$ plotted against the replicated ones;
- a visualization of the match-ordered goal differences along with their 50\% and 95\% credible intervals.
``` r
## PP checks: aggregated goal's differences and ordered goal differences
pp_foot(
object = fit1_stan, data = italy_2000,
type = "aggregated"
)
#> $pp_plot
```
```
#>
#> $pp_table
#> goal diff. Bayesian p-value
#> 1 -3 0.326
#> 2 -2 0.870
#> 3 -1 0.988
#> 4 0 0.041
#> 5 1 0.265
#> 6 2 0.109
#> 7 3 0.924
pp_foot(
object = fit1_stan, data = italy_2000,
type = "matches"
)
#> $pp_plot
```
```
#>
#> $pp_table
#> 1-alpha emp. coverage
#> 1 0.95 0.944
```
The aggregated goal difference frequencies seem to be decently captured by the model's replications: in the first plot, blue horizontal lines denote the observed goal differences frequencies registered in the dataset, whereas yellow jittered points denote the correspondent replications. Goal-difference of 0, corresponding to the draws occurrences, is only slightly underestimated by the model. However, in general there are no particular clues of model's misfit.
In the second plot, the ordered observed goal differences are plotted against their replications (50% and 95% credible intervals), and from this plot also we do not have particular signs of model's misfits.
Other useful PP checks, such as the overlap between data density and replicated data densities to check eventual inconsistencies, can be obtained through the standard use of the ```bayesplot``` package, in this case providing an approximation to a continuous distribution using an input kernel choice (```bw = 0.5``` in the ```ppc_dens_overlay``` used below):
``` r
## PPC densities overlay with the bayesplot package
# extracting the replications
draws_raw <- fit1_stan$fit$draws()
draws <- posterior::as_draws_rvars(draws_raw)
sims <- list()
sims$y_rep <- posterior::draws_of(draws[["y_rep"]])
goal_diff <- italy_2000$home_goals - italy_2000$away_goals
# plotting data density vs replications densities
ppc_dens_overlay(goal_diff, sims$y_rep[, , 1] - sims$y_rep[, , 2], bw = 0.5) +
theme_bw()
```
From this plot above we have the empirical confirmation that the goal difference is well captured by the static bivariate Poisson model.
# Predictions and predictive accuracy
## Posterior out-of-sample probabilities
The hottest feature in sports analytics is to obtain future predictions. By considering the **posterior predictive distribution** for future and observable data $\tilde{\mathcal{D}}$, we acknowledge the whole model's prediction uncertainty (which propagates from the posterior model's uncertainty) and we can then generate observable values $\tilde{D}$ conditioned on the posterior model's parameters estimates:
$$p(\tilde{\mathcal{D}}| \mathcal{D}) = \int p(\tilde{\mathcal{D}}| \boldsymbol{\theta}) \pi(\boldsymbol{\theta}| \mathcal{D}) d\boldsymbol{\theta}.$$
We may then predict test set matches by using the argument ```predict``` of the ```stan_foot``` function, for instance considering the last four weeks of the 2000/2001 season as the test set, and then computing posterior-results probabilities using the function ```foot_prob``` for two teams belonging to the test set, such as Reggina Calcio and AC Milan:
``` r
### Fit Stan models
## weekly dynamics, predictions of last four weeks
## 2 chains 'n_iter' iterations each
fit4_stan <- stan_foot(
data = italy_2000,
model = "biv_pois",
predict = 36,
dynamic_type = "weekly",
# parallel_chains = 2,
chains = 2,
iter_sampling = n_iter
) # biv poisson
```
``` r
foot_prob(
object = fit4_stan, data = italy_2000,
home_team = "Reggina Calcio",
away_team = "AC Milan"
)
#> $prob_table
#> home_team away_team prob_h prob_d prob_a mlo
#> 1 Reggina Calcio AC Milan 0.323 0.248 0.43 1-1 (0.101)
#>
#> $prob_plot
```
Darker regions are associated with higher posterior probabilities, whereas the red square corresponds to the actual observed result, 2-1 for Reggina Calcio. According to the posterior-results probabilities, this final observed result had in principle about a 5% probability to happen! (Remember, football is about rare events...).
## Home-win posterior probabilities
We can also use the out-of-sample posterior-results probabilities to compute some aggregated **home/draw/loss** probabilities (based then on the $S$ draws from the MCMC method) for a given match:
\begin{align}
p_{\text{home}}= &\ \textrm{Pr}(X>Y) = \frac{1}{S}\sum_{s=1}^S| \tilde{x}^{(s)}> \tilde{y}^{(s)}|\\
p_{\text{draw}} = &\ \textrm{Pr}(X=Y) = \frac{1}{S}\sum_{s=1}^S| \tilde{x}^{(s)}= \tilde{y}^{(s)}|\\
p_{\text{loss}} = &\ \textrm{Pr} (X $round_table
#> Home Away Home_prob Observed
#> 1 Hellas Verona Bologna FC 0.338 -
#> 2 Inter Bologna FC 0.438 -
#> 3 Udinese Calcio SSC Napoli 0.495 -
#> 4 ACF Fiorentina SSC Napoli 0.535 -
#> 5 US Lecce Lazio Roma 0.235 -
#> 6 Inter Lazio Roma 0.295 -
#> 7 Brescia Calcio AS Bari 0.627 -
#> 8 Reggina Calcio AS Bari 0.590 -
#> 9 Udinese Calcio Vicenza Calcio 0.392 -
#> 10 Brescia Calcio Vicenza Calcio 0.510 -
#> 11 AS Roma Parma AC 0.392 -
#> 12 US Lecce Parma AC 0.220 -
#> 13 AS Roma AC Milan 0.510 -
#> 14 Reggina Calcio AC Milan 0.318 -
#> 15 Hellas Verona AC Perugia 0.365 -
#> 16 Juventus AC Perugia 0.568 -
#> 17 ACF Fiorentina Atalanta 0.425 -
#> 18 Juventus Atalanta 0.527 -
#> 19 SSC Napoli AS Roma 0.230 -
#> 20 AS Bari AS Roma 0.255 -
#> 21 Lazio Roma Udinese Calcio 0.728 -
#> 22 Atalanta Udinese Calcio 0.555 -
#> 23 Bologna FC US Lecce 0.557 -
#> 24 Vicenza Calcio US Lecce 0.530 -
#> 25 AC Milan Brescia Calcio 0.510 -
#> 26 AC Perugia Brescia Calcio 0.390 -
#> 27 SSC Napoli Hellas Verona 0.488 -
#> 28 Parma AC Hellas Verona 0.640 -
#> 29 AC Perugia Reggina Calcio 0.525 -
#> 30 Atalanta Reggina Calcio 0.438 -
#> 31 Lazio Roma ACF Fiorentina 0.637 -
#> 32 AC Milan ACF Fiorentina 0.632 -
#> 33 Bologna FC Juventus 0.320 -
#> 34 Vicenza Calcio Juventus 0.275 -
#> 35 AS Bari Inter 0.288 -
#> 36 Parma AC Inter 0.620 -
#>
#> $round_plot
```
Red cells denote more likely home-wins (close to 0.6), such as: Lazio Roma - Fiorentina (observed result: 3-0, home win), Lazio Roma - Udinese (observed result: 3-1, home win), Juventus - AC Perugia (observed result: 1-0, home win), Brescia Calcio - AS Bari (observed result: 3-1, home win). Conversely, lighter cells denote more likely away wins (close to 0.6), such as: AS Bari - AS Roma (observed result: 1-4, away win), AS Bari - Inter (observed result: 1-2, away win).
## Rank-league reconstruction
Statisticians and football amateurs are much interested in the final rank-league predictions. However, predicting the final rank position (along with the teams' points) is often assimilated to an *oracle*, rather than a coherent statistical procedure.
We can provide here:
- **rank-league reconstruction** for the first model ```fit1_stan``` by using the **in-sample** replications $\mathcal{D}^{\text{rep}}$ (yellow ribbons for the credible intervals, solid blue lines for the observed cumulated points):
``` r
## Rank league reconstruction
# aggregated plot
foot_rank(object = fit1_stan, data = italy_2000)
```
``` r
# team-specific plot
foot_rank(
object = fit1_stan, data = italy_2000,
visualize = "individual"
)
```
- **rank-league prediction** (aggregated or at team-level) for the fourth model ```fit4_stan``` by using the **out-of-sample** replications $\tilde{\mathcal{D}}$ (yellow ribbons for the credible intervals, solid blue lines for the observed cumulated points):
``` r
## Rank predictions for individual teams
# aggregated plot
foot_rank(object = fit4_stan, data = italy_2000)
```
``` r
# team-specific plot
foot_rank(
object = fit4_stan, data = italy_2000,
teams = c("AC Milan", "AS Roma"),
visualize = "individual"
)
```
``` r
foot_rank(
object = fit4_stan, data = italy_2000,
visualize = "individual"
)
```
# Model comparisons
## Comparing Football Prediction Models using ```compare_foot```
Evaluating the performance of these models is crucial to understand their predictive power and reliability. The ```compare_foot``` function provides a comprehensive way to compare different models or probability matrices using metrics like accuracy, Brier score, ranked probability score (RPS), Pseudo $R^2$, and average coverage probability (ACP). It also offers the option to compute confusion matrices for a detailed performance analysis.
``` r
italy_2020_2021_fit <- italy %>%
dplyr::select(Season, home, visitor, hgoal, vgoal) %>%
dplyr::filter(Season == "2020" | Season == "2021") %>%
dplyr::mutate(periods = dplyr::case_when(
dplyr::row_number() <= 190 ~ 1,
dplyr::row_number() <= 380 ~ 2,
dplyr::row_number() <= 570 ~ 3,
TRUE ~ 4
)) %>% # Assign periods based on match number
dplyr::select(periods,
home_team = home,
away_team = visitor, home_goals = hgoal, away_goals = vgoal
)
fit_comp_1 <- stan_foot(
data = italy_2020_2021_fit,
model = "biv_pois",
home_effect = TRUE,
predict = 190,
dynamic_type = "season",
# parallel_chains = 4,
iter_sampling = n_iter
)
fit_comp_2 <- stan_foot(
data = italy_2020_2021_fit,
model = "double_pois",
home_effect = TRUE,
predict = 190,
dynamic_type = "season",
# parallel_chains = 4,
iter_sampling = n_iter
)
italy_2020_2021_test <- italy %>%
dplyr::select(Season, home, visitor, hgoal, vgoal) %>%
dplyr::filter(Season == "2014" | Season == "2015") %>%
dplyr::mutate(periods = dplyr::case_when(
dplyr::row_number() <= 190 ~ 1,
dplyr::row_number() <= 380 ~ 2,
dplyr::row_number() <= 570 ~ 3,
TRUE ~ 4
)) %>%
dplyr::filter(dplyr::row_number() > 570) %>%
dplyr::select(periods,
home_team = home,
away_team = visitor,
home_goals = hgoal,
away_goals = vgoal
)
```
``` r
compare_results_models <- compare_foot(
source = list(
biv_pois = fit_comp_1,
double_pois = fit_comp_2
),
test_data = italy_2020_2021_test,
metric = c("accuracy", "brier", "ACP", "pseudoR2", "RPS"),
conf_matrix = TRUE
)
print(compare_results_models, digits = 3)
#> Predictive Performance Metrics
#> Model RPS accuracy brier pseudoR2 ACP
#> biv_pois 0.255 0.353 0.716 0.311 0.350
#> double_pois 0.256 0.368 0.715 0.309 0.361
#>
#> Confusion Matrices
#> Model: biv_pois
#>
#> Actual
#> Predicted Home Win Draw Away Win
#> Home Win 51 28 32
#> Draw 0 0 0
#> Away Win 38 25 16
#>
#> Model: double_pois
#>
#> Actual
#> Predicted Home Win Draw Away Win
#> Home Win 57 33 35
#> Draw 0 0 0
#> Away Win 32 20 13
```
## LOOIC and WAIC
Comparing statistical models in terms of some **predictive information criteria** should conclude many analysis and can be carried out by using the Leave-one-out cross-validation criterion (**LOOIC**) and the Watanabe Akaike Information criterion (**WAIC**) performed by using the ```loo``` package. For more details about LOOIC and WAIC, read the paper @vehtari2016practical.
The general formulation for the predictive information criteria is the following:
$$ \text{crit}=-2 \widehat{\text{elpd}} = -2 (\widehat{\text{lpd}}- \text{parameters penalty})$$
- $\widehat{\text{elpd}}$: estimate of the expected log predictive density of the fitted model.
- $\widehat{\text{lpd}}$ is a measure of the log predictive density of the fitted model.
- parameters penalty is a penalization accounting for the effective number of parameters of the fitted model.
The interpretation is the following: the lower is the value for an information criterion, and the better is the estimated model's predictive accuracy. Moreover, if two competing models share the same value for the log predictive density, the model with less parameters is favored.
This is the **Occam's Razor** occurring in statistics:
> "Frustra fit per plura quod potest fieri per pauciora"
>
We can perform Bayesian model comparisons using the ``$loo()`` method, which computes an approximate LOO-CV using the ```loo``` package. To use this method, you must compute and store the pointwise log-likelihood in your Stan program. We will compare the static and weekly dynamic models for the 2000/2001 Italian Serie A:
``` r
### Model comparisons
## LOOIC, loo function
# compute loo
loo1 <- fit1_stan$fit$loo()
loo1_t <- fit1_stan_t$fit$loo()
loo2 <- fit2_stan$fit$loo()
loo3 <- fit3_stan$fit$loo()
loo3_t <- fit3_stan_t$fit$loo()
# compare three looic
loo_compare(loo1, loo1_t, loo2, loo3, loo3_t)
#> elpd_diff se_diff
#> model1 0.0 0.0
#> model2 -0.6 0.4
#> model3 -1.2 2.5
#> model4 -5.7 4.0
#> model5 -6.4 4.1
```
According to the above model LOOIC comparisons, the weekly-dynamic double Poisson models attain the lowest LOOIC values and are then the favored models in terms of predictive accuracy. The static model's ```fit1_stan``` final looic is suggesting that the assumption of static team-specific parameters is too restrictive and oversimplified to capture teams' skills over time and make reliable predictions. Anyway, from model checking we have the suggestion that even the static model has a reliable goodness of fit and could be used for some simplified analysis not requiring complex dynamic patterns.
# Extensions in next versions
Extensions and to-do list for the next package's versions:
1. **Data Web-scraping**: automatic routine to scrape data from internet;
2. **More numerical outputs**: posterior probabilities, betting strategies, etc.;
3. **Diagnostics, pp checks** designed for football;
4. **Teams' statistics**
5. **More covariates** to be included in the model (possibly by users).
6. **More priors choices**
# References
