Basic Insight
We look at some basic statistics. In our sample, nine countries have
participated each year.
participate <- colSums(is.finite(y))
names(participate)[participate == t]
#> [1] "CAN" "CHE" "CZE" "FIN" "LVA" "RUS" "SVK" "SWE" "USA"
The following countries hosted the championships, some of them
multiple times.
host <- sapply(x, FUN = sum)
host[host > 0L]
#> AUT BLR CAN CHE CZE DEU FIN FRA LVA NOR RUS SVK SWE
#> 1 1 1 2 2 3 3 1 1 1 3 2 3
In the years under analysis, the gold medals were awarded to the
following countries.
gold <- colSums(y == 1L)
gold[gold > 0L]
#> CAN CZE FIN RUS SVK SWE
#> 5 5 2 4 1 5
Model Estimation
The gasmodel package provides a single distribution on
rankings – the Plackett–Luce distribution.
distr(filter_type = "ranking")
#> distr_title param_title distr param type dim orthog default
#> 32 Plackett-Luce Worth pluce worth ranking multi FALSE TRUE
It is a convenient and simple probability distribution on rankings
utilizing a worth parameter for each item to be ranked. It originates
from Luce’s choice axiom and is also related to the Thurstone’s theory
of comparative judgment, see Luce (1977) and Yellott (1977). For more
details on this distribution, see Plackett (1975), Stern (1990), and
Critchlow et al. (1991).
We consider three different model specifications. We incorporate
x as an explanatory variable in our model to capture
possible home advantage. For each model specification, we assume a
panel-like structure where each worth parameter has its own intercept,
while the regression and dynamics parameters remain the same for all
worth parameters. In the gasmodel package, this structure
can be achieved using the coef_fix_value and
coef_fix_other arguments. Alternatively, for convenience,
the value panel_structure can be included in the
coef_fix_special argument. It is important to note that the
worth parameters in the Plackett–Luce distribution are not identifiable,
and it is common practice to impose a standardizing condition. In our
model, we enforce the condition that the sum of all \(\omega_i\) is 0. This can be accomplished
by including the value zero_sum_intercept in the
coef_fix_special argument.
First, we estimate the static model where there are no dynamics
involved. In this case, we set both the autoregressive and score orders
to zero. Either a single integer can be provided to determine the order
for all parameters, or a vector of integers can be supplied to specify
the order for individual parameters.
est_static <- gas(y = y, x = x, distr = "pluce", p = 0, q = 0,
coef_fix_special = c("zero_sum_intercept", "panel_structure"))
Second, we estimate the standard mean-reverting GAS model of order
one. In order to expedite the numerical optimization process, we
incorporate starting values based on the static model.
est_stnry <- gas(y = y, x = x, distr = "pluce",
coef_fix_special = c("zero_sum_intercept", "panel_structure"),
coef_start = as.vector(rbind(est_static$fit$par_unc / 2, 0, 0.5, 0.5)))
Third, we estimate the random walk model. In other words, we set the
autoregressive coefficient to 1. The easiest way to specify this is by
including the value random_walk in the
coef_fix_special argument. In our random walk model, we
consider the initial values of the worth parameters to be parameters to
be estimated. While the par_init argument does not directly
support this, we can set regress = "sep" and use cumulative
sums of exogenous variables to achieve this initialization for this
particular model. However, it is generally not recommended to estimate
initial parameter values as it introduces additional variables, lacks
reasonable asymptotics, and can lead to overfitting in finite samples.
It is important to approach the random walk model with caution, as it is
not stationary and the standard maximum likelihood asymptotics are not
valid.
est_walk <- gas(y = y, x = lapply(x, cumsum), distr = "pluce", regress = "sep",
coef_fix_special = c("zero_sum_intercept", "panel_structure", "random_walk"),
coef_start = as.vector(rbind(est_static$fit$par_unc, 0, 0.5, 1)))
To avoid redundancy, we will omit the output of the
gas() function, which contains rows for each coefficient of
each worth parameter. Since most coefficients are the same due to the
assumed panel structure, it is unnecessary to display them all. Instead,
we print only one set of the home advantage and dynamics
coefficients.
cbind(est_static = c("beta1" = unname(coef(est_static)[2]), "alpha1" = 0, "phi1" = 0),
est_stnry = coef(est_stnry)[2:4],
est_walk = coef(est_walk)[2:4])
#> est_static est_stnry est_walk
#> beta1 0.1707329 0.2274380 0.09873333
#> alpha1 0.0000000 0.3919431 0.34300134
#> phi1 0.0000000 0.5062478 1.00000000
In all three models, coefficient \(\alpha_1\) representing the home advantage
is positive but not significant.
cbind(est_static = c("beta1" = unname(est_static$fit$coef_pval)[2], "alpha1" = 0, "phi1" = 0),
est_stnry = est_stnry$fit$coef_pval[2:4],
est_walk = est_walk$fit$coef_pval[2:4])
#> est_static est_stnry est_walk
#> beta1 0.514887 3.772811e-01 5.995826e-01
#> alpha1 0.000000 2.141363e-06 2.634614e-09
#> phi1 0.000000 6.463390e-04 0.000000e+00
We compare the models using the Akaike information criterion (AIC).
The gas class allows for generic function
AIC(). In terms of AIC, the mean-reverting model
outperformed the remaining two by a wide margin.
AIC(est_static, est_stnry, est_walk)
#> df AIC
#> est_static 24 1299.600
#> est_stnry 26 1274.391
#> est_walk 25 1300.851
Time-Varying Worth Parameters
Additionally, we can examine the evolution of the worth parameters
for individual teams over the years. The point estimates of time-varying
parameter values can be directly obtained from the gas()
function. Using the generic plot() function allows us to
visualize the time-varying parameters of individual models. When
multiple parameters are time-varying, as in our scenario, the function
plots them in sequence. For the purpose of this document, we will only
display figures specific to the Canada team.
plot(est_static, which = 3)
plot(est_stnry, which = 3)
plot(est_walk, which = 3)
However, it is important to note that these estimates are subject to
uncertainty. To capture the uncertainty, we can utilize simulations by
leveraging the gas_filter() function, which accepts the
output of the gas() function as an argument. This allows us
to obtain the standard deviations and quantiles for the worth parameter
estimates, providing a more comprehensive understanding of the parameter
dynamics over time.
set.seed(42)
flt_stnry <- gas_filter(est_stnry)
To visualize time-varying parameters with confidence band, we can use
the plot() on the gas_filter object.
plot(flt_stnry, which = 3)
References
Alvo, M. and Yu, P. L. H. (2014). Statistical Methods for Ranking
Data. Springer. doi: 10.1007/978-1-4939-1471-5.
Critchlow, D. E., Fligner, M. A., and Verducci, J. S. (1991).
Probability Models on Rankings. Journal of Mathematical
Psychology, 35(3), 294–318. doi: 10.1016/0022-2496(91)90050-4.
Holý, V. and Zouhar, J. (2022). Modelling Time-Varying Rankings with
Autoregressive and Score-Driven Dynamics. Journal of the Royal
Statistical Society: Series C (Applied Statistics),
71(5), 1427–1450. doi: 10.1111/rssc.12584.
Luce, R. D. (1977). The Choice Axiom after Twenty Years. Journal
of Mathematical Psychology, 15(3), 215–233. doi:
10.1016/0022-2496(77)90032-3.
Plackett, R. L. (1975). The Analysis of Permutations. Journal of
the Royal Statistical Society: Series C (Applied Statistics),
24(2), 193–202. doi: 10.2307/2346567.
Stern, H. (1990). Models for Distributions on Permutations.
Journal of the American Statistical Association,
85(410), 558–564. doi: 10.1080/01621459.1990.10476235.
Yellott, J. I. (1977). The Relationship Between Luce’s Choice Axiom,
Thurstone’s Theory of Comparative Judgment, and the Double Exponential
Distribution. Journal of Mathematical Psychology,
15(2), 109–144. doi: 10.1016/0022-2496(77)90026-8.
Yu, P. L. H., Gu, J., and Xu, H. (2019). Analysis of Ranking Data.
Wiley Interdisciplinary Reviews: Computational Statistics,
11(6), e1483, 1–26. doi: 10.1002/wics.1483.
