Distributed lag models to identify the cumulative effects of training and recovery in athletes using multivariate ordinal wellness data

05/18/2020 ∙ by Erin M Schliep, et al. ∙ 0

Subjective wellness data can provide important information on the well-being of athletes and be used to maximize player performance and detect and prevent against injury. Wellness data, which are often ordinal and multivariate, include metrics relating to the physical, mental, and emotional status of the athlete. Training and recovery can have significant short- and long-term effects on athlete wellness, and these effects can vary across individual. We develop a joint multivariate latent factor model for ordinal response data to investigate the effects of training and recovery on athlete wellness. We use a latent factor distributed lag model to capture the cumulative effects of training and recovery through time. Current efforts using subjective wellness data have averaged over these metrics to create a univariate summary of wellness, however this approach can mask important information in the data. Our multivariate model leverages each ordinal variable and can be used to identify the relative importance of each in monitoring athlete wellness. The model is applied to athlete daily wellness, training, and recovery data collected across two Major League Soccer seasons.



There are no comments yet.


page 35

page 39

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The rapid increase in data collection technology in sports over the previous decade has led to an evolution of player monitoring with regard to player performance, assessment, and injury detection and prevention (Bourdon et al., 2017; Akenhead and Nassis, 2016; De Silva et al., 2018). Within both individual and team sports, objective and subjective player monitoring data are commonly being used to assess the acute and chronic effects of training and recovery in order to maximize the current and future performance of athletes (Mujika, 2017; Saw et al., 2016; Thorpe et al., 2015, 2017; Buchheit et al., 2013; Meeusen et al., 2013).

Customary metrics of an athlete’s training response include objective measures of performance, physiology, or biochemistry, such as heart rate, heart rate variability, blood pressure, or oxygen consumption (Borresen and Lambert, 2009). Recently, there has been strong emphasis on also including subjective, or perceptual, measures of well-being in athlete assessment (Akenhead and Nassis, 2016; Tavares et al., 2018). Subjective measures, which are often self-reported by the athletes using self-assessment surveys, include wellness profiles to quantify physical, mental, and emotion state, sleep quality, energy levels and fatigue, and measures of perceived effort during training. Subjective measures have been reported to be more sensitive and consistent than objective measures in capturing acute and chronic training loads (Saw et al., 2016; Tavares et al., 2018). In particular, Saw et al. (2016) found subjective well-being to negatively respond to acute increases in training load as well as to chronic training, whereas acute decreases in training load led to an increases in subjective well-being. Subjective rate of perceived effort has also been reported to provide a reasonable assessment of training load compared to objective heart-rate based methods (Borresen and Lambert, 2008), however accuracy varied as a function of the amount of high- and low-intensity exercise. Tavares et al. (2018) investigated the effects of training and recovery on rugby athletes based on an objective measure of fatigue using countermovement jump tests as well as perceptual muscle-group specific measures of soreness for six days prior to a rugby match and two days following. They found that perceptual muscle soreness tended to diminish after a recovery day while lower body muscle soreness remained above the baseline. In addition, even though running load, as measured by total distance and high metabolic load distance, did not differ between the match and training sessions, the change in perceptual measures of muscle soreness indicated greater physical demands on the athletes during rugby matches. The highest soreness scores for all muscle groups were reported the morning after the match, and these scores remained high the subsequent morning. Interestingly, no significant difference in countermovement jump scores was detected across days. Since no one metric is best at quantifying the effects of training on fitness and fatigue and predicting performance (Bourdon et al., 2017), subjective and objective measures are often used in conjunction to guide training programs to improve athlete performance.

Self-reported wellness measures are comprised of responses to a set of survey questions regarding the athlete’s mental, emotional, and physical well-being. Most commonly, these survey responses are recorded on ordinal (or Likert) scales. Collectively, these ordinal response variables inform on the athlete’s overall well-being. Previous studies looked at the average of a set of ordinal wellness response variables and treated the average as a continuous response variable in a multiple linear regression model

(Gallo et al., 2017). Whereas this approach can be used as an exploratory tool to identify important training variables (work load, duration of matches/games, recovery) on individual wellness, it has three major shortcomings. First and foremost, it throws away important information by reducing the multiple wellness metrics into one value. Second, it assumes that each ordinal response variable is equally important (i.e., assigns weight for each variable). Since there is likely variation in the sensitivity of some wellness variables in an athlete’s response to training and recovery, failing to differentiate between these variables could mask indicators of potential poor performance or negative health outcomes. Lastly, by modeling the average of the wellness metrics, the model is unable to identify important variable-specific relationships between training and wellness.

With the increase in collection of self-reported wellness measures and their identified significance in monitoring player performance, we need more advanced statistical methods and models that leverage the information across all wellness metrics in order to obtain a better understanding of the subjective measures of training and wellness. These models could then be used to synthesize these data in order to guide training programs and player evaluations.

There are many statistical challenges in modeling subjective wellness data. First, the data are multivariate, with each variable representing a particular aspect of wellness (e.g., energy levels, mental state, etc.). Second, the data are ordinal, requiring more advanced generalized linear models of which are not often included in customary statistical programming packages (although see the R package mvord Hirk et al., 2019). Lastly, the subjective wellness data are individual-specific. The day-to-day variation in wellness scores reported by an individual will vary greatly not only between individuals but also across wellness variable for an individual. In addition, training and recovery can have varying short- and long-term effects on athlete wellness, of which are also known to vary across individuals. As such, statistical modeling of individual wellness needs to be able account for the variation across individual as a response to the cumulative effects of training and recovery.

