Football_Data_Bayesian
2023-10-08
#install.packages(“devtools”)
#install_github("LeoEgidi/footBayes")
library(devtools)## Loading required package: usethislibrary(footBayes)
library(bayesplot)## This is bayesplot version 1.10.0## - Online documentation and vignettes at mc-stan.org/bayesplot## - bayesplot theme set to bayesplot::theme_default()## * Does _not_ affect other ggplot2 plots## * See ?bayesplot_theme_set for details on theme settinglibrary(loo)## This is loo version 2.6.0## - Online documentation and vignettes at mc-stan.org/loo## - As of v2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session.library(ggplot2)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.0
## ✔ readr 2.1.5## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsdata("italy")
italy <- as.data.frame(italy)
italy_2000 <- subset(italy[, c(2,3,4,6,7)], Season =="2000")
head(italy_2000)## Season home visitor hgoal vgoal
## 23421 2000 Udinese Calcio Brescia Calcio 4 2
## 23422 2000 AS Roma Bologna FC 2 0
## 23423 2000 AC Perugia US Lecce 1 1
## 23424 2000 Reggina Calcio Inter 2 1
## 23425 2000 SSC Napoli Juventus 1 2
## 23426 2000 AC Milan Vicenza Calcio 2 0# Fit Stan models
# no dynamics, no predictions
# 4 Markov chains, 'n_iter' iterations each
n_iter <- 200
fit1_stan <- stan_foot(data = italy_2000, model="biv_pois",chains = 4,iter = n_iter)## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.000217 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.17 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: WARNING: There aren't enough warmup iterations to fit the
## Chain 1: three stages of adaptation as currently configured.
## Chain 1: Reducing each adaptation stage to 15%/75%/10% of
## Chain 1: the given number of warmup iterations:
## Chain 1: init_buffer = 15
## Chain 1: adapt_window = 75
## Chain 1: term_buffer = 10
## Chain 1:
## Chain 1: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 1: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 1: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 1: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 1: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 1: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 1: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 1: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 1: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 1: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 1: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 1: Iteration: 200 / 200 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.328 seconds (Warm-up)
## Chain 1: 0.63 seconds (Sampling)
## Chain 1: 0.958 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0.000125 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.25 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: WARNING: There aren't enough warmup iterations to fit the
## Chain 2: three stages of adaptation as currently configured.
## Chain 2: Reducing each adaptation stage to 15%/75%/10% of
## Chain 2: the given number of warmup iterations:
## Chain 2: init_buffer = 15
## Chain 2: adapt_window = 75
## Chain 2: term_buffer = 10
## Chain 2:
## Chain 2: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 2: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 2: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 2: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 2: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 2: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 2: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 2: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 2: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 2: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 2: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 2: Iteration: 200 / 200 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.513 seconds (Warm-up)
## Chain 2: 0.497 seconds (Sampling)
## Chain 2: 1.01 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0.000122 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.22 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: WARNING: There aren't enough warmup iterations to fit the
## Chain 3: three stages of adaptation as currently configured.
## Chain 3: Reducing each adaptation stage to 15%/75%/10% of
## Chain 3: the given number of warmup iterations:
## Chain 3: init_buffer = 15
## Chain 3: adapt_window = 75
## Chain 3: term_buffer = 10
## Chain 3:
## Chain 3: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 3: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 3: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 3: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 3: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 3: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 3: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 3: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 3: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 3: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 3: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 3: Iteration: 200 / 200 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.506 seconds (Warm-up)
## Chain 3: 0.316 seconds (Sampling)
## Chain 3: 0.822 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0.000121 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.21 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: WARNING: There aren't enough warmup iterations to fit the
## Chain 4: three stages of adaptation as currently configured.
## Chain 4: Reducing each adaptation stage to 15%/75%/10% of
## Chain 4: the given number of warmup iterations:
## Chain 4: init_buffer = 15
## Chain 4: adapt_window = 75
## Chain 4: term_buffer = 10
## Chain 4:
## Chain 4: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 4: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 4: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 4: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 4: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 4: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 4: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 4: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 4: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 4: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 4: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 4: Iteration: 200 / 200 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.24 seconds (Warm-up)
## Chain 4: 0.348 seconds (Sampling)
## Chain 4: 0.588 seconds (Total)
## Chain 4:## Warning: There were 3 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: Examine the pairs() plot to diagnose sampling problems## Warning: The largest R-hat is 2.36, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#r-hat## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess# print of model summary for parameters:
# home, sigma_att, sigma_def
print(fit1_stan, pars =c("home", "rho", "sigma_att","sigma_def"))## Inference for Stan model: anon_model.
## 4 chains, each with iter=200; warmup=100; thin=1;
## post-warmup draws per chain=100, total post-warmup draws=400.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## home 0.28 0.00 0.07 0.13 0.24 0.28 0.33 0.41 389 1.01
## rho -1.71 0.02 0.32 -2.34 -1.87 -1.72 -1.47 -1.18 290 1.00
## sigma_att 0.21 0.01 0.07 0.10 0.16 0.20 0.25 0.38 63 1.08
## sigma_def 0.23 0.01 0.07 0.12 0.18 0.23 0.27 0.40 180 1.02
##
## Samples were drawn using NUTS(diag_e) at Sun May 19 23:51:56 2024.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).# Marginal posterior with bayesplot
posterior1 <- as.matrix(fit1_stan)
mcmc_areas(posterior1, regex_pars=c("home", "rho","sigma_att", "sigma_def"))
# Model's code extraction
fit1_stan@stanmodel## S4 class stanmodel 'anon_model' coded as follows:
## functions{
## real bipois_lpmf(int[] r , real mu1,real mu2,real mu3) {
## real ss;
## real log_s;
## real mus;
## int miny;
## miny = min(r[1], r[2]);
## ss = poisson_lpmf(r[1] | mu1) + poisson_lpmf(r[2] | mu2) -
## exp(mu3);
## if(miny > 0) {
## mus = -mu1-mu2+mu3;
## log_s = ss;
## for(k in 1:miny) {
## log_s = log_s + log(r[1] - k + 1) + mus
## + log(r[2] - k + 1)
## - log(k);
## ss = log_sum_exp(ss, log_s);
## }
## }
## return(ss);
## }
## }
## data{
## int N; // number of games
## int y[N,2];
## int nteams;
## int team1[N];
## int team2[N];
## real ranking[nteams];
## int<lower=0, upper=1> ind_home;
## // priors part
## int<lower=1,upper=4> prior_dist_num; // 1 gaussian, 2 t, 3 cauchy, 4 laplace
## int<lower=1,upper=4> prior_dist_sd_num; // 1 gaussian, 2 t, 3 cauchy, 4 laplace
## real hyper_df;
## real hyper_location;
## real hyper_sd_df;
## real hyper_sd_location;
## real hyper_sd_scale;
## }
## parameters{
## vector[nteams] att_raw;
## vector[nteams] def_raw;
## real<lower=0> sigma_att;
## real<lower=0> sigma_def;
## real beta;
## real rho;
## real home;
## real gamma;
## }
## transformed parameters{
## vector[nteams] att;
## vector[nteams] def;
## vector[3] theta[N];
## for (t in 1:nteams){
## att[t] = att_raw[t]-mean(att_raw);
## def[t] = def_raw[t]-mean(def_raw);
## }
## for (n in 1:N){
## theta[n,1] = exp(home*ind_home+att[team1[n]]+def[team2[n]]+
## (gamma/2)*(ranking[team1[n]]-ranking[team2[n]]));
## theta[n,2] = exp(att[team2[n]]+def[team1[n]]-
## (gamma/2)*(ranking[team1[n]]-ranking[team2[n]]));
## theta[n,3] = exp(rho);
## }
## }
## model{
## // log-priors for team-specific abilities
## for (t in 1:(nteams)){
## if (prior_dist_num == 1){
## target+= normal_lpdf(att_raw[t]|hyper_location, sigma_att);
## target+= normal_lpdf(def_raw[t]|hyper_location, sigma_def);
## }
## else if (prior_dist_num == 2){
## target+= student_t_lpdf(att_raw[t]|hyper_df, hyper_location, sigma_att);
## target+= student_t_lpdf(def_raw[t]|hyper_df, hyper_location, sigma_def);
## }
## else if (prior_dist_num == 3){
## target+= cauchy_lpdf(att_raw[t]|hyper_location, sigma_att);
## target+= cauchy_lpdf(def_raw[t]|hyper_location, sigma_def);
## }
## else if (prior_dist_num == 4){
## target+= double_exponential_lpdf(att_raw[t]|hyper_location, sigma_att);
## target+= double_exponential_lpdf(def_raw[t]|hyper_location, sigma_def);
## }
## }
## // log-hyperpriors for sd parameters
## if (prior_dist_sd_num == 1 ){
## target+=normal_lpdf(sigma_att|hyper_sd_location, hyper_sd_scale);
## target+=normal_lpdf(sigma_def|hyper_sd_location, hyper_sd_scale);
## }
## else if (prior_dist_sd_num == 2){
## target+=student_t_lpdf(sigma_att|hyper_sd_df, hyper_sd_location, hyper_sd_scale);
## target+=student_t_lpdf(sigma_def|hyper_sd_df, hyper_sd_location, hyper_sd_scale);
## }
## else if (prior_dist_sd_num == 3){
## target+=cauchy_lpdf(sigma_att|hyper_sd_location, hyper_sd_scale);
## target+=cauchy_lpdf(sigma_def|hyper_sd_location, hyper_sd_scale);
## }
## else if (prior_dist_sd_num == 4){
## target+=double_exponential_lpdf(sigma_att|hyper_sd_location, hyper_sd_scale);
## target+=double_exponential_lpdf(sigma_def|hyper_sd_location, hyper_sd_scale);
## }
## // log-priors fixed effects
## target+=normal_lpdf(rho|0,1);
## target+=normal_lpdf(home|0,5);
## target+=normal_lpdf(gamma|0,1);
## // likelihood
## for (n in 1:N){
## target+=poisson_lpmf(y[n,1]|theta[n,1]+theta[n,3]);
## target+=poisson_lpmf(y[n,2]|theta[n,2]+theta[n,3]);
## // target+=bipois_lpmf(y[n,]| theta[n,1],
## // theta[n,2], theta[n,3]);
## }
## }
## generated quantities{
## int y_rep[N,2];
## vector[N] log_lik;
## int diff_y_rep[N];
## //in-sample replications
## for (n in 1:N){
## y_rep[n,1] = poisson_rng(theta[n,1]+theta[n,3]);
## y_rep[n,2] = poisson_rng(theta[n,2]+theta[n,3]);
## diff_y_rep[n] = y_rep[n,1] - y_rep[n,2];
## log_lik[n] = poisson_lpmf(y[n,1]|theta[n,1]+theta[n,3])+
## poisson_lpmf(y[n,2]|theta[n,2]+theta[n,3]);
## //log_lik[n] =bipois_lpmf(y[n,]| theta[n,1],
## // theta[n,2], theta[n,3]);
## }
## }# Fitting MLE models
# no dynamics, no predictions
# Wald intervals
fit1_mle <- mle_foot(data = italy_2000,model="biv_pois",interval = "Wald")
fit1_mle$home## 2.5% mle 97.5%
## [1,] 0.2 0.3 0.39fit1_mle$sigma_att## NULLfit1_mle$sigma_def## NULL# Fitting 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 = student_t(4,0,NULL),
prior_sd = laplace(0,1),
iter = n_iter)##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.000179 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.79 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: WARNING: There aren't enough warmup iterations to fit the
## Chain 1: three stages of adaptation as currently configured.
## Chain 1: Reducing each adaptation stage to 15%/75%/10% of
## Chain 1: the given number of warmup iterations:
## Chain 1: init_buffer = 15
## Chain 1: adapt_window = 75
## Chain 1: term_buffer = 10
## Chain 1:
## Chain 1: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 1: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 1: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 1: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 1: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 1: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 1: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 1: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 1: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 1: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 1: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 1: Iteration: 200 / 200 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.555 seconds (Warm-up)
## Chain 1: 0.648 seconds (Sampling)
## Chain 1: 1.203 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0.000129 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.29 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: WARNING: There aren't enough warmup iterations to fit the
## Chain 2: three stages of adaptation as currently configured.
## Chain 2: Reducing each adaptation stage to 15%/75%/10% of
## Chain 2: the given number of warmup iterations:
## Chain 2: init_buffer = 15
## Chain 2: adapt_window = 75
## Chain 2: term_buffer = 10
## Chain 2:
## Chain 2: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 2: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 2: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 2: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 2: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 2: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 2: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 2: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 2: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 2: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 2: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 2: Iteration: 200 / 200 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.408 seconds (Warm-up)
## Chain 2: 0.615 seconds (Sampling)
## Chain 2: 1.023 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0.000123 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.23 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: WARNING: There aren't enough warmup iterations to fit the
## Chain 3: three stages of adaptation as currently configured.
## Chain 3: Reducing each adaptation stage to 15%/75%/10% of
## Chain 3: the given number of warmup iterations:
## Chain 3: init_buffer = 15
## Chain 3: adapt_window = 75
## Chain 3: term_buffer = 10
## Chain 3:
## Chain 3: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 3: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 3: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 3: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 3: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 3: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 3: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 3: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 3: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 3: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 3: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 3: Iteration: 200 / 200 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.628 seconds (Warm-up)
## Chain 3: 0.574 seconds (Sampling)
## Chain 3: 1.202 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0.000118 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.18 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: WARNING: There aren't enough warmup iterations to fit the
## Chain 4: three stages of adaptation as currently configured.
## Chain 4: Reducing each adaptation stage to 15%/75%/10% of
## Chain 4: the given number of warmup iterations:
## Chain 4: init_buffer = 15
## Chain 4: adapt_window = 75
## Chain 4: term_buffer = 10
## Chain 4:
## Chain 4: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 4: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 4: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 4: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 4: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 4: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 4: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 4: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 4: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 4: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 4: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 4: Iteration: 200 / 200 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.418 seconds (Warm-up)
## Chain 4: 0.568 seconds (Sampling)
## Chain 4: 0.986 seconds (Total)
## Chain 4:## Warning: There were 4 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: Examine the pairs() plot to diagnose sampling problems## Warning: The largest R-hat is 2.41, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#r-hat## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess# comparing posteriors
posterior1_t <- as.matrix(fit1_stan_t)
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))) ## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
#DYNAMIC FIT
# Fit Stan models
# seasonal dynamics, no predictions
# 4 Markov chains, 'n_iter' iterations each
fit2_stan <- stan_foot(data = italy_2000,
model="biv_pois",
dynamic_type ="weekly",
iter = n_iter) ## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.001772 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 17.72 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: WARNING: There aren't enough warmup iterations to fit the
## Chain 1: three stages of adaptation as currently configured.
## Chain 1: Reducing each adaptation stage to 15%/75%/10% of
## Chain 1: the given number of warmup iterations:
## Chain 1: init_buffer = 15
## Chain 1: adapt_window = 75
## Chain 1: term_buffer = 10
## Chain 1:
## Chain 1: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 1: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 1: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 1: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 1: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 1: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 1: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 1: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 1: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 1: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 1: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 1: Iteration: 200 / 200 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 33.155 seconds (Warm-up)
## Chain 1: 33.78 seconds (Sampling)
## Chain 1: 66.935 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0.001452 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 14.52 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: WARNING: There aren't enough warmup iterations to fit the
## Chain 2: three stages of adaptation as currently configured.
## Chain 2: Reducing each adaptation stage to 15%/75%/10% of
## Chain 2: the given number of warmup iterations:
## Chain 2: init_buffer = 15
## Chain 2: adapt_window = 75
## Chain 2: term_buffer = 10
## Chain 2:
## Chain 2: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 2: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 2: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 2: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 2: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 2: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 2: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 2: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 2: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 2: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 2: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 2: Iteration: 200 / 200 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 32.981 seconds (Warm-up)
## Chain 2: 46.999 seconds (Sampling)
## Chain 2: 79.98 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0.001432 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 14.32 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: WARNING: There aren't enough warmup iterations to fit the
## Chain 3: three stages of adaptation as currently configured.
## Chain 3: Reducing each adaptation stage to 15%/75%/10% of
## Chain 3: the given number of warmup iterations:
## Chain 3: init_buffer = 15
## Chain 3: adapt_window = 75
## Chain 3: term_buffer = 10
## Chain 3:
## Chain 3: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 3: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 3: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 3: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 3: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 3: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 3: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 3: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 3: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 3: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 3: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 3: Iteration: 200 / 200 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 36.318 seconds (Warm-up)
## Chain 3: 36.439 seconds (Sampling)
## Chain 3: 72.757 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0.001435 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 14.35 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: WARNING: There aren't enough warmup iterations to fit the
## Chain 4: three stages of adaptation as currently configured.
## Chain 4: Reducing each adaptation stage to 15%/75%/10% of
## Chain 4: the given number of warmup iterations:
## Chain 4: init_buffer = 15
## Chain 4: adapt_window = 75
## Chain 4: term_buffer = 10
## Chain 4:
## Chain 4: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 4: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 4: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 4: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 4: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 4: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 4: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 4: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 4: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 4: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 4: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 4: Iteration: 200 / 200 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 29.618 seconds (Warm-up)
## Chain 4: 35.691 seconds (Sampling)
## Chain 4: 65.309 seconds (Total)
## Chain 4:## Warning: There were 8 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
## https://mc-stan.org/misc/warnings.html#bfmi-low## Warning: Examine the pairs() plot to diagnose sampling problems## Warning: The largest R-hat is NA, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#r-hat## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-essprint(fit2_stan, pars =c("home", "rho", "sigma_att","sigma_def"))## Inference for Stan model: anon_model.
## 4 chains, each with iter=200; warmup=100; thin=1;
## post-warmup draws per chain=100, total post-warmup draws=400.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## home 0.25 0.01 0.07 0.10 0.20 0.25 0.30 0.39 60 1.07
## rho -1.73 0.02 0.31 -2.43 -1.94 -1.70 -1.48 -1.22 415 1.00
## sigma_att 0.06 0.01 0.02 0.04 0.05 0.06 0.08 0.11 7 2.61
## sigma_def 0.08 0.01 0.02 0.05 0.07 0.09 0.10 0.12 5 1.92
##
## Samples were drawn using NUTS(diag_e) at Sun May 19 23:58:08 2024.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).fit3_stan <- stan_foot(data = italy_2000,
model="double_pois",
dynamic_type = "weekly",
#cores = 4,
iter = n_iter) # double poisson## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.001808 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 18.08 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: WARNING: There aren't enough warmup iterations to fit the
## Chain 1: three stages of adaptation as currently configured.
## Chain 1: Reducing each adaptation stage to 15%/75%/10% of
## Chain 1: the given number of warmup iterations:
## Chain 1: init_buffer = 15
## Chain 1: adapt_window = 75
## Chain 1: term_buffer = 10
## Chain 1:
## Chain 1: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 1: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 1: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 1: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 1: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 1: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 1: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 1: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 1: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 1: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 1: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 1: Iteration: 200 / 200 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 32.04 seconds (Warm-up)
## Chain 1: 31.935 seconds (Sampling)
## Chain 1: 63.975 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0.001399 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 13.99 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: WARNING: There aren't enough warmup iterations to fit the
## Chain 2: three stages of adaptation as currently configured.
## Chain 2: Reducing each adaptation stage to 15%/75%/10% of
## Chain 2: the given number of warmup iterations:
## Chain 2: init_buffer = 15
## Chain 2: adapt_window = 75
## Chain 2: term_buffer = 10
## Chain 2:
## Chain 2: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 2: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 2: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 2: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 2: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 2: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 2: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 2: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 2: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 2: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 2: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 2: Iteration: 200 / 200 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 28.352 seconds (Warm-up)
## Chain 2: 44.15 seconds (Sampling)
## Chain 2: 72.502 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0.001363 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 13.63 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: WARNING: There aren't enough warmup iterations to fit the
## Chain 3: three stages of adaptation as currently configured.
## Chain 3: Reducing each adaptation stage to 15%/75%/10% of
## Chain 3: the given number of warmup iterations:
## Chain 3: init_buffer = 15
## Chain 3: adapt_window = 75
## Chain 3: term_buffer = 10
## Chain 3:
## Chain 3: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 3: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 3: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 3: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 3: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 3: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 3: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 3: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 3: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 3: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 3: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 3: Iteration: 200 / 200 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 28.482 seconds (Warm-up)
## Chain 3: 37.232 seconds (Sampling)
## Chain 3: 65.714 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0.001363 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 13.63 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: WARNING: There aren't enough warmup iterations to fit the
## Chain 4: three stages of adaptation as currently configured.
## Chain 4: Reducing each adaptation stage to 15%/75%/10% of
## Chain 4: the given number of warmup iterations:
## Chain 4: init_buffer = 15
## Chain 4: adapt_window = 75
## Chain 4: term_buffer = 10
## Chain 4:
## Chain 4: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 4: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 4: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 4: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 4: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 4: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 4: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 4: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 4: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 4: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 4: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 4: Iteration: 200 / 200 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 24.905 seconds (Warm-up)
## Chain 4: 43.815 seconds (Sampling)
## Chain 4: 68.72 seconds (Total)
## Chain 4:## Warning: There were 4 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
## https://mc-stan.org/misc/warnings.html#bfmi-low## Warning: Examine the pairs() plot to diagnose sampling problems## Warning: The largest R-hat is NA, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#r-hat## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-essprint(fit3_stan, pars =c("home", "sigma_att",
"sigma_def"))## Inference for Stan model: anon_model.
## 4 chains, each with iter=200; warmup=100; thin=1;
## post-warmup draws per chain=100, total post-warmup draws=400.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## home 0.40 0.00 0.05 0.30 0.37 0.40 0.43 0.49 299 1.02
## sigma_att 0.06 0.01 0.02 0.04 0.04 0.06 0.08 0.09 3 2.50
## sigma_def 0.07 0.01 0.01 0.05 0.06 0.07 0.08 0.10 6 1.78
##
## Samples were drawn using NUTS(diag_e) at Mon May 20 00:03:29 2024.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).# Fit Stan models
# weekly dynamics, no predictions
# 4 chains, 'n_iter' iterations each
fit3_stan_t <- stan_foot(data = italy_2000,
model="double_pois",
prior = student_t(4,0, NULL), # 4 df
prior_sd = cauchy(0,25),
dynamic_type = "weekly",
#cores = 4,
iter = n_iter) # double poisson##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.001506 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 15.06 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: WARNING: There aren't enough warmup iterations to fit the
## Chain 1: three stages of adaptation as currently configured.
## Chain 1: Reducing each adaptation stage to 15%/75%/10% of
## Chain 1: the given number of warmup iterations:
## Chain 1: init_buffer = 15
## Chain 1: adapt_window = 75
## Chain 1: term_buffer = 10
## Chain 1:
## Chain 1: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 1: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 1: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 1: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 1: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 1: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 1: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 1: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 1: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 1: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 1: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 1: Iteration: 200 / 200 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 27.235 seconds (Warm-up)
## Chain 1: 34.593 seconds (Sampling)
## Chain 1: 61.828 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0.001406 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 14.06 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: WARNING: There aren't enough warmup iterations to fit the
## Chain 2: three stages of adaptation as currently configured.
## Chain 2: Reducing each adaptation stage to 15%/75%/10% of
## Chain 2: the given number of warmup iterations:
## Chain 2: init_buffer = 15
## Chain 2: adapt_window = 75
## Chain 2: term_buffer = 10
## Chain 2:
## Chain 2: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 2: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 2: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 2: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 2: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 2: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 2: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 2: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 2: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 2: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 2: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 2: Iteration: 200 / 200 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 28.863 seconds (Warm-up)
## Chain 2: 34.23 seconds (Sampling)
## Chain 2: 63.093 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0.001419 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 14.19 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: WARNING: There aren't enough warmup iterations to fit the
## Chain 3: three stages of adaptation as currently configured.
## Chain 3: Reducing each adaptation stage to 15%/75%/10% of
## Chain 3: the given number of warmup iterations:
## Chain 3: init_buffer = 15
## Chain 3: adapt_window = 75
## Chain 3: term_buffer = 10
## Chain 3:
## Chain 3: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 3: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 3: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 3: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 3: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 3: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 3: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 3: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 3: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 3: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 3: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 3: Iteration: 200 / 200 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 25.124 seconds (Warm-up)
## Chain 3: 33.65 seconds (Sampling)
## Chain 3: 58.774 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0.001406 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 14.06 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: WARNING: There aren't enough warmup iterations to fit the
## Chain 4: three stages of adaptation as currently configured.
## Chain 4: Reducing each adaptation stage to 15%/75%/10% of
## Chain 4: the given number of warmup iterations:
## Chain 4: init_buffer = 15
## Chain 4: adapt_window = 75
## Chain 4: term_buffer = 10
## Chain 4:
## Chain 4: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 4: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 4: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 4: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 4: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 4: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 4: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 4: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 4: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 4: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 4: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 4: Iteration: 200 / 200 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 24.557 seconds (Warm-up)
## Chain 4: 32.239 seconds (Sampling)
## Chain 4: 56.796 seconds (Total)
## Chain 4:## Warning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
## https://mc-stan.org/misc/warnings.html#bfmi-low## Warning: Examine the pairs() plot to diagnose sampling problems## Warning: The largest R-hat is NA, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#r-hat## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-essprint(fit3_stan_t, pars =c("home", "sigma_att",
"sigma_def"))## Inference for Stan model: anon_model.
## 4 chains, each with iter=200; warmup=100; thin=1;
## post-warmup draws per chain=100, total post-warmup draws=400.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## home 0.39 0.00 0.05 0.29 0.36 0.39 0.43 0.48 200 1.01
## sigma_att 0.06 0.01 0.02 0.03 0.04 0.05 0.07 0.09 3 1.79
## sigma_def 0.06 0.00 0.01 0.03 0.05 0.06 0.06 0.09 14 1.37
##
## Samples were drawn using NUTS(diag_e) at Mon May 20 00:07:41 2024.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).# Plotting abilities: credible and confidence 95% intervals
foot_abilities(object = fit1_stan, data = italy_2000, cex.var = 1)
foot_abilities(object = fit1_mle, data = italy_2000, cex.var = 1)
# Plotting abilities: credible and confidence 95% intervals
foot_abilities(fit2_stan, italy_2000)
# PP checks: aggregated goal's differences and ordered goal differences
pp_foot(data = italy_2000, object = fit1_stan,
type = "aggregated")## Warning: Continuous limits supplied to discrete scale.
## ℹ Did you mean `limits = factor(...)` or `scale_*_continuous()`?## $pp_plot
##
## $pp_table
## goal diff. Bayesian p-value
## 1 -3 0.338
## 2 -2 0.870
## 3 -1 0.980
## 4 0 0.040
## 5 1 0.208
## 6 2 0.080
## 7 3 0.895pp_foot(data = italy_2000, object = fit1_stan,
type = "matches")## $pp_plot
##
## $pp_table
## 1-alpha emp. coverage
## 1 0.95 0.951# PPC densities overlay with the bayesplot package
# extracting the replications
sims <-rstan::extract(fit1_stan)
goal_diff <- italy_2000$hgoal-italy_2000$vgoal
# plotting data density vs replications densities
ppc_dens_overlay(goal_diff, sims$y_rep[,,1]-sims$y_rep[,,2], bw = 0.5)
# Fit Stan models
## weekly dynamics, predictions of last four weeks
## 4 chains 'n_iter' iterations each
fit4_stan <- stan_foot(data = italy_2000,
model="biv_pois",
predict = 36,
dynamic_type = "weekly",
iter = n_iter) ## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 15.0.0 (clang-1500.3.9.4)’
## using SDK: ‘MacOSX14.4.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.001858 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 18.58 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: WARNING: There aren't enough warmup iterations to fit the
## Chain 1: three stages of adaptation as currently configured.
## Chain 1: Reducing each adaptation stage to 15%/75%/10% of
## Chain 1: the given number of warmup iterations:
## Chain 1: init_buffer = 15
## Chain 1: adapt_window = 75
## Chain 1: term_buffer = 10
## Chain 1:
## Chain 1: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 1: Exception: poisson_rng: Rate parameter is 2.21409e+11, but must be less than 1.07374e+09 (in 'anon_model', line 170, column 8 to column 73)
## Chain 1: Exception: poisson_rng: Rate parameter is 1.91258e+10, but must be less than 1.07374e+09 (in 'anon_model', line 170, column 8 to column 73)
## Chain 1: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 1: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 1: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 1: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 1: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 1: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 1: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 1: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 1: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 1: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 1: Iteration: 200 / 200 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 42.979 seconds (Warm-up)
## Chain 1: 38.392 seconds (Sampling)
## Chain 1: 81.371 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0.001403 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 14.03 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: WARNING: There aren't enough warmup iterations to fit the
## Chain 2: three stages of adaptation as currently configured.
## Chain 2: Reducing each adaptation stage to 15%/75%/10% of
## Chain 2: the given number of warmup iterations:
## Chain 2: init_buffer = 15
## Chain 2: adapt_window = 75
## Chain 2: term_buffer = 10
## Chain 2:
## Chain 2: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 2: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 2: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 2: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 2: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 2: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 2: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 2: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 2: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 2: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 2: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 2: Iteration: 200 / 200 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 35.529 seconds (Warm-up)
## Chain 2: 47.955 seconds (Sampling)
## Chain 2: 83.484 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0.001431 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 14.31 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: WARNING: There aren't enough warmup iterations to fit the
## Chain 3: three stages of adaptation as currently configured.
## Chain 3: Reducing each adaptation stage to 15%/75%/10% of
## Chain 3: the given number of warmup iterations:
## Chain 3: init_buffer = 15
## Chain 3: adapt_window = 75
## Chain 3: term_buffer = 10
## Chain 3:
## Chain 3: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 3: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 3: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 3: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 3: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 3: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 3: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 3: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 3: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 3: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 3: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 3: Iteration: 200 / 200 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 39.777 seconds (Warm-up)
## Chain 3: 36.245 seconds (Sampling)
## Chain 3: 76.022 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0.001424 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 14.24 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: WARNING: There aren't enough warmup iterations to fit the
## Chain 4: three stages of adaptation as currently configured.
## Chain 4: Reducing each adaptation stage to 15%/75%/10% of
## Chain 4: the given number of warmup iterations:
## Chain 4: init_buffer = 15
## Chain 4: adapt_window = 75
## Chain 4: term_buffer = 10
## Chain 4:
## Chain 4: Iteration: 1 / 200 [ 0%] (Warmup)
## Chain 4: Iteration: 20 / 200 [ 10%] (Warmup)
## Chain 4: Iteration: 40 / 200 [ 20%] (Warmup)
## Chain 4: Iteration: 60 / 200 [ 30%] (Warmup)
## Chain 4: Iteration: 80 / 200 [ 40%] (Warmup)
## Chain 4: Iteration: 100 / 200 [ 50%] (Warmup)
## Chain 4: Iteration: 101 / 200 [ 50%] (Sampling)
## Chain 4: Iteration: 120 / 200 [ 60%] (Sampling)
## Chain 4: Iteration: 140 / 200 [ 70%] (Sampling)
## Chain 4: Iteration: 160 / 200 [ 80%] (Sampling)
## Chain 4: Iteration: 180 / 200 [ 90%] (Sampling)
## Chain 4: Iteration: 200 / 200 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 37.825 seconds (Warm-up)
## Chain 4: 34.637 seconds (Sampling)
## Chain 4: 72.462 seconds (Total)
## Chain 4:## Warning in validityMethod(object): The following variables have undefined
## values: y_rep[1,1],The following variables have undefined values:
## y_rep[2,1],The following variables have undefined values: y_rep[3,1],The
## following variables have undefined values: y_rep[4,1],The following variables
## have undefined values: y_rep[5,1],The following variables have undefined
## values: y_rep[6,1],The following variables have undefined values:
## y_rep[7,1],The following variables have undefined values: y_rep[8,1],The
## following variables have undefined values: y_rep[9,1],The following variables
## have undefined values: y_rep[10,1],The following variables have undefined
## values: y_rep[11,1],The following variables have undefined values:
## y_rep[12,1],The following variables have undefined values: y_rep[13,1],The
## following variables have undefined values: y_rep[14,1],The following variables
## have undefined values: y_rep[15,1],The following variables have undefined
## values: y_rep[16,1],The following variables have undefined values:
## y_rep[17,1],The following variables have undefined values: y_rep[18,1],The
## following variables have undefined values: y_rep[19,1],The following variables
## have undefined values: y_rep[20,1],The following variables have undefined
## values: y_rep[21,1],The following variables have undefined values:
## y_rep[22,1],The following variables have undefined values: y_rep[23,1],The
## following variables have undefined values: y_rep[24,1],The following variables
## have undefined values: y_rep[25,1],The following variables have undefined
## values: y_rep[26,1],The following variables have undefined values:
## y_rep[27,1],The following variables have undefined values: y_rep[28,1],The
## following variables have undefined values: y_rep[29,1],The following variables
## have undefined values: y_rep[30,1],The following variables have undefined
## values: y_rep[31,1],The following variables have undefined values:
## y_rep[32,1],The following variables have undefined values: y_rep[33,1],The
## following variables have undefined values: y_rep[34,1],The following variables
## have undefined values: y_rep[35,1],The following variables have undefined
## values: y_rep[36,1],The following variables have undefined values:
## y_rep[37,1],The following variables have undefined values: y_rep[38,1],The
## following variables have undefined values: y_rep[39,1],The following variables
## have undefined values: y_rep[40,1],The following variables have undefined
## values: y_rep[41,1],The following variables have undefined values:
## y_rep[42,1],The following variables have undefined values: y_rep[43,1],The
## following variables have undefined values: y_rep[44,1],The following variables
## have undefined values: y_rep[45,1],The following variables have undefined
## values: y_rep[46,1],The following variables have undefined values:
## y_rep[47,1],The following variables have undefined values: y_rep[48,1],The
## following variables have undefined values: y_rep[49,1],The following variables
## have undefined values: y_rep[50,1],The following variables have undefined
## values: y_rep[51,1],The following variables have undefined values:
## y_rep[52,1],The following variables have undefined values: y_rep[53,1],The
## following variables have undefined values: y_rep[54,1],The following variables
## have undefined values: y_rep[55,1],The following variables have undefined
## values: y_rep[56,1],The following variables have undefined values:
## y_rep[57,1],The following variables have undefined values: y_rep[58,1],The
## following variables have undefined values: y_rep[59,1],The following variables
## have undefined values: y_rep[60,1],The following variables have undefined
## values: y_rep[61,1],The following variables have undefined values:
## y_rep[62,1],The following variables have undefined values: y_rep[63,1],The
## following variables have undefined values: y_rep[64,1],The following variables
## have undefined values: y_rep[65,1],The following variables have undefined
## values: y_rep[66,1],The following variables have undefined values:
## y_rep[67,1],The following variables have undefined values: y_rep[68,1],The
## following variables have undefined values: y_rep[69,1],The following variables
## have undefined values: y_rep[70,1],The following variables have undefined
## values: y_rep[71,1],The following variables have undefined values:
## y_rep[72,1],The following variables have undefined values: y_rep[73,1],The
## following variables have undefined values: y_rep[74,1],The following variables
## have undefined values: y_rep[75,1],The following variables have undefined
## values: y_rep[76,1],The following variables have undefined values:
## y_rep[77,1],The following variables have undefined values: y_rep[78,1],The
## following variables have undefined values: y_rep[79,1],The following variables
## have undefined values: y_rep[80,1],The following variables have undefined
## values: y_rep[81,1],The following variables have undefined values:
## y_rep[82,1],The following variables have undefined values: y_rep[83,1],The
## following variables have undefined values: y_rep[84,1],The following variables
## have undefined values: y_rep[85,1],The following variables have undefined
## values: y_rep[86,1],The following variables have undefined values:
## y_rep[87,1],The following variables have undefined values: y_rep[88,1],The
## following variables have undefined values: y_rep[89,1],The following variables
## have undefined values: y_rep[90,1],The following variables have undefined
## values: y_rep[91,1],The following variables have undefined values:
## y_rep[92,1],The following variables have undefined values: y_rep[93,1],The
## following variables have undefined values: y_rep[94,1],The following variables
## have undefined values: y_rep[95,1],The following variables have undefined
## values: y_rep[96,1],The following variables have undefined values:
## y_rep[97,1],The following variables have undefined values: y_rep[98,1],The
## following variables have undefined values: y_rep[99,1],The following variables
## have undefined values: y_rep[100,1],The following variables have undefined
## values: y_rep[101,1],The following variables have undefined values:
## y_rep[102,1],The following variables have undefined values: y_rep[103,1],The
## following variables have undefined values: y_rep[104,1],The following variables
## have undefined values: y_rep[105,1],The following variables have undefined
## values: y_rep[106,1],The following variables have undefined values:
## y_rep[107,1],The following variables have undefined values: y_rep[108,1],The
## following variables have undefined values: y_rep[109,1],The following variables
## have undefined values: y_rep[110,1],The following variables have undefined
## values: y_rep[111,1],The following variables have undefined values:
## y_rep[112,1],The following variables have undefined values: y_rep[113,1],The
## following variables have undefined values: y_rep[114,1],The following variables
## have undefined values: y_rep[115,1],The following variables have undefined
## values: y_rep[116,1],The following variables have undefined values:
## y_rep[117,1],The following variables have undefined values: y_rep[118,1],The
## following variables have undefined values: y_rep[119,1],The following variables
## have undefined values: y_rep[120,1],The following variables have undefined
## values: y_rep[121,1],The following variables have undefined values:
## y_rep[122,1],The following variables have undefined values: y_rep[123,1],The
## following variables have undefined values: y_rep[124,1],The following variables
## have undefined values: y_rep[125,1],The following variables have undefined
## values: y_rep[126,1],The following variables have undefined values:
## y_rep[127,1],The following variables have undefined values: y_rep[128,1],The
## following variables have undefined values: y_rep[129,1],The following variables
## have undefined values: y_rep[130,1],The following variables have undefined
## values: y_rep[131,1],The following variables have undefined values:
## y_rep[132,1],The following variables have undefined values: y_rep[133,1],The
## following variables have undefined values: y_rep[134,1],The following variables
## have undefined values: y_rep[135,1],The following variables have undefined
## values: y_rep[136,1],Th## Warning: There were 12 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
## https://mc-stan.org/misc/warnings.html#bfmi-low## Warning: Examine the pairs() plot to diagnose sampling problems## Warning: The largest R-hat is NA, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#r-hat## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-essfoot_prob(object = fit4_stan, data = italy_2000,
home_team = "Reggina Calcio",
away_team= "AC Milan")## Warning in rbind(deparse.level, ...): number of columns of result, 6, is not a
## multiple of vector length 4 of arg 1## $prob_table
## home_team away_team prob_h prob_d prob_a mlo
## 1 Reggina Calcio AC Milan 0.291 0.242 0.467 1-1 (0.109)
##
## $prob_plot
## Home win out-of-sample probabilities
foot_round_robin(data = italy_2000, object = fit4_stan)## $round_plot
##
## $round_table
## # A tibble: 36 × 4
## Home Away Home_prob Observed
## <chr> <chr> <chr> <chr>
## 1 Hellas Verona Bologna FC 0.338 -
## 2 Inter Bologna FC 0.418 -
## 3 Udinese Calcio SSC Napoli 0.527 -
## 4 ACF Fiorentina SSC Napoli 0.545 -
## 5 US Lecce Lazio Roma 0.2 -
## 6 Inter Lazio Roma 0.275 -
## 7 Brescia Calcio AS Bari 0.682 -
## 8 Reggina Calcio AS Bari 0.63 -
## 9 Udinese Calcio Vicenza Calcio 0.402 -
## 10 Brescia Calcio Vicenza Calcio 0.488 -
## # ℹ 26 more rows## Rank league reconstruction
# aggregated plot
foot_rank(data = italy_2000, object = fit1_stan)## $rank_table
## teams mid obs lo hi
## 1 AS Roma 61 75 55 66
## 2 Juventus 60 73 54 66
## 3 Lazio Roma 58 69 52 64
## 4 Parma AC 54 56 49 60
## 5 AC Milan 50 49 44 58
## 6 Atalanta 48 44 41 54
## 7 Brescia Calcio 48 44 42 54
## 8 Inter 47 51 40 53
## 9 ACF Fiorentina 47 43 41 52
## 10 AC Perugia 45 42 39 51
## 11 Bologna FC 45 43 40 51
## 12 Udinese Calcio 42 38 36 49
## 13 SSC Napoli 41 36 35 46
## 14 Vicenza Calcio 41 36 35 47
## 15 US Lecce 41 37 35 47
## 16 Reggina Calcio 40 37 34 46
## 17 Hellas Verona 39 37 34 45
## 18 AS Bari 33 20 28 39
##
## $rank_plot
# team-specific plot
foot_rank(data = italy_2000, object = fit1_stan, visualize = "individual")## $rank_plot
## Rank predictions for individual teams
# aggregated plot
foot_rank(data = italy_2000, object = fit4_stan)## $rank_table
## teams mid obs lo hi
## 1 AS Roma 74 75 72 76
## 2 Lazio Roma 71 69 69 72
## 3 Juventus 68 73 67 70
## 4 Parma AC 58 56 56 60
## 5 AC Milan 54 49 52 56
## 6 Bologna FC 48 43 46 49
## 7 Atalanta 48 44 46 49
## 8 Inter 48 51 47 51
## 9 AC Perugia 45 42 43 47
## 10 Brescia Calcio 43 44 40 45
## 11 ACF Fiorentina 43 43 42 45
## 12 Udinese Calcio 38 38 37 40
## 13 Vicenza Calcio 37 36 35 39
## 14 US Lecce 36 37 34 37
## 15 Reggina Calcio 34 37 32 36
## 16 SSC Napoli 32 36 30 34
## 17 Hellas Verona 32 37 31 34
## 18 AS Bari 23 20 21 25
##
## $rank_plot
# team-specific plot
foot_rank(italy_2000, fit4_stan,
team_sel = c("AC Milan", "AS Roma"),
visualize = "individual")## $rank_plot
foot_rank(italy_2000, fit4_stan,
visualize = "individual")## $rank_plot
### Model comparisons
## LOOIC, loo function
# extract pointwise log-likelihood
log_lik_1 <- extract_log_lik(fit1_stan)
log_lik_1_t <- extract_log_lik(fit1_stan_t)
log_lik_2 <- extract_log_lik(fit2_stan)
log_lik_3_t <- extract_log_lik(fit3_stan_t)
# compute loo
loo1 <- loo(log_lik_1)## Warning: Relative effective sample sizes ('r_eff' argument) not specified.
## For models fit with MCMC, the reported PSIS effective sample sizes and
## MCSE estimates will be over-optimistic.## Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.loo1_t <- loo(log_lik_1_t)## Warning: Relative effective sample sizes ('r_eff' argument) not specified.
## For models fit with MCMC, the reported PSIS effective sample sizes and
## MCSE estimates will be over-optimistic.
## Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.loo2 <- loo(log_lik_2)## Warning: Relative effective sample sizes ('r_eff' argument) not specified.
## For models fit with MCMC, the reported PSIS effective sample sizes and
## MCSE estimates will be over-optimistic.
## Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.loo3_t <- loo(log_lik_3_t)## Warning: Relative effective sample sizes ('r_eff' argument) not specified.
## For models fit with MCMC, the reported PSIS effective sample sizes and
## MCSE estimates will be over-optimistic.## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.# compare three looic
loo_compare(loo1, loo1_t, loo2, loo3_t)## elpd_diff se_diff
## model1 0.0 0.0
## model2 -1.1 0.5
## model3 -3.9 2.8
## model4 -7.3 4.1