Plackett–Luce modeling with trajectory models for measuring athlete strength

2.1 Plackett–Luce modeling

For competition k within time period t, assume athletes 1, 2, …, m _kt compete. Assuming no ties, let R _i , i = 1, …, m _kt be the rank of athlete i from the results of the competition. That is, R _i = j is equivalent to athlete i coming in jth place. Suppose { i 1 , i 2 , … , i m k t } is a permutation of {1, 2, …, m _kt }. Then the Plackett–Luce model (Luce 1959; Plackett 1975) asserts that

where θ t = ( θ 1 t , … , θ m k t t ) , and θ _it for i = 1, …, m _kt is the strength parameter for athlete i at time period t. The model in (1) is essentially a telescoping product of multinomial logit probabilities in which the first term is the probability player i₁ outperforms everyone else in the competition; the second term is the probability that player i₂ outperforms everyone else in the competition excluding player i₁; the third term is the probability that player i₃ outperforms everyone else in the competition excluding players i₁ and i₂, and so on.

The Plackett–Luce model can be derived through a latent variable framework. This representation will be helpful in Section 4.1 when evaluating our fitted model. Suppose the performance of athlete i at time t with strength parameter θ _it follows a specified probability distribution. In particular, let Y _itk be the unobserved latent performance by competitor i at time t during competition k. If we assume a cumulative distribution function

that is, the Y _itk follows a Gumbel distribution, then the probability athlete i₁ outperforms i₂ who outperforms i₃, and so on, can be shown to be

This result is due to Plackett (1975), and can be easily shown using basic probability calculus. The likelihood contribution for multiple competitions is simply the product of (1) over the competitions.

2.2 Trajectory model

A growth curve, to which we refer in this manuscript as a trajectory model, is any parametric function that describes how a quantity changes over time. Shapes can range widely and may take the form of linear, exponential, logistic, or S-shaped models. Trajectory models can also include more complicated relationships like mixed-effects or nonlinear regression (Panik 2014; Pinheiro and Bates 2000). Growth-curve models are well-proven and have been utilized for prediction in many fields, including sports-analytics literature. However, the trajectory models commonly used to predict athlete strengths are too simple to account for complicated changes in performance over time, such as flexible career trajectories impacted by injuries. Models typically use a quadratic function to find the time of the athletes’ optimal performance, such as Brander, Egan, and Yeung (2014) with National Hockey League data, Malcata, Hopkins, and Pearson (2014) with triathlete data. Examples of more general trajectory models include Moudud et al. (2008) who use the nonlinear logistic trajectory model with a mixed-effects trajectory model to predict the speed of youth cross country skiers, and Bradlow and Fader (2001) who use the generalized gamma curve to model the ranking of songs on Billboard charts. Both of these trajectory models use an exponential decay term with the expectation that, over time, performance converges to a specific value or decays altogether, a characteristic we expect to be shared with an athlete’s career trajectory.

The main alternative to parametric trajectory models is to assume a stochastic process on the θ _it . This is the approach taken by GH. Specifically, the GH model assumes a Gaussian random walk on the θ _t of the form

where τ² is the innovation variance of the Gaussian random walk. Other assumptions on the stochastic process in (4) may be made, such as the inclusion of an autoregression factor, but as shown in Glickman and Hennessy (2015) the basic random walk works well in practice. A challenge using the GH approach is in the difficulty of forecasting beyond several time periods. The model assumes that the expected ability of an athlete at time t + 1 is centered at their ability at time t.

Instead of including a stochastic process component, we can extend the Plackett–Luce model by including a trajectory model component. Trajectory models are useful for modeling career trajectories because they can be constructed in a flexible manner. While a quadratic function describes the typical athlete trajectory, increasing performance until reaching their peak and then declining in performance, it does not allow for more complex trajectories such as careers with several peaks due to injuries or adoptions of new strategies. To encourage flexibility in the model, the trajectory model in Equation (5) allows for a polynomial of varying degrees. We assume a novel parametric growth-curve model on θ _it as a function of the current period t and the period in which the athlete first competed t_0i. The model we assume on θ _it is

where the bracketed term is a pth order polynomial. Intercepts and coefficients α _i and β _i = {β_i0, β_i1, β_i2, …, β _ip } vary per individual, and the rate, ω > 0 is fixed across individuals. The β _i and α _i are competitor-specific random effects. The intercepts α _i not only adjust for different starting strength, but they also represent the limit of decay over a long career. The coefficients β_i0, β_i1, …, β _ip and intercept α _i parameters differ by athlete, since we do not expect each athlete to have the same career trajectory.