Generalized latent variable models, including multilevel models and structural equation models, have been proposed as a comprehensive approach for modeling multivariate ordinal data and for capturing complex dependencies both between variables and within variables across time and space (Skrondal and Rabe-Hesketh, 2004). To capture flexible, nonlinear relationships, DeYoreo and Kottas (2018)

developed a Bayesian nonparametric approach for multivariate ordinal regression based on mixture modeling for the joint distribution of latent responses and covariates.

Schliep and Hoeting (2013) developed a multi-level spatially-dependent latent factor model to assess the biotic condition of wetlands across a river basin using five ordinal response metrics. A modified approach by Cagnone and Viroli (2018)

used a latent Markov model to model temporal dynamics in the latent factor for multivariate longitudinal ordinal data.

Cagnone et al. (2009) propose a latent variable model with both a common latent factor and auto-regressive random effect to capture dependencies between the variables dynamically through time. Within the class of generalized linear multivariate mixed models, Chaubert et al. (2008) proposed a dynamic multivariate ordinal probit model by placing an auto-regressive structure on the coefficients and threshold parameters of the probit regression model. Multivariate ordinal data are also common in the educational testing literature where mixed effects models (or item response theory models) are used to compare ordinal responses across individual (Lord, 2012). Liu and Hedeker (2006) developed a multi-level item response theory regression model for longitudinal multivariate ordinal data to study the substance use behaviors over time in inner-city youth.

Drawing on this literature, we propose using the latent factor model approach to capture marginal dependence between the multivariate ordinal wellness variables. We extend the current approaches by modeling the latent factors using distributed lag models to allow for functional effects of training and recovery. Distributed lag models (also known as dynamic regression models) stem from the time series literature (see Hyndman and Athanasopoulos, 2018, Chapter 9) and offer an approach for identifying the dynamics relating two time series (Haugh and Box, 1977). Distributed lag models can be written as regression models in which a series of lagged explanatory variables accounts for the temporal variability in the response process. The coefficients of these lagged variables can be used to infer short- and long-term effects of important explanatory variables on the response. Due to the dependent relationship between the lagged values of the process, constraints are often imposed on the coefficients to induce shrinkage but maintain interpretability. Often the constraints impose a smooth functional relationship between the explanatory variables and response. The smoothed coefficients can then be interpreted as a time series of effect sizes and provide insights to the significance of the explanatory variables at various lags. In epidemiological research, distributed lag models have been used to capture the lag time between exposure and response. For example, Schwartz (2000) studied air pollution exposure on adverse health outcomes in humans and found up to a 5 day lag in exposure effects. Gasparrini et al. (2010) developed the family of distributed lag non-linear models to study non-linear relationships between exposure and response. The added flexibility of their model enables the shape and temporal lag of the relationship to be captured simultaneously.

In the ecological context, Ogle et al. (2015) and Itter et al. (2019) presented a special case of distributed lag models to capture so-called ‘ecological memory’. Like distributed lag models, ecological memory models assign a non-negative measure, or weight, to multiple previous time points to capture the possible short- and long-term effects of environmental variables (e.g., climate variables, such as precipitation or temperature) on various environmental processes (e.g., species occupancy). These models are able to identify not only significant past events on the environmental process of interest but also the length of “memory” these environmental process have with respect to the environmental variables. An important distinction between ecological memory models and distributed lag models is that ecological memory models assume the short- and long-term effects to be consistent. That is, with non-negative weights, the relationship between the explanatory variable and the response is the same at all lags governed by the sign of the coefficient. In modeling the short- and long-term effects of training and recovery on athlete wellness, we note that this limitation may be overly restrictive.

We propose a joint multivariate ordinal response latent factor distributed lag model to tackle the challenges outlined above. The joint model specification enables the borrowing of strength across athletes with regard to capturing the general trends in training response. However, athlete-specific model components allow for individualized effects and relationships between wellness measures, training, and recovery. The latent factors capture the temporal variability in athlete wellness as a function of training and recovery. In addition, the multivariate model allows for the latent factors to be informed by each of the self-reported wellness metrics, which alleviates the pre-processing of computing the average and treating it as a univariate response. We model the latent factors using distributed lag models in order to capture the cumulative effects of training load and recovery on athlete wellness.

The remainder of the paper is outlined as follows. In Section 2 we describe the ordinal athlete wellness data, including the metrics of workload and recovery. Summaries of the data as well as exploratory data analysis is included. The multivariate ordinal response latent factor distributed lag model is developed in Section 3. Details with regard to the model specification, model inference, and important identifiability constraints are included. The model is applied to the athlete data in Section 4 and important results and inference are discussed. We conclude with a summary and discussion of future work in Section 5.

2 Athlete wellness, training, and recovery data

Daily wellness, training, and recovery data were obtained for 20 professional referees during the 2015 and 2016 seasons of Major League Soccer (MLS), which spanned from approximately February 1 to October 30 of each year. Each referee followed a training program established by the Professional Referees Organization and attended bi-weekly training camps. Over the span of these two seasons, the number of matches officiated by each referee ranged from 3 to 44, with an average of 28 matches.

