Athlete rating in multi-competitor games with scored outcomes via monotone transformations

A standard DLM for rating athletes models each athlete’s observed scores as normally distributed around their latent ability, which evolves over time. The model likelihood for the score observed from a single athlete competing during time $t$ takes the form:

where $y_{t}$ is the observed score, $\theta_{t}$ is the athlete’s latent ability parameter at time $t$, and $\sigma^{2}$ is the observation variance. In this paper, we allow the latent ability parameter to evolve between time points as a normal random walk:

where $w$ is an additional parameter that controls how much latent abilities may vary over time relative to the observation variance. Other stochastic processes may also be considered for the evolution of latent ability parameters, such as a mean-preserving random walk (Glickman and Hennessy, 2015) or an autoregressive process (Glickman and Stern, 1998). In this paper, we use a simple normal random walk rather than an autoregressive process to acknowledge that weak athletes’ latent abilities are not likely to improve, i.e., revert to the mean, when they avoid playing games over several time periods. Unlike these other stochastic processes, the normal random walk diverges over time; this increasing variance represents the increasing uncertainty about each athlete’s latent ability if they do not play any games. In practice, we place an upper limit on the variance of athlete latent abilities to prevent them from growing too large (see Appendix A.1 for details).

In a general sporting setting, we divide time into discrete rating periods, with multiple games within each rating period and multiple athletes within each game. The choice of rating period should reflect the structure of the data. Our model assumes that athletes’ latent abilities remain constant within each rating period, but that they may vary between rating periods. Sports typically have regularly spaced seasons which serve as natural rating period breakpoints; for example, the biathlon season spans November through March, so we may assign the full season to a single rating period or split each season into two rating periods to allow mid-season changes in ability. Shorter rating periods generally reduce prediction bias, since more recent data is used for prediction, but they increase variance, since fewer games are used to infer latent abilities within each period. Appendix A.2 provides further guidance on this bias-variance tradeoff and illustrates how different choices of rating periods marginally affect results on real data.

Let $p$ denote the total number of athletes in the data. Let $T$ denote the total number of discrete rating periods, indexed by $t=1,\dots,T$. Finally, let $n_{t}$ denote the total number of observed scores within rating period $t$. Athletes may compete in any number (including zero) of games within each rating period. The data model (Equation 1) and innovation model (Equation 2) may now be rewritten in multivariate form as:

Here, $\mathbf{y}_{t}$ is the $n_{t}\times 1$ column vector of observed scores, ${\bm{\theta}}_{t}$ is the $p\times 1$ column vector of athlete latent abilities for rating period $t$, and $I_{k}$ denotes the $k\times k$ identity matrix. The $n_{t}\times p$ model matrix $X_{t}$ simply matches each athlete’s observed score(s) to their latent ability parameter.²²2Temporarily suppressing the time subscript $t$, the matrix $X$ has a single nonzero entry per row. Per row $r$, if entry $r$ of the column vector $\mathbf{y}$ corresponds to a score earned by athlete $a$, then $X_{ra}$ is set to equal 1.

In practice, athlete-rating DLMs need to account for conditions that affect game scores in a manner unrelated to latent athlete abilities. For example, a hot day might make all athletes in a race run more slowly, but the increased race times do not indicate weaker athletes. One approach to incorporate game-specific variation is to assume game-specific intercepts as part of the outcome model. This approach has been used, for example, by Glickman and Stern (1998).

Instead of assuming game-specific intercepts, we pre-process the data by subtracting game-specific means from the observed scores. This approach avoids concerns relating to the arbitrary specification of priors for the intercepts, which could affect downstream results. Each game-centered score $\tilde{y}_{tg}$ for an individual athlete in game $g$ within rating period $t$ may then be modeled using a game-centered latent ability, as:

The value $\bar{\theta}_{tg}$ denotes the average latent ability across all of the players in game $g$ within rating period $t$. Subtracting $\bar{\theta}_{tg}$ from each athlete’s latent ability adjusts for the fact that competing against a stronger pool of opponents in a game naturally results in worse scores relative to the competition.