The stabilization factor of ω shows the performance decaying over time. This concept is intuitive because of the physical nature of sports, an athlete’s performance decays with age. The parameter ω remains fixed across athletes because we believe the overall stabilization factor of an athlete’s performance should be consistent between athletes in the same sport, on average. The polynomial component provides enough flexibility to account for variation in stabilization factors on an individual level.

The order of the polynomial is a free parameter p that needs to be chosen via model selection as described in Section 3. To reduce correlation among polynomial orders of time, we use orthogonal polynomials of t − t_0i = 0, …, T − t_0i, denoted as t 1 * , t 2 * , … , t p * . Orthogonal polynomials are calculated across the entire data set by creating an orthogonal basis from a QR decomposition via Gram–Schmidt orthogonalization. The orthogonal polynomials are calculated through time T, and for forecasts we can continue the Gram–Schmidt orthogonalization on times T + 1, T + 2, …. See Appendix A for details on the creation of the orthogonal polynomials in fitting and prediction. Raw time is used in the exponent since orthogonal polynomials are not guaranteed to be non-negative. To ensure identifiability of the strength parameters, we center the θ _it at every time t by requiring ∑ i = 1 n θ i t = 0 .

The intercept and coefficient parameters are modeled as random effects, allowing for shrinkage toward a population mean for the parameters of athletes with few observations. Malcata, Hopkins, and Pearson (2014) chose to model performances of triathletes as a function of age with coefficients for random effects so that each athlete had a different set of parameters. Bell et al. (2016) used multilevel modeling on Formula 1 race results to account for team and driver effects, and identified the best racers over time. For our model, the polynomial parameters are assumed to come from a common normal distribution. The α _i are normally distributed with a mean of 0 and a variance of σ α 2 . The β _ib coefficients are drawn from independent normal distributions with different means η β b and variances σ β b 2 corresponding to the different polynomial orders b = 0, 1, …, p.

We choose to keep priors non-informative or weakly informative. For the intercept and coefficient parameters, the support of the prior must be a continuous distribution that spans all real numbers, e.g., a normal distribution. The prior distribution for ω should span all non-negative or positive numbers, and could take the form of an exponential or gamma distribution.

2.3 Adding covariates

In some cases, (e.g., the case study in Section 4) we may be able to collect more player-specific information. We illustrate two examples of how such covariates may be incorporated into the growth-curve model described in Equation (5).

If available, variables such as weight, height, and equipment specifications that change over time could be useful in predicting performance. For example, youth athletes who grow taller more quickly may increase in performance more quickly than expected compared to their slower-growing counterparts. In this example, time-varying, player-specific covariates X _it can be included in the trajectory model as an additive term,

where γ is the vector of coefficients on the covariates describing linear changes over time.

Another variable we might include is player age. The older an athlete, the less time they have in their career and the more quickly their performance stabilizes. However, the impact varies from sport to sport. For example, in target shooting, skill is usually a far superior predictor of performance than age, which is why these sports boast the oldest Olympic athletes with the longest careers. Conversely, gymnastics is a sport that requires an enormous amount of flexibility and dexterity that tends to be only attainable by young athletes. We assume that the coefficient for age is fixed within a sport. Our model incorporates the age of an athlete’s first professional appearance, z _i , as

where ω₁ ≥ 0, ω₂ ≥ 0 are the linear coefficients describing the stabilization factor. This form of trajectory model is used to fit women’s luge data in Section 4. Another option may be to make our trajectory model a function of age rather than a function of time. Although directly incorporating age into the trajectory model would be straightforward, we choose to index the trajectory as a function of time since the start of the athlete’s professional career and include age as a covariate. Indexing by length of career instead of age makes it easier to compare career trajectories across athletes of different ages with more overlap.

4.1 Evaluation

To evaluate the predictive ability of our model, we use a weighted average of the Spearman rank correlation (Spearman 1904) between our projected latent performance and the true race results. This method is also used to evaluate the GH model in Glickman and Hennessy (2015). The reason for using Spearman rank correlation instead of the posterior predictive is to have an evaluation metric that is robust against the large variation in the posterior predictive due to a large number of competitors and a small number of competitions.