Upon waking up each morning, the referee (hereafter, “athlete”) was prompted on his smart phone to complete a wellness survey. The survey questions entailed assigning a value to each wellness variable (metric) on a scale from 1 to 10 where 1 is low/worst and 10 is high/best. These wellness variables include energy, tiredness, motivation, stress, mood, and appetite. Prior to the start of data collection, the referees were trained on the ordinal scoring for these variables such that a high value of all variables corresponds to being energized, fully rested, feeling highly motivated, with low stress, a positive attitude, and being hungry.

Due to the limited number of observations in each of the 10 categories for most individuals and wellness variables, we transformed the raw ordinal data to a 5 category scale. The transformation we elected to use was individual specific, but uniform across metrics as we assumed each individual’s ability to differentiate between ordinal values was consistent across variables. For each individual, we combined the observed ordinal data and conducted -means clustering using 5 clusters. The -means clustering technique minimizes within group variability and identifies a center value for each cluster. The midpoints between the ordered cluster centers were used as cut-points to assign ordinal values from 1 to 5. Each of the six ordinal wellness variables for the individual was transformed to the 5 category ordinal scale using the same set of cut-points. We investigated various transformation (e.g., basic combining of classes 1-2, 3-4, etc, and re-scaling the data to (0,1) and then applying a threshold based on percentiles 0.2, 0.4, etc.) but found model inference in general to be robust to these choices.

Distributions of the 5-category ordinal data for the two metrics, tiredness and stress, are shown in Figure 1 for four athletes denoted Athlete A, B, C, and D. The distributions of wellness scores vary quite drastically both across metrics for a given athlete as well as across athletes for a given metric. The distributions of the other wellness variables for these four athletes are included in the Figure A1 of the Supplementary Material.

In addition to the wellness variables, each athlete also reported information on training and recovery. These data consisted of the duration (hours) and rate of perceived effort (RPE; scored 0-10) of the previous days workout as well as sleep quantity (hours) and quality (scored 1-10) of the previous nights sleep. Distributions of RPE and workout duration are shown in Figure 2 for the four athletes. RPE has been found to be a valid method for quantifying training across a variety of types of exercise (Foster et al., 2001; Haddad et al., 2017). The distribution of RPE varies across each individual both in terms of center and spread as well as with regard to the frequency of “0” RPE days (i.e., rest days). For example, over 40% of the days are reported as rest days for Athlete B. For non-rest days (those with RPE ), Athlete D has a much higher average RPE compared to the other three athletes. The distribution of workout duration tends to be multi-modal for each athlete. In particular, we see spikes around 90 minutes, which is consistent with the length of MLS games. The distribution of shorter workouts for Athlete C tends to be more uniform between 10 and 70 minutes with fewer short workouts than the other athletes.

Similar summaries of the number of hours slept and the quality of sleep are shown in Figure 3

. In general, the average number of hours of sleep for an individual ranges between 6-8 hours, however the frequency of nights with 5 or less hours and 10 or more hours varies significantly across individuals. Sleep quality is most concentrated on values between 6 and 8 for each individual, with the distributions being skewed towards lower values.

We define two important measures of training and recovery to use in our modeling; namely, workload and recovery. Workload, which is also commonly referred to as training load, is quantified as the product of the rate of perceived effort and the duration for the training period (Foster et al., 1996, 2001; Brink et al., 2010)

. Therefore, a training session of moderate to high intensity and average to long duration will result in a large workload value. Recovery is defined using the reported quality and quantity of sleep. We applied principle component analysis on the two sleep metrics for each individual and found that one loading vector captured between 64% and 94% of the variation. As a result, recovery for each athlete was defined using the first principle component, where large values correspond to overall good recovery (high quality and long duration sleep).

3 Joint multivariate ordinal response latent factor distributed lag model

We model the ordinal wellness data using a multivariate ordinal response latent factor distributed lag model. We first describe the multivariate ordinal response model in Section 3.1 and then offer two latent factor model specifications using distributed lag models in Section 3.2. Identifiability constraints are discussed in Section 3.3

, and full prior specifications for Bayesian inference are given in Section

3.4. Section 3.5 introduces important inference measures for addressing questions regarding athlete wellness, workload, and recovery.

3.1 Multivariate ordinal response model

Let denote individual, denote wellness variable (which we refer to as metric), and denote time (day). Then, define to be the ordinal value for individual and wellness metric on day . Without loss of generality, let for each and such that each wellness metric is ordinal taking integer values for each individual.

We model the ordinal response variables using a cumulative probit regression model. We utilize the efficient parameterization of Albert and Chib (1993), and define the latent metric parameter such that


Here, and denote the lower and upper thresholds of ordinal value , for individual and wellness metric , where . Under the general probit regression specification,



. In the Bayesian framework with inference obtained using Markov chain Monte Carlo, this parameterization enables efficient Gibbs updates of the model parameters,

, , and , for all , , and (Albert and Chib, 1993). Posterior samples of the threshold parameters, , require a Metropolis step. More details with regard to the sampling algorithm are given in Section 3.3.

3.2 Latent factor models

In modeling , we propose both a univariate and multivariate latent factor model specification to generate important, distinct inferential measures. We begin with the univariate latent factor model for . Let denote the latent factor at time for individual . The assumption of this model is that is driving the multivariate response for each individual at each time point. That is, for each , , and , we define