The variation in scores may also differ across games within a sport. For example, biathlon races vary in length. Longer races give strong athletes more time to build a larger lead and weak athletes more time to fall behind, increasing the observed variation in scores. Ideally, the transformation introduced in the following section would address this issue by shrinking extreme game-centered scores toward the mean. If the variation in scores significantly differs from game to game in a particular application, however, it may be reasonable to additionally pre-process the scores e.g., by log-transforming them, prior to game-centering, or (for large multicompetitor sports) by game-centering and then scaling each game’s scores to have variance 1.

To simplify notation for vector of pre-processed, game-centered scores $\tilde{\mathbf{y}}_{t}$, we write:

where $\tilde{\mathbf{y}}_{t}\equiv\begin{bmatrix}\tilde{\mathbf{y}}_{t1}\\ \vdots\\ \tilde{\mathbf{y}}_{tq}\end{bmatrix}$ and $\bar{X}_{t}\equiv\begin{bmatrix}H_{n_{t1}}X_{t1}\\ \vdots\\ H_{n_{tg_{t}}}X_{tg_{t}}\end{bmatrix}$, for centering matrices $H_{k}\equiv I_{k}-\mathbf{1}_{k}\mathbf{1}_{k}^{T}$ with $\mathbf{1}_{k}$ denoting the $k\times 1$ column vector with entries equal to one. For games $g=1,\dots,g_{t}$ within rating period $t$, the $n_{g}\times 1$ vector of scores and $n_{g}\times p$ model matrix for game $g$ are written as $\tilde{\mathbf{y}}_{tg}$ and $X_{tg}$, respectively.

In many games, we might also suspect that athletes’ game-centered scores are not normally distributed around their game-centered latent abilities as assumed by Equation 3. Instead, we may assume that some transformation of the athletes’ scores is normally distributed around their latent abilities, so that the resulting model for transformed outcomes is:

The transformation $\tau_{\bm{\lambda}}(\cdot)$ is a function parameterized by a vector-valued parameter ${\bm{\lambda}}$. While the above assumption is likely reasonable for most sports, we note that it may be difficult to justify for sports with very low scores, e.g., soccer or hockey.

In this paper, we use the monotone spline transformation from Ramsay (1988), but any monotone transformation with a computable Jacobian would work as well. The monotone spline transformation is a polynomial spline built from an I-spline basis (see Ramsay (1988) for more detail). For a given polynomial order $d$ and knot sequence $\mathbf{k}$, an I-spline basis consists of $B$ fixed, monotonically increasing basis functions $I_{b}(y\mid d,\mathbf{k})$, $b=1,\ldots,B$. The monotone spline transformation is then constructed as a linear combination of these basis functions:

where $\lambda_{0}$ represents an intercept term and $\lambda_{1},\dots,\lambda_{B}$ determine the shape of the transformation. We treat $\lambda_{0}$ as fixed, and define the transformation parameter ${\bm{\lambda}}$ as ${\bm{\lambda}}\equiv\begin{bmatrix}\lambda_{1}\ \dots\ \lambda_{B}\end{bmatrix}^{T}$. Note that because each basis function $I_{b}(\cdot)$ is monotone increasing, constraining the parameters $\lambda_{1},\dots,\lambda_{B}$ to be non-negative ensures that the resulting spline transformation is also monotone increasing. Also, the sum $\sum_{b=1}^{B}\lambda_{b}$ determines the range of the monotone spline transformation function, so we constrain it to equal a constant $c$ for identifiability.

Choosing the number and placement of the spline knots is generally a challenging problem. In this paper, we simply place three interior knots at evenly spaced quantiles (James et al., 2013, Section 7.4.4). Adding more knots increases the flexibility of the transformation, but also increases the number of parameters that will need to be optimized in our model-fitting procedure. In our application, we found that three interior knots provided sufficient flexibility while maintaining computational efficiency and stability. Appendix A.1 provides further guidance about knot choice.

The full DLM with transformations is therefore as follows:

where we define ${\bm{\psi}}_{t}\equiv\tau_{\bm{\lambda}}(\tilde{\mathbf{y}}_{t})$, the transformed, game-centered observations, suppressing dependence on ${\bm{\lambda}}$ to simplify notation. To complete the model specification, we specify prior distributions for ${\bm{\theta}}_{1}$, $\sigma^{2}$, $w$, and ${\bm{\lambda}}$:

The hyperparameters $v_{0}$, $a_{0}$, $b_{0}$, $s_{w}$, $c$, and $\bm{\alpha}$ may be set to reflect prior beliefs about ${\bm{\theta}}_{1}$, $\sigma^{2}$, $w$, and ${\bm{\lambda}}$. In the absence of available prior information, we set $v_{0}=10$, $a_{0}=b_{0}=0.1$, and $s_{w}=1$ to keep the priors fairly uninformative about ${\bm{\theta}}_{1}$, $\sigma^{2}$, and $w$ (Gelman et al., 1995). We set $\bm{\alpha}$ to be proportional to the $\lambda$ parameters corresponding to the identity transformation, but keep $\sum_{b=1}^{B}\alpha_{b}=1$ to keep the prior diffuse. If more shrinkage toward the identity transformation is desired, $\sum_{b=1}^{B}\alpha_{b}$ can be set to a higher value. Finally, we set $\lambda_{0}$ to equal the lowest score in the data and $c$ to equal the range of the scores in the data so that the learned transformation roughly preserves the scale of the original data.

The Dirichlet prior on ${\bm{\lambda}}$ constrains the transformation parameter ${\bm{\lambda}}$ to have nonnegative components that sum to $c$. In practice, we find that using a weakly regularizing unconstrained prior for ${\bm{\lambda}}$ (Kowal and Canale, 2020) often results in better performance:

where $N^{+}$ indicates a normal distribution truncated below at 0. Theoretically, the shape of the monotone spline transformation is unidentifiable without a constraint on $\sum_{b=1}^{B}\lambda_{b}$; empirically, the mild regularization induced by the truncated normal priors effectively addresses this issue. For this unconstrained model, we set $\bm{\alpha}$ equal to the $\lambda$ parameters corresponding to the identity transformation and control shrinkage toward the identity transformation using $s^{2}_{\bm{\lambda}}$. To allow maximum flexibility, we set $s^{2}_{\bm{\lambda}}$ to a large value.

Figure 1 displays the graphical model representing the relationships between the model parameters $\{{\bm{\theta}}_{t}\}_{t=1,\dots,T}$, $\sigma^{2}$, $w$, and ${\bm{\lambda}}$ and the transformed data $\{{\bm{\psi}}_{t}\}_{t=1,\dots,T}$.

Figure 1 displays the conditional independences implied by Equations 4 and 5, which allow us to factor the joint distribution of our untransformed data and model parameters as:

where $J({\bm{\psi}}_{1:T}\to\mathbf{y}_{1:T})$ is the Jacobian of the inverse transformation, $\tau_{\bm{\lambda}}^{-1}(\cdot)$. The Jacobian term is vital to appropriately account for how the transformation rescales the data. For example, the Jacobian of the inverse monotone spline transformation is:

The I-spline basis functions used for the monotone spline are constructed by integrating M-spline basis functions; here, $M_{b}(\cdot)$ are the M-spline basis functions corresponding to their respective I-splines. Figure 2 displays the seven spline basis functions constructed for the biathlon relay training dataset we study in Section 4. Figure 2 shows how each I-spline basis function is constructed as the integral of an M-spline basis function.

While the general multi-competitor setup introduced above can be used for games with any number of players, we introduce a slightly simpler setup for the special case of head-to-head games. In head-to-head games, it is generally more natural to consider monotone transformations of score differences, which have no direct analogues in the multi-competitor setting. Instead of modeling the $2g_{t}\times 1$ vector of athlete scores $\mathbf{y}_{t}$, we can model the $g_{t}\times 1$ vector of score differences $\mathbf{z}_{t}$, where we subtract the second athlete’s score from the first athlete’s score. The resulting model likelihood is:

Note that we have simply replaced the model matrix $\bar{X}_{t}$ with the model matrix $Z_{t}$, which we define as a matrix with two nonzero entries per row: 1 in the column corresponding to the first athlete and -1 in the column corresponding to the second athlete. Each observed score difference is taken as the first athlete’s score minus the second athlete’s score, so that each transformed score difference is normally distributed around the latent ability parameter difference between the two competing athletes.

We estimate the model parameters $({\bm{\theta}}_{1:T},\sigma^{2},w,{\bm{\lambda}})$ via a two-step procedure. In our application of interest, we would like to be able to quickly update athlete ratings (i.e., latent ability estimates) shortly after the results from a game. While we could estimate the full posterior distribution of all of the model parameters after each game, this may be an unnecessarily slow and computationally expensive process. Instead, we fit the model using the following procedure:

1. 1.

   Estimate $w$ and ${\bm{\lambda}}$: using a training subset consisting of the first $T_{train}$ rating periods in the dataset, obtain estimates $\hat{w}$ and $\hat{{\bm{\lambda}}}$.

2. 2.

   Estimate ${\bm{\theta}}_{1:T}$ and $\sigma^{2}$: given $\hat{w}$ and $\hat{{\bm{\lambda}}}$, use the full dataset to obtain estimates $\hat{{\bm{\theta}}}_{1:T}$ and $\hat{\sigma}^{2}$.

While step 1 is computationally expensive, step 2 can be implemented efficiently. When the results from a new game become available, we can just run step 2 using the previously learned values of $\hat{w}$ and $\hat{{\bm{\lambda}}}$ to quickly update ratings.

We first describe step 2, the final estimation of ${\bm{\theta}}_{1:T}$ and $\sigma^{2}$. Given fixed values for $w$ and ${\bm{\lambda}}$, we can transform the full dataset using $t_{{\bm{\lambda}}}(\cdot)$ and then use standard Kalman filter equations for a model with an unknown constant variance parameter to estimate the posterior distributions of ${\bm{\theta}}_{1:T}$ and $\sigma^{2}$ (Prado and West, 2010, Section 4.3.2). The posterior distributions of ${\bm{\theta}}_{t}$ and $\sigma^{2}$ after each time step $t$ can be expressed as:

where $\mathbf{m}_{t}$, $V_{t}$ $a_{t}$, and $b_{t}$ are computed using the transformed data as:

where we initialize $V_{0}=v_{0}\cdot I_{p}$ and $\mathbf{m}_{0}=\mathbf{0}$.

It is worth noting that the $V_{t}$ matrix allows the variance of ${\bm{\theta}}_{t}$ to vary from athlete to athlete as a function of sample size. Athletes with many observed scores in a particular rating period will have less posterior uncertainty around their latent ability parameter for that period. The observation variance parameter $\sigma^{2}$, on the other hand, is driven by residual variation across all athletes, as demonstrated by the update equation for $b_{t}$. Finally, we note that under this framework, if an athlete does not compete in a given time period, the estimated mean of their latent ability parameter does not change. The variance of their latent ability parameter, however, still increases by $w\cdot\sigma^{2}$ to reflect the increased uncertainty about their current ability.

To estimate $w$ and ${\bm{\lambda}}$ for step 1, we use the marginal posterior density of $w$ and ${\bm{\lambda}}$ where ${\bm{\theta}}_{1:T}$ and $\sigma^{2}$ are integrated out, which can be expressed in closed form (up to a normalizing constant) as:

for the priors on $w$ and ${\bm{\lambda}}$ and posterior predictive densities:

where the $V_{t}$, $\mathbf{m}_{t}$, $a_{t}$, and $b_{t}$ values are all computed using the Kalman filter as described above. We can obtain samples of $w$ and ${\bm{\lambda}}$ from $p(w,{\bm{\lambda}}\mid\mathbf{y}_{1:T_{\text{train}}})$ using standard Markov Chain Monte Carlo (MCMC) methods. We implement our model in Stan (Stan Development Team, 2021), which uses Hamiltonian Monte Carlo and a no-U-turn sampler (Hoffman et al., 2014).