The weighted average of the Spearman rank correlation at predicted time T is:

where K _T is the number of events at time T and m _kT is the number of athletes competing in event k = 1, …, K _T . We evaluate ρ _W at each iteration of the posterior draw. To estimate ρ ̂ k we first evaluate θ _iT at desired time T, for each player that appears in event k using the posterior draws of the trajectory model parameters. We then draw the latent performance Y _iT from the Gumbel distribution with location θ _iT . We can then calculate the Spearman rank correlation between the true event results and the sampled latent variables to obtain ρ ̂ k .

Table 2 shows the evaluation for the women’s luge example in several settings. First, we remove only the final year of results for each athlete and fit the growth-curve model. We use this model built with missing data to predict the results for matches in the two out of sample observed periods of that year.^[1] We then remove the final 2 years of results to fit the data and predict the results for matches in the four observed periods across both years. Figure 2 shows the estimated weighted Spearman correlation and 95 % credible interval across the four periods within the 2 years of withheld information. Notice that the correlation does not drop, even evaluated out greater than four periods, showing that our methods are not just valuable for current predictions, but can extrapolate far beyond the end of the observed data.

Figure 2:

The estimate of the weighted Spearman correlation and 95 % CI. The square points show the performance of the GH model.

We compare the Spearman correlation of our trajectory model to that of the GH model calculated on the same data. Table 2 shows the Spearman correlation on the test data with 2 years (four periods) held out. The Spearman correlation from using the GH model is shown by the gold boxes in Figure 2). Our model performs similarly across all predictions, even up to 8 periods where the GH model has more variation from period to period as well as has a drop in performance in the last two periods compared to the first two. This highlights the advantage of using the trajectory model with Plackett–Luce, which allows us to include not only the estimated performance at a given time, but how we predict that might change as time progresses.

4.2 Prediction summaries

Figure 3 shows the fitted player strength parameters with 95 % credible intervals for four of the top athletes in luge: Alex Gough (CAN), Natalie Geisenberger (GER), Summer Britcher (USA) and Tatjana Hüfner (GER). These four athletes started competing at different dates so they were at different stages of their career. Hüfner, the oldest, was beginning to decline in performance, whereas Britcher was the youngest athlete and had a sharp increase in performance over the past 3 years. Geisenberger has had a steady career and always performed at the top, whereas Gough took several years to reach a similar strength.

Figure 3:

Growth-curve model fit with 95 % credible intervals to the observed time range for four of the top women’s luge athletes. Plotting the trajectory models against one another lets us compare the different levels of performance between the athletes at a given time.

Figure 4 shows the same trajectories but aligned by years since the start of each athlete’s career t − t_0i. This view allows us to compare trajectories between athletes while they are at the same points in their career. As we can see, the top luge athletes share similarly shaped trajectories. One point of interest for coaches or team managers would be comparing younger players’ career trajectories to the trajectories of their more established counterparts. In Figure 5, we show Britcher’s career projected out 7 years compared to Geisenberger’s current career trajectory. By plotting the trajectories with their respective 95 % credible intervals side-by-side, aligned at the beginning of their careers, we can see that they have similarly shaped trajectories, but that Britcher never reaches the level of performance of Geisenberger.

Figure 4:

The same trajectory model fits with 95 % credible intervals as in Figure 3 but aligned at the beginning of the career.

Figure 5:

Career trajectory and 95 % credible interval of Natalie Geisenberger next to the projected career trajectory and 95 % credible interval of Summer Britcher aligned at the start of their careers. The gray vertical line shows the current point, whereas to the right is the predicted strength for Britcher.

We can also use our model fit to understand the general trajectory of athletes in a particular sport. To do so, we take advantage of the hierarchical structure of the polynomial coefficients and construct a curve using the estimated means,

In the women’s luge case we incorporate the average starting age across all players z ̄ ,

Figure 6:

The mean curves for women’s luge (left) and men’s slalom (right). By looking at the mean curves, we can compare the typical trajectory of athletes between two sports. Note that the y-axis differs between the two images.

4.3 Clustering career trajectories