where is a metric-specific intercept term and is a metric-specific coefficient of the latent factor individual .

We can extend (3) to an -variate latent factor model where we now assume that the multivariate response might be a function of multiple latent factors. Let denote the latent factors at time for individual . Then, we define as


where captures the metric-specific effect of each latent factor.

We investigate the univariate and multivariate latent factors models in modeling the multivariate ordinal wellness data. The two important covariates of interest identified above that are assumed to be driving athlete wellness include workload and recovery. Therefore, we model the latent factors as functions of these variables. Specifically, we model the latent factors using distributed lag models such that we are able to capture the cumulative effects of workload and recovery on athlete wellness.

Let and denote the workload and recovery variables for individual and time , respectively. Starting with the univariate latent factor model, we model as a linear combination of these lagged covariates. We write the distributed lag model for as


where and are coefficients for the lagged covariates and , and is an error term. Here, we assume . The distributed lag model is able to capture the covariate-specific cumulative effects at lags ranging from . The benefit of the univariate approach is that latent factor offers a univariate summary for individual on day as a function of both wellness and recovery. We can easily compare these univariate latent factors across days in order to identify anomalies in wellness across time for an individual.

The distributed lag model specification can also be utilized in the multivariate latent factor model. Having two important covariates of interest, we specify a bivariate latent factor model with factors and . Here, is modeled using a distributed lag model with covariate , and is modeled using a distributed lag model with covariate . The benefit of this approach is that we can infer about the separate metric-specific relationships with each of the lagged covariates for each individual. That is, we can compare the relationships across metrics within an individual as well as within metric across individuals. For , let


where , and are analogous to above, and is the error term for factor . Again, we assume .

It is worth mentioning that under certain parameter constraints, distributed lag models are equivalent to the ecological memory models proposed by Ogle et al. (2015). That is, ecological memory models are a special case of distributed lag models where the lagged coefficients are assigned non-negative weights that sum to 1. For example, with , the ecological memory model is such that and for all and . Under this approach, is modeled using a Dirichlet distribution. The drawback of this approach is that this forces the relationship between the latent wellness metric and each element of the vector to be the same (e.g., all positive or all negative according to the sign of ). In our application, we desire the flexibility of having both positive and negative short- and long-term effects of training and recovery on athlete wellness. For example, we might expect high-intensity training sessions to have immediate negative effects on wellness, but they could have positive impacts on wellness at longer time scales given proper recovery.

Under either the univariate or multivariate latent factor model, we can borrow strength across individuals by incorporating shared effects. Here, we include shared distributed lag coefficients. Recall that in (5) and (6), denotes the lagged coefficient for variable , individual , and lag . We model where is the global mean coefficient of covariate at lag and represents the variability across individuals for this effect. We can obtain inference with respect to these global parameters to provide insight into the general effects of the covariates at various lags as well as the measures of variability across individuals.

3.3 Identifiability constraints

We begin with a general depiction of the important identifiability constraints of the model parameters assuming one athlete (i.e., ). As such, we drop the dependence on in the following. Additionally, we note that there is more than just one set of parameter constraints that will result in an identifiable model, and will therefore justify our choices with regard to desired inference when necessary.

First, as is customary in probit regression models, the first threshold parameter, for each (Chib and Greenberg, 1998). This enables the identification of the intercept terms, . Then, to identify the lag coefficients, , for , and , without loss of generality, we set the factor coefficients for the first metric equal to 1 (Cagnone et al., 2009). In the univariate latent factor model, this results in and in the bivariate latent factor model, this results in . With the number of metrics , we also must specify a common in order to identify the metric-specific latent factor coefficients, .

Another common identifiability constraint for probit regression models imposes a fixed variance for the latent continuous metrics,

(Chib and Greenberg, 1998). This is the variance of from (2) which is denoted . With metric specific threshold parameters , we drop the dependence on such that . One option is to fix and model the variance parameters of the latent factors, , in the univariate latent factor model, and and in the bivariate latent factor model (Cagnone et al., 2009; Cagnone and Viroli, 2018). However, since we are modeling the latent factors using distributed lag models, this approach can mask some of the effects of the lagged covariates as well as the relationships between the latent factors and the ordinal wellness metrics. Therefore, we opt to work with the marginal variance of , which is equal to

in the univariate factor model and

in the bivariate factor model. For , this reduces to and . We set and and use a Dirichlet prior with two and three categories, respectively. Details regarding this prior are given below.

In extending to modeling multiple athletes, we add subscript to each of the parameters and latent factors. That is, we have individual specific threshold parameters, intercepts, factor coefficients, latent factors, and variances. In addition, we introduce the global mean coefficients, , and variances, . By imposing the same set of constraints above for each individual, the model parameters are identifiable.

Fitting this model to the referee ordinal wellness data discussed in Section 2 requires one additional modification. In looking at the ordinal response distributions (Figures 1 and 11), notice that for some individuals and some metrics, some ordinal values have few, if any, counts (e.g., Athlete A: Stress). In such a case, there is no information in the data to inform about the cut points for these individual and metric combinations. Therefore, we drop the individual specific threshold parameters to leverage information across athletes for each metric. With this modification, we can relax the constraint on to allow for metric specific thresholds,