In practice, we can instead take a maximum a posteriori (MAP) approach using standard optimization routines to greatly reduce the computational burden. Note that each evaluation of $p({\bm{\psi}}_{t}\mid{\bm{\psi}}_{1:t-1},w,{\bm{\lambda}})$ for $t=1,\dots,T_{\text{train}}$ in the marginal posterior density $p(w,{\bm{\lambda}}\mid\mathbf{y}_{1:T_{\text{train}}})$ requires a Kalman filter to be run on the full training dataset. MCMC sampling requires many such likelihood evaluations, so even though each Kalman filter run is relatively efficient, MCMC sampling ultimately requires significant computation. For most applications, however, we are only interested in estimating reasonable values for $w$ and ${\bm{\lambda}}$, not in conducting full Bayesian inference on them; as such, it often makes sense to simply find the posterior mode of $p(w,{\bm{\lambda}}\mid\mathbf{y}_{1:T_{\text{train}}})$, which requires fewer evaluations of $p({\bm{\psi}}_{t}\mid{\bm{\psi}}_{1:t-1},w,{\bm{\lambda}})$ for $t=1,\dots,T_{\text{train}}$. Importantly, integrating ${\bm{\theta}}_{1:T_{\textbf{train}}}$ and $\sigma^{2}$ out of the posterior instead of maximizing $p(w,{\bm{\lambda}},{\bm{\theta}}_{1:T_{\text{train}}},\sigma^{2}\mid\mathbf{y}_{1:T_{\text{train}}})$ with respect to each parameter produces a better-informed posterior mode of $w$ and ${\bm{\lambda}}$.

Any nonlinear optimization algorithm can be used to obtain MAP estimates. For the constrained optimization (with a Dirichlet prior on ${\bm{\lambda}}$), we choose to use the Augmented Lagrangian Adaptive Barrier Minimization Algorithm (Varadhan, 2015). For the unconstrained optimization (with normal priors on the $\lambda_{b}$ parameters), we find that the Nelder-Mead algorithm (Nelder and Mead, 1965) generally gives the most stable results, though the L-BFGS algorithm (Liu and Nocedal, 1989) is much faster in practice. While these optimization-based approaches may theoretically get stuck in local modes, we find that in practice they produce reasonable and effective results.

We implement the full model-fitting procedure in the dlmt package in R. Full implementation details for our model-fitting procedure can be found in Appendix A.1, and a simulation study confirming appropriate parameter recovery can be found in Appendix B.

The Kalman filter equations produce estimates of the filtered latent ability parameter distributions $p({\bm{\theta}}_{t}\mid{\bm{\psi}}_{1:t},\sigma^{2},{\bm{\lambda}})$. In the second stage of the model-fitting procedure, we may also wish to calculate the smoothed latent ability parameter distributions $p({\bm{\theta}}_{t}\mid{\bm{\psi}}_{1:T},\sigma^{2},{\bm{\lambda}})$, where the full dataset informs each latent ability parameter estimate. The Rauch-Tung-Striebel smoother (Rauch, Tung and Striebel, 1965) provides a simple algorithm for doing so. We can compute the smoother updates as:

for:

for scaling matrix $S_{t}=V_{t}(V_{t}+wI)^{-1}$, where the $\mathbf{m}_{t}$ and $V_{t}$ values are the original $\mathbf{m}_{t}$ and $V_{t}$ values computed using the Kalman filter equations. The smoother updates are computed starting from rating period $t=T$ backwards to rating period $t=1$, starting from $\mathbf{m}_{T}^{s}=\mathbf{m}_{T}$ and $V_{T}^{s}=V_{T}$.

We illustrate our model on a variety of Olympic sport datasets provided to us by the US Olympic and Paralympic Committee. The data roughly span from 2004 to 2019 and include the score outcomes from selected national and international competitions. We briefly describe each dataset below:

We fit our unconstrained model to the biathlon, biathlon relay, diving, fencing, and rugby datasets. We conduct full MCMC (four chains, 1000 burn-in iterations, 1000 samples) as well as MAP estimation using the Nelder-Mead and L-BFGS optimization algorithms. To assess the convergence of our MCMC estimates of the posterior distributions of $w$ and ${\bm{\lambda}}$, we check trace plots and the $\hat{R}$ diagnostic for Hamiltonian Monte Carlo. Visual inspection of the trace plots does not suggest any evidence of non-convergence, and $\hat{R}$ is nearly equal to one for all model parameters. The Nelder-Mead and L-BFGS algorithms converge under default convergence tolerances.