In Section 4, we compare sports by estimating the mean trajectory model. Here, we perform an exploratory analysis to understand how athletes share similarities to one another. An exploratory way of finding similarities is to cluster career trajectories. Bartolucci and Murphy (2015) fits a finite mixture model on runners in a 24-hour-long endurance race with the intent of identifying strategies amongst clusters of runners based on their performance. Such an analysis can be useful in finding athletes with a potentially similar style or technique, or help understand what makes athletes more successful during a certain career stage. To enable the exploration of the answers to these questions, we cluster the estimated career trajectories of the athletes. Instead of clustering by the predicted rating, we cluster the estimated trajectory model parameters. Clustering based on estimated trajectory model parameters has the advantage of being able to describe the athlete’s entire career trajectory, rather than just the observed years.

Many techniques for clustering trajectories can be used, including two-stage, nonparametric, and model-based methods (Jacques and Preda 2014). In particular, we focus on two-stage methods which first reduce the dimensionality of the data by providing a functional basis, and then use a non-parametric clustering approach to cluster the observations using this simplified basis. Abraham et al. (2003) first fits B-splines to the longitudinal data set and uses k-means clustering to the corresponding parameters. In an application to sports, Miller and Bornn (2017) cluster segments of the National Basketball Association athletes’ movements during possession by mapping them to a set of characteristic actions.

Instead of clustering against point estimates of the parameters, we want to utilize the entire distribution of these parameters to retain the uncertainty of the estimates. The following definition of the Wasserstein distance is a good metric for measuring the distance between two empirical distributions, which in our case are samples from the posterior. For simplicity, we use the distance between the marginal coefficients, although we recognize we ignore correlation. If we define B _ib as the empirical distribution of the posterior of β _ib , for athlete i and B _jb to be the equally sized empirical distribution of β _jb for athlete j then the distance between the two distributions is defined as a function of the ordered posterior draws,

where β^(k) indicates the kth ordered draw. The total distance between two athletes is defined as

where P is the order selected through model selection. The β _ib are centered by η _p and standardized by σ β b . The intercept terms β₀, α are ignored since we are more interested in the shape of the trajectories in this clustering exercise rather than the value of the rating.

We perform agglomerative hierarchical clustering using the distance matrix created via the Wasserstein distance (Ward 1963). Hierarchical clustering methods are advantageous because they enable us to perform clustering before determining the number of partitions. We use the average silhouette method to determine the number of clusters (Rousseeuw 1987). The optimal number of clusters for the women’s luge trajectories is two clusters, but we will investigate the second most optimal number of clusters, four clusters. Table 3 records the number of athletes within each cluster for both clustering schemes. Details about the clustering methods and finding the optimal number of clusters can be found in Appendix B.

Figure 7 shows the mean and interquartile range (IQR) of the estimates of the career trajectories within each cluster at each period, taken from the fitted trajectory models from all men’s slalom athletes. Immediately, we can see the separation between groups of athletes by simply clustering on three parameters. The clustering method appears to divide the trajectories into athletes with mostly increasing ratings (cluster 1, 3), athletes with ratings that increase and then decrease (cluster 2), and athletes that remain the same or increase slightly (cluster 4). Clusters 1 and 3 contain athletes whose ratings increase steadily throughout their career, yet cluster 3 appears to have athletes who increase at a quicker rate. The athletes in cluster 2 at first increase as quickly as those in cluster 3, but then at some point begin to level off or decrease in rating. Cluster 4 appears to have athletes who have a lower rating than athletes in the other clusters and don’t change much over time.

Figure 7:

Mean and IQR of the estimated career trajectory for men’s slalom athletes within each of the four clusters.

Figure 8 shows the section of the dendrogram and corresponding projected career trajectories for the athletes in cluster 2. The trajectories feature a similar shape, increasing ratings in the earlier part of their career and then leveling off or decreasing in rating after about 5 years. We can use the dendrogram to find athletes who are more alike by following the branches from the top (larger distances) to the bottom (smaller distances). For example, Manzenreiter and Oberstolz-Antonova are the last branch in the dendrogram and have roughly the same trajectory shape. Hüfner was the final athlete to be included in the cluster, and has a more dramatic upward trajectory in the first 5 years and downward trajectory in the second 5 years than other athletes in cluster 2. Hüfner and Britcher, whose trajectories we compared in Section 4.2 are both assigned to this cluster and show similar shape.

Figure 8:

The section of the dendrogram representing the path of the clustering and the distance between each of women’s luge athletes within cluster 2 (top). The estimated trajectories and 95 % CI for the members of cluster 2 (bottom).