. The shared metric-specific threshold approach across individuals is preferred over having individual threshold parameters that are shared across metrics for two reasons. First, some referees have very small counts for some ordinal values, even when aggregated across metrics, resulting in challenges in estimating these parameters. When aggregating across referees for a given metric, the distribution of observations across ordinal values is much more uniform. Second, by retaining the metric-specific threshold parameters, we can more easily compare the metric-specific relationships with the latent factor(s). That is, we can directly compute correlations between the latent factor(s) and the latent continuous wellness metrics as discussed below. Due to ordinal data not having an identifiable scale, specifying individual threshold parameters that are shared across metrics requires computing more complex functions of the model parameters in order to obtain this important inference.

3.4 Model inference and priors

Model inference was obtained in a Bayesian framework. Prior distributions are assigned to each model parameter and non-informative and conjugate priors were chosen when available. Each global mean lagged coefficient parameter is assigned an independent, conjugate hyper prior where

for all and . The variance parameters are assigned independent Inverse-Gamma(0.01, 0.01) priors. The latent factor coefficients are assigned independent normal priors, where for in the univariate factor model and in the bivariate factor model.

Given the identifiability constraints above for the variance parameters, we specify a two category Dirichlet prior for in the univariate latent factor model and a three category Dirichlet prior for in the bivariate model. Both Dirichlet priors are defined with concentration parameter 10 for each category.

The threshold parameters were modeled on a transformed scale due to their order restriction where . To ensure these inequalities hold true, with for all and , we define for and model . This transformation improves mixing and convergence when using MCMC for model inference (Higgs and Hoeting, 2010). Sampling the threshold parameters requires a Metropolis step within the MCMC algorithm.

3.5 Posterior inference

Important posterior inference includes estimates of the model parameters as well as correlation and relative importance measures for each metric. Dropping the dependence on for ease of notation, let define the correlation between latent ordinal response metric and the univariate latent factor , computed as


For the multivariate latent factor model, we can define analogous correlations for each factor where


These correlations provide a measure for which to compare the importance of each wellness metric in capturing the variation in the latent factor. To compare across metric, we compute the relative importance of each metric for the univariate latent factor model as


and as


for the multivariate factor model. Metrics with higher relative importance indicate that they are more important in capturing the variation in the latent factor. For the univariate latent factor model, these relative importance scores could be used as weights in computing an overall wellness score for each individual. Then, these latent factors could be monitored through time to identify possible changes in each athlete’s wellness in response to training and recovery throughout the season. For example, we can investigate the variation in the latent factors for each athlete by comparing match days to the days leading up to and following the match. Similarly, for the multivariate metric, these weights could identify which metrics are more or less important in explaining the variation in each particular latent factor, and again, can be monitored throughout the season to identify potential spikes in wellness as a response to training and recovery. We can obtain full posterior distributions, including estimates of uncertainty, of all correlations and relative importance metrics post model fitting.

4 Application: modeling athlete wellness

We apply our model to the subjective ordinal wellness data and workload and recovery data for 20 MLS referees collected during the 2015 and 2016 seasons. We investigated the effects of the workload and recovery variables on wellness for lags up to 10 day. Therefore, we limited the analysis to daily data for which at least 10 prior days of workload and recovery data were available. The number of observations for each individual ranged from 170 to 467 days.

The univariate and multivariate latent factor models were each fitted to the data. Model inference was obtained using Markov chain Monte Carlo and a hybrid Metropolis-within-Gibbs sampling algorithm. The chain was run for 100,000 iterations, and the first 20,000 were discarded as burn-in. Traceplots of the chain for each parameter were investigated for convergences and no issues were detected.

Boxplots of the posterior distributions of the global lagged coefficients of the univariate latent factor model are shown in Figure 4

for the workload and recovery covariates. Also shown are the upper and lower limits of the central 95% credible intervals. In general, workload is negatively related with athlete wellness, and the previous day’s workout (lag equal to 1) is the most significant. This implies that, in general, workload has an acute effect on player wellness, such that a heavy workload on the previous day tends to lead to a decrease in wellness on the following day. The lagged coefficients of the recovery variable show a positive relationship between recovery and wellness, where a large recovery value corresponds to high sleep quality and quantity. The lagged coefficients for this variable are significant for lags 1 through 5 as indicated by the 95% credible intervals not including 0. These significant lagged coefficients suggest that sleep quality and quantity may have a longer lasting effect on athlete wellness.

Posterior distributions of the individual-specific lagged coefficients are shown for Athlete A, B, C, and D in Figures 5 and 6 for the workload and recovery variables, respectively. In general, there is a lot of variation between individuals with regard to the lagged effects of the two variables. For example, Athlete A and B experience significant negative effects of workload at both lags 1 and 2, whereas Athlete C and D do not experience such negative effects. In fact, a heavy workload on the previous day has a positive relationship with wellness for Athlete D. For all four athletes, we see positive effects of workload at longer lags (e.g., lag 9 for Athlete A, lags 7-9 for Athlete B). The individual-specific lagged coefficients of the recovery variable show that the previous nights sleep quantity and quality have a very significant positive relationship with wellness for each athlete (Figure 6). However, we detect a more short-term effect of recovery for Athlete A and B (2 days) relative to Athlete C and D (3+ days) than the average shown in Figure 4.