Before examining the DLM results, we first confirm that the MAP estimates produce similar results as full MCMC. Figure 3 compares the posterior mean of the transformations learned using full MCMC to the transformations learned using MAP with the Nelder-Mead and L-BFGS algorithms.

We see that the learned transformations are essentially the same across all three algorithms for all five sports. The learned $w$ parameters (not shown) are also very similar. As a result, using MAP methods rather than full MCMC does not significantly impact model performance, while it can lead to significant speedups. For example, for the biathlon dataset (~6000 observations from ~700 athletes over ~60 events), full MCMC took roughly 8 hours, but the Nelder-Mead optimization took only 30 minutes, and L-BFGS converged in less than one minute (all on a laptop with an Intel Core i7-8550U CPU). To simplify assessment, visualization, and discussion of our results, we will focus on the transformations estimated using full MCMC moving forward, though results are naturally similar for the transformations estimated using MAP.

Next, we assess the fit of the DLM on our transformed datasets. To do so, we use the first two-third rating periods as a training set to learn an appropriate transformation and leave the last one-third rating periods as a test set. We then visualize $({\bm{\psi}}_{t}-\bar{X}_{t}\mathbf{m}_{t-1})$, i.e., the one-step prediction residuals, on the test set. If the learned transformation is effective, the residuals should be approximately normally distributed. Figure 4 shows Q-Q plots of standardized test-set residuals against standard normal quantiles.

While the residuals show some slight outliers at the extremes, they are largely normally distributed. The learned transformations thus appear to produce reasonable fits of the standard DLM to these datasets.

We use the biathlon and rugby datasets to illustrate the results from our model. For both datasets, we learn an unconstrained monotone spline transformation using MCMC, transform the dataset, and run the Kalman filter on the transformed data.

Figure 5 shows the smoothed biannual latent ability point estimates for the 25 biathletes with the most biathlons entered. In the biathlon, low race times are better, so negative ability parameters indicate strong athletes.

Of the visualized latent ability trajectories, two stand out, and Figure 5 additionally shows their 90% central posterior probability intervals. The trajectory in blue belongs to the “King of Biathlon,” Ole Einar Bjørndalen, the winningest biathlete of all time at the Olympics, Biathlon World Championships, and the Biathlon World Cup tour. Though the scope of our dataset does not include the start of his career in the 1990s, the model clearly notes his dominance in the early 2000s. The trajectory in red belongs to Martin Fourcade, who began serious international competition in 2008 (hence the wide probability intervals prior to 2008) and proceeded to put together a record-breaking string of seven overall World Cup titles in a row from 2011 to 2018. Bjørndalen and Fourcade are the two names considered in discussions of the greatest male biathlete of all time. Our model notes Fourcade’s quick rise to prominence, projecting that he would begin to outperform Bjørndalen as early as 2009, though interestingly it never projects Fourcade’s latent ability to exceed Bjørndalen’s peak latent ability in 2005.

Figure 6 shows the smoothed quarterly latent ability point estimates and corresponding 90% posterior intervals for the national rugby teams of Fiji and New Zealand, two countries popularly known for their strong rugby teams. Here, more positive latent ability estimates represent stronger teams.

While the two teams have similar estimated strengths during the time span of the data, we see that the model briefly estimates New Zealand to be significantly stronger in 2007 and 2013, and Fiji to be significantly stronger in 2015.

Another result of interest is the monotone spline transformation learned by the DLM with transformations. Figure 7 displays 100 posterior transformation samples for each of the five sports, under the unconstrained optimization.