The latent factors, , provide a univariate measure of wellness for each individual on each day. Figure 7 shows boxplots of the posterior mean estimates of for the four athletes for match days, denoted “M,” compared to the 3 days leading up to and following the match. Due to possible variation throughout the seasons, the estimates are centered within each match week by subtracting the 7-day average. Estimates of the latent factors vary throughout the 7-day period for each athlete. In general, the wellness of Athlete A is highest on match day relative to the days leading up to and following the match. The wellness estimates for Athlete B and C show less variation across days, although wellness for Athlete B is lower, on average, the day following the match relative to the match day. Wellness for Athlete D appears similar across all days except for the day following the match, in which wellness is higher.

We computed the correlation between the vectors of the latent continuous response metric, and the univariate latent factor, , for each athlete and metric. Boxplots of the posterior distributions for these correlations for the four athletes are shown in Figure 8, indicating variation both within metric across individuals and across metrics within individual. The majority of the significant correlations between the ordinal wellness metric and latent factor are positive, although the correlation was negative for Athlete B for the appetite metric.

To compare the significance of the different metrics within an individual in relation to the latent factor, we compute the relative importance statistics defined in (9). The relative importance statistics give a measure of the ability of each wellness metric at capturing the variation in the latent factor. The posterior mean estimates of for each of the four athletes across the six metrics are shown in Figure 9. The relative importance of each ordinal wellness metric varies across the athletes. Note that a value of 1/6 for each metric would correspond to an equal weighting. The most notable similarity between the four athletes is the high relative importance of energy, with each greater than 1/6. The relative importance of motivation and mood vary a lot between athletes. The estimates for Athlete C and D closely resemble an equal weighting scheme across the six metrics, whereas A and B each have unequal relative importance estimates with emphasis on mood, energy and tiredness for Athlete A, and motivation, energy, and tiredness for Athlete B. These results clearly depict a difference between computing the average across all metrics and the utility of the multivariate model in leveraging the individual wellness measures. Plots of the correlation and relative importance metrics for all 20 athletes are included in Figures A2 and A3 of the Supplementary Material.

The multivariate latent factor model resulted in similar global and individual lagged coefficient estimates as the univariate model. (See Figures A4 - A6 of the Supplementary Material). In addition, the variation in the estimates of the two latent factors across days leading up to and following each match also appeared similar to the univariate model (Figures A7 and A8). Important inference from the multivariate model consists of the factor-specific correlations and relative importance estimates for each wellness metric. Posterior distributions of the correlation estimates are shown in Figure 10 for both workload and recovery variables for the same four athletes, Athlete A, B, C, and D. (A similar figure with all athletes is given in Figure A9 of the Supplementary Material). This figure shows some important similarities and differences for each of the wellness metrics in terms of the correlations with the two latent factors. Both energy and tiredness show stronger correlations with recovery than workload for Athlete B, C, and D. The motivation metric is significantly more correlated with workload than recovery for Athlete A and B, whereas it is similar for Athlete C and D. Stress and mood both appear more strongly correlated with workload, whereas appetite appears more strongly correlated with recovery for each athlete.

Posterior mean estimates of relative importance for each latent factor for the four athletes are shown in Figure 11. Some wellness metrics that appeared insignificant in the univariate latent factor model now show significance when the workload and recovery latent factors are considered separately. For example, motivation appears significant for both workload and recovery for Athlete A whereas it was the least important metric in the univariate latent model. The relative importance of mood on the workload latent factor is greater than 1/6 for each athlete, and is the highest relative importance for Athlete A. Energy and tiredness appear to capture the majority of the variation in the recovery latent factor for Athlete B. Interestingly, Athlete C retains a fairly equal weighting scheme across the six metrics for both workload and recovery. The relative importance of stress is high for workload and low for recovery for Athlete D, whereas tiredness is low for workload and high for recovery. Additional comparisons between all athletes with respect to the relative importance for each metric and latent factor can be made looking at Figure A10 of the Supplementary Material. Interestingly, none of the metrics appear to be uniformly insignificant across the 20 athletes, providing justification in each component of the self-assessment survey.

5 Discussion

We develop a joint multivariate latent factor model to study the relationships between athlete wellness, training, and recovery using subjective and objective measures. Importantly, the multivariate response model incorporates the information from each of the subjective ordinal wellness variables for each individual. Additionally, the univariate and bivariate latent factors are modeled using distributed lag models to identify the short- and long-term effects of training and recovery. Individual-specific parameters enable individual-level inference with respect to these effects. The relative importance indices provide individual-specific estimates of the sensitivity of each ordinal wellness metric to the variation in training and recovery. The joint modeling approach enables the sharing of information across individuals to strengthen the results.

We applied our model to daily wellness, training, and recovery data collected across two MLS seasons. While the results show important similarities and differences across athletes with regard to the training and recovery effects on wellness and the importance of each of the wellness variables, the model could provide new and insightful information when applied to competing athletes. When referencing physical performance, our findings align with important known differences in training programs for referees and players. Training programs aim at maximizing the physical performance of players on match days, whereas referee training does not place the same significance on these days. This suggests interesting comparisons could be made using the results of this type of analysis between athletes who compete in different sports with differing levels of intensity and periods of recovery. For example, football has regular weekly game schedules at the college and professional levels, whereas soccer matches are scheduled typically twice per week, and hockey leagues often play two and three-game series with games on back-to-back nights. Maximizing player performance under these different competition schedules and levels of intensity using subjective wellness data is an open area for future work.

Distributed lag models, like those applied in this work, make the assumption that the lagged coefficients are constant in time. That is, the effect of variable at lag on the response at time is captured by and is the same for all . In terms of training and recovery for athletes, one could argue that these effects might vary throughout a season as a function of fitness and fatigue. For example, if an athlete’s fitness level is low during the early part of the season, a hard training session might have a longer lasting effect on wellness than it would mid-season when the athlete is at peak fitness. Alternatively, as the season wears on, an athlete might require a longer recovery time in order to return to their maximum athletic performance potential. As future work, we plan to incorporate time-varying parameters into the distributed lag models. This will require strategic model development in in order to minimize the number of additional parameters and retain computation efficiency in model fitting. The scope of this future work spans beyond sports, as the lagged effects of environmental processes could also have important time-varying features.


  • R. Akenhead and G. P. Nassis (2016) Training load and player monitoring in high-level football: current practice and perceptions. International Journal of Sports Physiology and Performance 11 (5), pp. 587–593. Cited by: §1, §1.
  • J. H. Albert and S. Chib (1993) Bayesian analysis of binary and polychotomous response data. Journal of the American Statistical Association 88 (422), pp. 669–679. Cited by: §3.1.
  • J. Borresen and M. I. Lambert (2008) Quantifying training load: a comparison of subjective and objective methods. International Journal of Sports Physiology and Performance 3 (1), pp. 16–30. Cited by: §1.
  • J. Borresen and M. I. Lambert (2009) The quantification of training load, the training response and the effect on performance. Sports Medicine 39 (9), pp. 779–795. Cited by: §1.
  • P. C. Bourdon, M. Cardinale, A. Murray, P. Gastin, M. Kellmann, M. C. Varley, T. J. Gabbett, A. J. Coutts, D. J. Burgess, W. Gregson, et al. (2017) Monitoring athlete training loads: consensus statement. International Journal of Sports Physiology and Performance 12 (s2), pp. S2–161. Cited by: §1, §1.
  • M. S. Brink, E. Nederhof, C. Visscher, S. L. Schmikli, and K. A. Lemmink (2010) Monitoring load, recovery, and performance in young elite soccer players. The Journal of Strength & Conditioning Research 24 (3), pp. 597–603. Cited by: §2.
  • M. Buchheit, S. Racinais, J. Bilsborough, P. Bourdon, S. Voss, J. Hocking, J. Cordy, A. Mendez-Villanueva, and A. Coutts (2013) Monitoring fitness, fatigue and running performance during a pre-season training camp in elite football players. Journal of Science and Medicine in Sport 16 (6), pp. 550–555. Cited by: §1.
  • S. Cagnone, I. Moustaki, and V. Vasdekis (2009) Latent variable models for multivariate longitudinal ordinal responses. British Journal of Mathematical and Statistical Psychology 62 (2), pp. 401–415. Cited by: §1, §3.3, §3.3.
  • S. Cagnone and C. Viroli (2018) Multivariate latent variable transition models of longitudinal mixed data: an analysis on alcohol use disorder. Journal of the Royal Statistical Society: Series C (Applied Statistics) 67 (5), pp. 1399–1418. Cited by: §1, §3.3.
  • F. Chaubert, F. Mortier, and L. Saint André (2008) Multivariate dynamic model for ordinal outcomes.

    Journal of Multivariate Analysis

    99 (8), pp. 1717–1732.
    Cited by: §1.
  • S. Chib and E. Greenberg (1998) Analysis of multivariate probit models. Biometrika 85 (2), pp. 347–361. Cited by: §3.3, §3.3.
  • V. De Silva, M. Caine, J. Skinner, S. Dogan, A. Kondoz, T. Peter, E. Axtell, M. Birnie, and B. Smith (2018) Player tracking data analytics as a tool for physical performance management in football: a case study from Chelsea Football Club Academy. Sports 6 (4), pp. 130. Cited by: §1.
  • M. DeYoreo and A. Kottas (2018) Bayesian nonparametric modeling for multivariate ordinal regression. Journal of Computational and Graphical Statistics 27 (1), pp. 71–84. Cited by: §1.
  • C. Foster, E. Daines, L. Hector, A. C. Snyder, and R. Welsh (1996) Athletic performance in relation to training load.. Wisconsin Medical Journal 95 (6), pp. 370–374. Cited by: §2.
  • C. Foster, J. A. Florhaug, J. Franklin, L. Gottschall, L. A. Hrovatin, S. Parker, P. Doleshal, and C. Dodge (2001) A new approach to monitoring exercise training. The Journal of Strength & Conditioning Research 15 (1), pp. 109–115. Cited by: §2, §2.
  • T. F. Gallo, S. J. Cormack, T. J. Gabbett, and C. H. Lorenzen (2017) Self-reported wellness profiles of professional Australian football players during the competition phase of the season. The Journal of Strength & Conditioning Research 31 (2), pp. 495–502. Cited by: §1.
  • A. Gasparrini, B. Armstrong, and M. G. Kenward (2010) Distributed lag non-linear models. Statistics in Medicine 29 (21), pp. 2224–2234. Cited by: §1.
  • M. Haddad, G. Stylianides, L. Djaoui, A. Dellal, and K. Chamari (2017) Session-rpe method for training load monitoring: validity, ecological usefulness, and influencing factors. Frontiers in Neuroscience 11, pp. 612. Cited by: §2.
  • L. D. Haugh and G. E. Box (1977) Identification of dynamic regression (distributed lag) models connecting two time series. Journal of the American Statistical Association 72 (357), pp. 121–130. Cited by: §1.
  • M. D. Higgs and J. A. Hoeting (2010) A clipped latent variable model for spatially correlated ordered categorical data. Computational Statistics & Data Analysis 54 (8), pp. 1999–2011. Cited by: §3.4.
  • R. Hirk, K. Hornik, L. Vana, and A. Genz (2019) Mvord: Multivariate Ordinal Regression Models. Note: R package version 0.3.6 External Links: Link Cited by: §1.
  • R. J. Hyndman and G. Athanasopoulos (2018) Forecasting: principles and practice. OTexts. Cited by: §1.
  • M. S. Itter, J. Vanhatalo, and A. O. Finley (2019) EcoMem: An R package for quantifying ecological memory. Environmental Modelling & Software. Cited by: §1.
  • L. C. Liu and D. Hedeker (2006) A mixed-effects regression model for longitudinal multivariate ordinal data. Biometrics 62 (1), pp. 261–268. Cited by: §1.
  • F. M. Lord (2012) Applications of item response theory to practical testing problems. Routledge. Cited by: §1.
  • R. Meeusen, M. Duclos, C. Foster, A. Fry, M. Gleeson, D. Nieman, J. Raglin, G. Rietjens, J. Steinacker, and A. Urhausen (2013) Prevention, diagnosis, and treatment of the overtraining syndrome: joint consensus statement of the European College of Sport Science and the American College of Sports Medicine.. Medicine and Science in Sports and Exercise 45 (1), pp. 186. Cited by: §1.
  • I. Mujika (2017) Quantification of training and competition loads in endurance sports: methods and applications. International Journal of Sports Physiology and Performance 12 (s2), pp. S2–9. Cited by: §1.
  • K. Ogle, J. J. Barber, G. A. Barron-Gafford, L. P. Bentley, J. M. Young, T. E. Huxman, M. E. Loik, and D. T. Tissue (2015) Quantifying ecological memory in plant and ecosystem processes. Ecology Letters 18 (3), pp. 221–235. Cited by: §1, §3.2.
  • A. E. Saw, L. C. Main, and P. B. Gastin (2016) Monitoring the athlete training response: subjective self-reported measures trump commonly used objective measures: a systematic review. British Journal of Sports Medicine 50 (5), pp. 281–291. Cited by: §1, §1.
  • E. M. Schliep and J. A. Hoeting (2013) Multilevel latent Gaussian process model for mixed discrete and continuous multivariate response data. Journal of Agricultural, Biological, and Environmental Statistics 18 (4), pp. 492–513. Cited by: §1.
  • J. Schwartz (2000) The distributed lag between air pollution and daily deaths. Epidemiology 11 (3), pp. 320–326. Cited by: §1.
  • A. Skrondal and S. Rabe-Hesketh (2004) Generalized latent variable modeling: multilevel, longitudinal, and structural equation models. Chapman and Hall/CRC. Cited by: §1.
  • F. Tavares, P. Healey, T. B. Smith, and M. Driller (2018) Short-term effect of training and competition on muscle soreness and neuromuscular performance in elite rugby athletes. Australian Journal of Strength and Conditioning 26 (1), pp. 11–7. Cited by: §1.
  • R. T. Thorpe, A. J. Strudwick, M. Buchheit, G. Atkinson, B. Drust, and W. Gregson (2015) Monitoring fatigue during the in-season competitive phase in elite soccer players. International Journal of Sports Physiology and Performance 10 (8), pp. 958–964. Cited by: §1.
  • R. T. Thorpe, A. J. Strudwick, M. Buchheit, G. Atkinson, B. Drust, and W. Gregson (2017) The influence of changes in acute training load on daily sensitivity of morning-measured fatigue variables in elite soccer players. International Journal of Sports Physiology and Performance 12 (s2), pp. S2–107. Cited by: §1.