The transformations generally reflect conventional wisdom about scores from these sports. For example, we sometimes see very slow race times in the biathlon and biathlon relay, which occur when an athlete makes a few shooting mistakes and takes penalties. We expect to see a negative feedback loop in the biathlon where taking penalties causes athletes to get more frustrated or tired from penalty laps, which causes more penalties. This means that extremely slow race times may not reflect extremely poor skill. The learned transformation shrinks the very slow race times to be less extreme, which intuitively helps to make them better reflect athlete skill in the normal DLM. On the other hand, the transformation learned for fencing magnifies extreme score differences. In a fencing bout, points are scored one-at-a-time; after each touch (i.e., point scored), the fencers reset to their starting positions. This makes very one-sided matches relatively rare, since even outmatched fencers can typically score some number of lucky points in a match to fifteen points. The transformation notes this and indicates that when a fencer wins by many points, they are much stronger than their opponent, even more so than the large score gap may suggest.

Table 2 shows the posterior means of the $\sqrt{w}$ and $\sigma$ parameters for each of the five sports. Recall that $\sigma^{2}$ represents the observation variance and $w$ represents the ratio of the innovation variance to the observation variance, so $\sigma$ and $\sigma\sqrt{w}$ would represent the observation and innovation standard deviations, respectively. For example, in the biathlon data, the model estimates that the innovation standard deviation is approximately 28% as large as the observation standard deviation of 98.5. We make two notes about the $\sqrt{w}$ values estimated in these datasets. First, we see that $\sqrt{w}<1$ in all five datasets, i.e., that the athletes’ period-to-period variability in their latent abilities is generally much smaller than the game-to-game variability in their scores. The competitors in our datasets are all professional-level athletes. As such, while they may still improve or worsen over time, their game scores are largely driven by how well they perform in each particular game rather than by major changes to their skill level. Secondly, we note that the $\sqrt{w}$ values vary across sports. In addition to capturing period-to-period variability, the value of $\sqrt{w}$ also reflects the extent to which skill, rather than random chance, appears to drive game results; for example, a game with results driven almost entirely by skill would have minimal game-to-game variation, i.e., $\sigma\approx 0$, so any period-to-period variation in athlete abilities would imply a large value of $\sqrt{w}$. The results in Table 2 therefore suggest that of the five sporting events, the results from diving competitions may be the most driven by skill. On the other hand, the random chance associated with particular performances appears to play a significant role in the outcomes of fencing bouts, which has been noted by other sources as well (e.g., Zappalà et al. (2022)).

To evaluate the DLM with transformations, we compared the accuracy of its predictions to predictions made by other models for multi-competitor and head-to-head athlete rating. For multi-competitor sports, we compared the DLM with transformations (LM-T) to the DLM without transformations (LM) and the dynamic rank-order logit model (ROL) from Glickman and Hennessy (2015). For head-to-head sports, we compare to the Glicko rating system (GLO; Glickman, 1999). We use the first two-third rating periods in each dataset as a training set to tune the model hyperparameters $w$ and ${\bm{\lambda}}$, fit the model on the full dataset, and finally evaluate its predictions for the test set (i.e., the last one-third rating periods). For multi-competitor games, we evaluate predictions using the Spearman correlations between the observed and predicted athlete rankings in each game. This approach was taken in Glickman and Hennessy (2015) to evaluate predictability on a test set. We summarize these game correlations $\rho_{tg}$ over the test set using a game-size-weighted average (Glickman and Hennessy, 2015):

Note that we limit ourselves to using rank-based metrics to facilitate comparison with the ROL model, which predicts athlete ranking probabilities rather than scores. For head-to-head games, we evaluate predictions using the average accuracy of winner predictions in the test set.

Tables 3 and 4 show the weighted Spearman correlation (Equation 6) of the ranking predictions for the multi-competitor sports and the accuracy of winner/loser predictions for the head-to-head sports.

Across the five datasets, we see some evidence that the LM-T model outperforms the other models in terms of predictive performance. While it is difficult to generally evaluate the relative empirical performance of the LM-T model with only five datasets, we see that it improves predictive accuracy across nearly all of the datasets. The only exception is the rugby dataset, where the LM-T model essentially learns an identity transformation, so we would not expect it to outperform the LM model. The results of athletic competitions are generally challenging to predict; small improvements in predictive performance, such as those demonstrated above, can therefore be fairly valuable for generating a competitive edge.