Appendix A Supplementary Material

Figure A1: Raw summaries of the ordinal wellness metrics energy, motivation, mood, and appetite for four athletes.
Figure A2: Boxplots of the posterior distributions of the correlation between the univariate latent factor, and , for all athletes and metric.
Figure A3: Posterior mean estimates of the relative importance statistics, , defined in (9) for all athletes and metric.
Figure A4: Distribution of the global lagged coefficients for the two factor model for workload (left) and recovery (right). indicates 95% credible interval.
Figure A5: Distribution of the individual specific lagged coefficients for the workload latent factor. indicates 95% credible interval.
Figure A6: Distribution of the individual specific lagged coefficients for the recovery latent factor. indicates 95% credible interval.
Figure A7: Posterior mean estimates of the workload latent factor, , for each individual and match (M), as well as the three days leading up to and following the match. The estimates shown are mean-centered for each match.
Figure A8: Posterior mean estimates of the recovery latent factor, , for each individual and match (M), as well as the three days leading up to and following the match. The estimates shown are mean-centered for each match.
Figure A9: Boxplots of the posterior distributions of the correlation between each latent continuous wellness metric, , and the latent factors for workload and recovery and , for each athlete and latent metric variable .
Figure A10: Posterior mean estimates of the relative importance statistics, , defined in (10) for each athlete and metric for the workload latent factor (top) and recovery latent factor (bottom).