Flexible modelling of time-varying exposures and recurrent events to analyse training load effects in team sports injuries

Abstract

We present a flexible modelling approach to analyse time-varying exposures and recurrent events in team sports injuries. The approach is based on the piece-wise exponential additive mixed model where the effects of past exposures (i.e. high-intensity training loads) may accumulate over time and present complex forms of association. In order to identify a relevant time window at which past exposures have an impact on the current risk, we propose a penalty approach. We conduct a simulation study to evaluate the performance of the proposed model, under different true weight functions and different levels of heterogeneity between recurrent events. Finally, we illustrate the approach with a case study application involving an elite male football team participating in the Spanish LaLiga competition. The cohort includes time-loss injuries and external training load variables tracked by Global Positioning System devices, during the seasons 2017–2018 and 2018–2019.

1 Introduction

In recent years, research on sports injuries has received much attention from a wide range of areas, including the statistical modelling field, given the ever-increasing amount of data now being collected (Casals & Finch, 2017; Sainani et al., 2021). Injuries are common, and modelling and understanding their occurrence would assist in developing tailored prevention programs. In this regard, it has been claimed that analysing changes in training load exposures over time is crucial when studying sports injury aetiology, as well as analysing recurrent injury, subsequent injury, or injury exacerbation (Nielsen et al., 2020). Athletes are repeatedly exposed to high competition demands that, in turn, increase the strain on their bodies and exposure to the risk of injury. Injury prevention depends on the athlete’s capacity to tolerate repeated exposures to injury risk.

The physical load from training and competition, hereafter termed as training load, is defined as ‘the cumulative stress placed on an individual from multiple training sessions and games over a period of time’ (Gabbett et al., 2014). Training load can thus be applied to the athlete over varying time periods and with varying magnitudes (Soligard et al., 2016). Nowadays, vast amounts of data are gathered (e.g. training and competition time, distance, speed, power output), primarily using Global Positioning System (GPS) devices, to quantify training load. Indeed, the study of training load is key to developing effective training plan strategies that enhance athletes’ performance while also reducing their risk of injury. The relationship between training load and injury, however, remains uncertain (Griffin et al., 2020; Windt et al., 2018).

Quantifying this relationship requires the development of an etiologically plausible time-varying exposure model, which estimates how previous training affects the injury risk. The effects of past exposures may cumulate over time and exhibit complex forms of association. Moreover, the model should account for the possibility that subsequent injuries may be associated within players. To take into account the latter, specifically dependencies induced by injury recurrence, as well as the intensity and duration of past exposures, we propose a piece-wise exponential additive mixed model (PAMM) (Bender et al., 2018) with weighted cumulative exposure (WCE) cumulative effects (Sylvestre & Abrahamowicz, 2009).

PAMMs are a semi-parametric extension of the piece-wise exponential model (M. Friedman, 1982; Holford, 1980; Laird & Olivier, 1981) that allow for penalized estimation of flexible survival models (cf. Argyropoulos & Unruh, 2015; Bender et al., 2018 for a thorough overview) with a wide range of covariate effects, such as nonlinear, time-varying effects, cumulative effects, and/or random effects. This framework has also been shown to support the estimation of cumulative effects of time-varying exposure histories. The WCE-type cumulative effect suggested by Sylvestre and Abrahamowicz (2009), a weighted sum of all past exposures over a relevant time window, is a common way to address this. The weight function assigns weights to past exposures based on the time elapsed since the exposure occurred, which, ideally, is determined according to the true underlying biological mechanism. They proposed to estimate the weight function using B-spline regression.

The PAMM methodology has been applied in many recent publications. Bender et al. (2019) explore complex exposure-lag-response associations and provide a general formulation of PAMMs that includes previous approaches for cumulative effects like the WCE model and the distributed lag nonlinear model (DLNM) (Gasparrini et al., 2017) as special cases. Ramjith et al. (2022) study the PAMM framework for recurrent events analysis and show that under the assumption of proportional hazards, PAMM and the shared frailty Cox model (McGilchrist & Aisbett, 1991) are equivalent. Danieli and Abrahamowicz (2019) and Li et al. (2022) introduce approaches to model cumulative effects of time-varying exposures with competing risks, via cause-specific hazards model and subdistribution hazards model for each competing event, respectively, by incorporating separate Cox WCE models with an event-specific weight function. Recently, in the field of sports medicine, Bache-Mathiesen et al. (2022) evaluated different methods to assess the cumulative effect of training load on the risk of injury in team sports and suggest the use of DLNM.

In this work, we extend the PAMM by incorporating WCE-type cumulative effects in the recurrent events setting combined with a method to identify a relevant time window in which past exposures have an effect. We further demonstrate the practical application of this model in the field of sports medicine.

The proposed modelling framework is detailed throughout Section 2, which first introduces the PAMM framework and then focuses on how we adapt it to flexibly model time-varying exposures and recurrent events in addition to how we penalize the model for identifying a relevant window. Section 3 illustrates a real-life application of this method to assess the cumulative effects of past training exposures on the hazard of subsequent injuries in a football team. Section 4 describes the simulation study carried out to evaluate the proposed models performance and examine their properties. The final section concludes the paper with a discussion.

2 Methods

PAMMs transform a survival task to a Poisson regression task by partitioning the follow-up period into a finite number of intervals and assuming that hazards are piece-wise constant in each of these intervals. The key idea—explained in detail in the following sections—is to flexibly model the baseline hazard and other time-varying effects using penalized splines. This approach alleviates the problem of an arbitrary selection of the cut-points, that partition the follow-up time, and thus avoids over- and underfitting as well as instability issues (Bender et al., 2018). In practice, one can simply use a relatively large number of cut-points and use spline basis functions evaluated at e.g. $κj$⁠, the right end of each interval, and penalize the wiggliness of the estimate via penalized splines (cf. Bender, 2018, for an empirical discussion on the placing of cut-points).

We denote for player l, where $l=1,…,L$ and L the total number of players, $Yil$ as the ith recurrent time, where $i=1,…,nl$ and $nl$ the maximum event number that individual l has been at risk for, and $Cl$ as the censoring time. Then, $Til=minYil,Cl$ corresponds to each follow-up time and $δil$ is a binary indicator for recurrent events (i.e. subsequent injuries) which is 0 if the player has not been injured and 1 if the ith injury time, $Yil$⁠, is observed. That is, $δil=1(Til=Yil)$⁠, where $1$ represents the indicator function. In the following, we will use t instead of $til$⁠, for ease of notation.

Therefore, as we proceed with the modelling of injuries, the hazard rate of the ith injury (event) of the lth player, given the player’s training exposure history $zl(t)={zl(tz):tz≤t}$⁠, is expressed as:

for all $t∈(κj−1,κj],t>0$ and $tj:=κj$⁠, where $κj$ are the $J+1$ cut points defining J intervals that partition the study follow-up time $(0,tmax]$⁠, specifically, $κ0=0$⁠, $κJ=tmax$ and $j=1,…,J−1$⁠.

In equation (1), the expression $β0+f0(tj)$ denotes the log-baseline hazard, where $f0(tj)$ is expressed as a smooth term of the form $∑m=1M′γ0mBm(tj)$⁠. The term $g(zl(t),t)=∫τ(t)h(t,tz,zl(tz))dtz$ denotes the cumulative effect of $zl(t)$ at time t, i.e. the past exposure effects of $zl(t)$ cumulating over a relevant time-window $τ(t)$⁠, resulting in a sum of weighted effects. The dependence induced by subsequent injuries is accounted for by $bl$⁠, a Gaussian random effect (i.e. frailty) associated with player l, which acts as a random intercept term for the lth player, i.e. $bl∼N(0,σb2)$⁠.

For a WCE-type effect, in equation (1), we consider time-varying exposure effects weighted by latency $t−tz$ and linear in $z(tz)$⁠. That is, the contribution of covariate $z(t)$ observed at time $tz$ with value $z(tz)$⁠, is defined by $h(t,tz,z(tz)):=h(t−tz)z(tz)$⁠, and called partial effect. Thus, the cumulative effect $g(z(t),t)$ at follow-up time t is the integral of these partial effects over exposure times $tz$ contained within the so-called lag-lead window, $τ(t)$⁠, which controls how many observations of z contribute to the cumulative effect at time t (the minimal requirement is that $tz≤t$⁠).

Let $z(t)={z(tz):tz≤t}={z(tz,1),…,$  $z(tz,Q)}$ be the set of all registered exposure variables up to time t. Then, $g(z(t),t)$ is estimated with penalized splines (e.g. by using P-splines Eilers & Marx, 1996, 2021 that penalize the differences of neighbouring basis coefficients) and with quadrature weights $Δq=tz,q−tz,q+1$ (and $tz,0=t$⁠), the time difference between two consecutive exposure measurements, for numerical integration, as follows:

with $Δ~q=z(tz)(tz,q−tz,q+1)iftz,q∈τ(t)$ and 0 otherwise; $Bm(⋅)$ B-spline basis functions and $γm$ the associated spline coefficients.

2.1 Penalization of the weight function

One of the challenging issues is to determine a relevant time window $τ(t)$⁠. Without solid prior knowledge, it can be defined as $τ(t)={tz:t≥tz}$⁠, so all past exposures, collected before actual time t, contribute to the cumulative effect $g(z(t),t)$⁠. Yet, it is plausible that the effects of exposure variables may not be everlasting. As time passes, the effect of exposures recorded long ago may smoothly decrease to zero and eventually disappear. The exact length of the window, however, is usually unknown.

Subsequently, adapting the method by Obermeier et al. (2015) to the PAMM framework, we present two approaches to penalize the weight function, allowing it to transition smoothly to zero at the right end of the support interval: (i) a constrained-effect approach and (ii) a ridge penalty approach.

2.1.1 Constrained-effect approach

In this approach, we force the weight function, and its first derivative, to reach the zero value at time $t−tz,Q$⁠, by imposing the last two coefficients in equation (2) to be equal to zero (Sylvestre & Abrahamowicz, 2009). Working in matrix notation, the right part of equation (2) is given by:

with

denoting a $Q×M$-dimensional basis matrix of B-splines of degree d, $γ$ denoting a $M×1$ column-vector of associated spline coefficients and $Δ~$ denoting a $Q×1$ column-vector having $Δ~q$ as elements for $q=1,…,Q$⁠.

Therefore, the value $z(tz,Q)$ is assumed to have no impact on the current risk at t, by constraining the two last spline coefficients to zero, i.e. $γ~Q−1=γ~Q=0$⁠.

2.1.2 Ridge-penalty approach

This approach consists of adding a shrinkage L2-penalty to penalize the last B-spline basis coefficients of $γ$ in equation (3), see Obermeier et al. (2015).

Hence, we seek that the last coefficient $γ~Q$ to be close to zero:

At any chosen time point $tz,q∈τ(t)$⁠, exactly $d+1$  d-degree B-spline basis functions are nonzero. Thus, equation (4) can be expressed as,

Notice also that we implicitly assume that $γ~Q+1$ is zero. Following, the last $γ~Q$ coefficient can be forced to shrink towards zero by penalizing the last $d+1$  γ-coefficients by an $M×M$ shrinkage matrix $Kr$⁠, a diagonal matrix with all elements equal to zero except the last $d+1$ which have value one, i.e. $Kr=diag(0M−(d+1)×1,1(d+1)×1)$⁠.

Then, an extra regularization parameter controls the shrinking of the last basis coefficients. Large values of this parameter imply strong shrinkage of the last $γ~Q$ coefficient and decrease the estimated values of it.

2.2 Estimation and inference

The estimation of the model coefficients $γ$ can be carried out by maximizing the penalized likelihood (Wood, 2011). In this case, for the Poisson GAMM, the model deviance $l(γ)=∑l=1L∑il=1nl{δillog(λil(t|zl(t),bl))−λil(t|zl(t),bl)}$⁠, with $δil∈{0,1}$ the ith event indicator of subject l, is penalized with a smoothing and a shrinkage matrix, $Kd$ and $Kr$⁠, giving rise to:

For the smoothing penalty matrix, we use second-order differences, i.e. $Kd=D2′D2$⁠, where $D2$ is the matrix representation of applying the Δ operator to α twice: $Δ2α=Δ(αj−αj−1)=αj−2αj−1+αj−2$⁠. The smoothing parameter $λd$ penalizes large differences in adjacent basis coefficients, while the regularization parameter $λr$ shrinks the last basis coefficient.

Therefore, the coefficients can be estimated via penalized iteratively reweighted least squares P-IRLS (Marx & Eilers, 1998; Wood, 2017). The P-IRLS consists of iteratively updating the coefficient estimates until convergence is reached using numerical optimization methods of the restricted maximum likelihood. This is implemented in the gam function from mgcv (Wood, 2017) R package.

2.3 Software specification

The analyses of the case study and the simulation study are coded in R version 4.2.2, on a 64-bit Unix platform (x86$_$64 linux-gnu) computer, as well as on a high-performance cluster system. The package msm 1.6.9 was used to draw piece-wise exponential survival times, pammtools 0.5.8 and mgcv 1.8–41 to fit the models, batchtools 0.9.16 to structure, write down and submit the simulation experiment in a convenient and reproducible fashion and the package injurytools 1.0.1 to structure and explore the football injury data set. The code to reproduce these analyses is available at https://github.com/lzumeta/flex-mod-training-loads-recu-injuries.

3 Application: football injury data

3.1 Data

We apply the proposed model to observational injury data from an elite male football team that competed in LaLiga during the 2017–2018 and 2018–2019 seasons. A total of $L=36$ players were monitored to assess their performance and health status through external training load variables (e.g. training and competition time, distance covered, speed, and heart rate, see Soligard et al., 2016) collected using tracking devices during each match and training session. These variables measure the physical exertion to which the players were exposed. A total of 72 noncontact time-loss injuries occurred among 23 players (64%) and 15 (65%, 15/23) players were reinjured during the follow-up, see Figure 1.

Our focus is on understanding the relationship between external training load and time-loss injuries (Fuller et al., 2006), particularly how the cumulative stress from multiple training sessions and matches over time impacts a player’s risk of sustaining a (subsequent) injury. We specifically focus on the variables average speed per session (Speed) and total distance covered (Dist) per session, as external training load variables that represent the intensity of each session. These metrics are chosen for their clear link to physical exertion and injury risk. Average speed reflects the session’s sustained intensity, while total distance indicates overall workload. Both are crucial in evaluating the cumulative stress imparted on players, potentially leading to injuries when consistently high.

3.2 Why fit a flexible time-varying exposures model?

In practical settings, a common question that practitioners face is ‘how much training load is too much among players with different characteristics?’ (Nielsen et al., 2018; Soligard et al., 2016). To address this question, researchers have studied changes in training load as a primary exposure to injury. For example, metrics such as week-to-week changes and the acute:chronic workload ratio (ACWR) have been heavily relied upon. This latter measure compares short-term (acute) training load to longer-term (chronic) load to estimate injury risk (see Section A.2 of the online supplementary material for a more comprehensive explanation). Numerous recent studies have investigated the relationship between training load and injury risk (Bache-Mathiesen et al., 2021, 2022; Carey et al., 2018; Hulin et al., 2016; Windt & Gabbett, 2017).

Before describing our modelling approach, we may first ask why a flexible time-varying exposure model might be preferred over relying on predetermined measures for calculating changes in training load. First, predefined metrics such as week-to-week changes, ACWR and its variants, or cumulative load measures, tend to simplify the complex relationships between training load and injury risk and may miss important time-varying patterns or immediate changes in load. Second, the effects of long-term exposure to training load may manifest days or even weeks later, but this time period is neither fixed nor previously known. Therefore, a method that can identify this time window for each specific context is highly valuable. Third, flexible models can capture nonlinear patterns learning directly from the data. For these reasons, we adopt a flexible time-varying exposure model with penalization to simultaneously account for potentially nonlinear, time-varying, and cumulative effects while identifying relevant time windows.

3.3 Model specification

We consider that the unit of the follow-up time t, as well as of the exposure time $tz$⁠, to be the nth number of session (i.e. match and training sessions). The analysis time zero is defined as the first session (match or training session), from 7 July 2017, in which the player has taken part in the team. Players are followed until an injury occurs, or until they are transferred to another team, the end of the contract or the end of the study (18 May 2019), whichever occurs first. Players fully recovered from an injury are followed again until one of the previously described events occurs, see Figure 1. We consider that external training load is applied to the player over varying time periods and with varying magnitude by considering cumulative effects and we adjust for the type of session, whether training or match session, since it has been suggested as one of the primary risk factors (Bahr, 2003; Bahr et al., 2020); and account for subsequent injuries adding a random effect (Gaussian frailty). Thus, the log-hazard rate of player l of the fitted model is expressed as:

where $β0+f0(tj)$ indicates the log-baseline hazard rate, $zltypesession(tj)$ the type of session undertaken by player l at $tj$⁠, $g1$ and $g2$ are nonlinear time-varying effects of the training load variables and $bl$ a Gaussian random intercept term associated to player l. The cumulative effects, $g1$ and $g2$⁠, are defined as $∫τSpeed(t)h(t−tz)zlSpeed(tz)dtz$ and $∫τDist(t)h(t−tz)zlDist(tz)dtz$⁠, and each lag-lead window, $τHRt(t)$ and $τHSR(t)$⁠, is chosen to be large enough to identify relevant past exposure effects by fitting a PAMM with a ridge penalization. Importantly, we add a minimum lag time of one session to minimize confounding by indication bias (Signorello et al., 2022), i.e. we exclude the current session to have an effect on the hazard of injury, $t>tz$⁠. The reason is that each session depends on the player’s physical condition and, presumably, sessions in which a player was injured had lower intensity in contrast to sessions in which he had no complaints. All smooth terms are estimated using P-Splines with second-order difference penalties.

3.4 Results

Table 1 shows the descriptive characteristics of the data, overall and by session type. The injury incidence was 4 injuries per 1,000 player-hours, and the injury burden was 90 days lost per 1,000 player-hours. The first calculates the rate at which new injury occurs (likelihood), whereas the second indicates how severe an injury is (consequences) (Bahr et al., 2020). Match injury incidence and burden were higher than those during training, 2.6 vs. 1.5 injuries per 1,000 player-hours and 60 vs. 33 days lost per 1,000 player-hours, respectively.

For the PAMM with a ridge penalization modelling approach, a sufficiently large lag-lead window must be set. We selected this window based on a reasonable range suggested by domain experts. Consequently, the estimated cumulative effects are computed considering that all Speed and Dist values recorded in the last 10 sessions prior to t (i.e. approximately three weeks before)) could have an effect on the hazard of injury at time t, i.e. $τ1(t)=τ2(t)={tz:t>tz∧t<tz+11}$⁠. The estimated partial effects for the Speed and Dist training load variables, $h^1(t−tz)z1(tz)$ and $h^2(t−tz)z2(tz)$⁠, are shown in Figure 2.

The results suggest that no more than seven sessions in the past are relevant for the cumulative effects of the Speed and Dist variables. Both cumulative effects are estimated to have a nonlinear decaying effect on the covariate $z$ with respect to latency and a linear effect on latency with respect to the covariate $z$⁠. Regarding the estimated partial effects of both variables, the values contributing the most to the hazard are those most recently recorded, while the contribution of the values recorded longer ago diminishes. Although there is not much difference in the trend, the greater the average speed and the total distance covered in recent sessions, the greater the impact on the resulting cumulative effect. See also Figure S1 of the online supplementary material.

The Gaussian frailty term (random intercepts) that accounts for the correlation between subsequent injuries from the same player is statistically significant (⁠$p-value<0.01$⁠), with $σ^b2=0.22$ the estimated variance. Players who suffered more injuries (e.g. Id04, Id28) have a higher baseline hazard of injury, as seen in Figure 3, which shows the estimated smooth baseline hazard, $λ^0(t)$⁠, together with the estimated player-specific smooth baseline hazard. With respect to the session type effect, match sessions have a higher risk of injury compared to training sessions (⁠$β^1=2.45$⁠, and $p-value<0.01$⁠). See also Table S1 and Figure S3 of the online supplementary material where the estimated linear and nonlinear effects are presented.

Lastly, we evaluated different methods for capturing cumulative effects by fitting several PAMM models. Each PAMM model incorporated a different approach commonly used in the sports medicine and exercise literature for estimating the training load effect: the ACWR using rolling averages, ACWR with exponentially weighted moving averages (EWMA), and an unweighted sum of the last seven days of exposure (cf. Section A.2 of the online supplementary material). Table 2 shows the comparison of the model’s performance based on likelihood-based measures. Our model shows a better overall performance compared to the other alternatives.

4 Simulation study

We conduct extensive simulation studies to evaluate the proposed models and to investigate their properties, aiming to validate their applicability in real-world settings. In particular, we aim to (i) assess the ability of the model to simultaneously estimate both, flexible WCE-type effects and heterogeneity resulting from recurrent events; and (ii) study the implementation of penalties on the basis coefficients to select the maximum length of the time window in which past exposures are cumulatively associated with the hazard.

4.1 Data generation

We draw survival times from the piece-wise exponential distribution. Let $nl$ be the maximum event number that individual l has been at risk for, where $l=1,…,L$⁠. Then, it suffices to specify a vector of piece-wise constant hazards $λ=(λ11,λ21,…,λn1,λ12,λ22,…,λn2,…,λ1L,λ2L,…,λnL)$⁠, in intervals defined by $J+1$ cut-points, i.e. by the vector of interval borders $κ=(0=κ0,…,κJ=tmax)$⁠. That is, $λil$ is composed of $(λil1,λil2,…,λilJ)$⁠, where each element $λilj$ is the hazard rate of ith event for individual l in the interval j, $il=1,…,nl$ and $j=1,…,J$⁠; and can be defined through a function of time t, current and past exposure covariates $z(t)$ and a random effect $bl$⁠, i.e. $λil,j(t|z(t),b)=f(t,z(t),b)=exp(const.+f0(t)+∫{tz:t≥tz}h(t−tz)z(tz)dt+bl)$⁠, evaluated at time $t=κj$⁠.

Then, we draw recurrent survival times from the piece-wise exponential distribution (PEXP), $t∼PEXP(λ,κ)$⁠, for which the algorithm is outlined in Table S3 of the online supplementary material. The hazard rate vector $λ$ is defined based on the simulation settings described in the following section. All further details on data generation are provided in Section B of the online supplementary material.

4.2 Scenarios and parameter settings

Basing on real-world application data for parameter settings, we simulate $Nsim=500$ times a cohort of L players, varying L in ${20,40,100}$⁠, with exposures recorded at $tz,1=1,tz,2=2,…,tz,Q=40=40$ days before the time at which we model the hazard, $zl(t)=(zl(tz,1),zl(tz,2)$⁠, $…,zl(tz,Q))$⁠, and draw survival times from the piece-wise exponential distribution under four different true weight functions, $h(t−tz)$⁠, (a) exponential decay, (b) bi-linear (c) early peak and (d) inverted U shapes, each defined over a $[0,tz,Q]$ interval (see the black curves in Figure 4); and under three different levels of heterogeneity between recurrent events, $σb∈{0.05,0.5,1}$⁠, indicating very low heterogeneity, low heterogeneity and high heterogeneity, respectively (see also Figure S4 in the online supplementary material). We then fit three different PAMMs with WCE-type cumulative effects: a model with no constraint (Uncons.), adding a constraint (Constr.) and adding a ridge penalty (Ridge). The performance of the models is evaluated by graphical inspection of the estimated $h^(t−tz)$ function in comparison to the true simulated $h(t−tz)$ function; the accuracy of these $h^(t−tz)$ WCE-type cumulative effects estimates are also evaluated via the mean RMSE, i.e. $RMSE¯$⁠, over all simulation runs, as:

where $Ntz=41$⁠, since $t−tz={0,1,2,…,40}$ takes 40 $+$1 number of different values.

We assess the accuracy of the standard deviation of the random effects, $σ^b$⁠, through:

We also evaluate the rate at which the estimated confidence interval of WCE-type cumulative effects estimates contains the true estimand $h(t−tz)$⁠, computing the mean coverage at the $1−α$ confidence level, i.e. $Coverage¯α$⁠, as:

where $ζq$ is the q-quantile of the standard normal distribution and $σ^h^$ the standard error of the estimated $h^(t−tz)$⁠.

4.3 Simulation results

Figure 5 shows boxplots of the distribution of the RMSE, the 95% point-wise coverage (across all time points) and the squared error of $σb$ across all simulation settings. In general, as the sample size (i.e. number of players) increases, the estimation errors for $h(t−tz)$ and $σb$ decrease and the mean point-wise coverage for $h(t−tz)$ increases, across all models. The Ridge model provides the most accurate estimate for $h(t−tz)$⁠, outperforming the other models considered (see the first row in Figure 5). On the other hand, the estimation of the standard deviation of the random effect, $σb$⁠, does not depend on the estimation of the shape of the weight function $h(t−tz)$⁠. All models provide similar estimates for $σb$⁠, but these estimates become less accurate for higher true $σb$ values (see the second row in Figure 5). The 95% coverage of $h(t−tz)$ for shapes (a), (c) and (d) shows underfitting, specifically for models that penalize the weight function. This underfitting occurs because the true wight function in these cases has features (e.g. sharp peaks) that are not well captured by the penalized models. In addition, the method used to calculate the coverage, i.e. point-wise, may also contribute to this observed underfitting (see the third row in Figure 5). The numerical results for all simulation settings are presented in Tables S4–S7 of the online supplementary material.

Results regarding the estimation of the weight function $h^(t−tz)$ for true weight functions (a)–(d), $L=100$ players and $σb=0.5$ are shown in Figure 4 (the rest of the scenarios are shown in Figures S8–S19 of the online supplementary material). In general, the model estimates capture well the underlying weight function. The models with an additional penalty (middle and bottom panels), Constr. and Ridge models, are more accurate for scenarios where the exposures that occurred relatively long ago (e.g. from the 20th lag on) have little impact on the risk at the actual time t, i.e. shapes (b) and (c). This is proved by the lower mean RMSE values in the Tables S4–S7 of the online supplementary material. The best estimation of $h(t−tz)$⁠, among the settings considered, is obtained for shape (b) bi-linear with ridge penalty, according to RMSE and 95% coverage, that has a mean RMSE of 0.010 and a mean coverage of 93.2% (see Tables S4–S7 of the online supplementary material).

5 Discussion

We extended and assessed the PAMM model class to the context of recurrent events with time-varying covariates, modelled as WCE-type cumulative effects. By introducing a ridge penalty to diminish the influence of past exposures registered long ago, we presented a method to determine a relevant time window based on data. Lastly, motivated by the research question regarding the association between external training load and (subsequent) time-loss injuries, we applied the proposed methodology in the sports medicine context.

Simulations indicate that PAMM with ridge penalization is the method yielding the greatest accuracy for the partial effects estimates, $h^(t,tz,z(tz))=h^(t−tz)z(tz)$⁠. The additional ridge penalization of the weight function enables us to identify the relevant window $τ(t)$ at which past exposures cumulatively affect the hazard at time t, as seen also in our application on football injury data. From a practical point of view, the presented modelling framework provides a suitable way to flexibly model training load exposures and analyse their effect on subsequent football injuries. Its flexibility in modelling the exposure window is a significant advantage, especially given the wide array of player load metrics available through GPS devices in professional sports. Being able to determine and analyse the time window in which load influences injury outcomes, without the need to fit a range of separate models (e.g. for 3,…, 6, or 7 days back), is highly appealing to practitioners in professional sports. We propose to use wide time windows to properly determine the exposure time at which, from that time on, the estimated effects are close to zero.

Moreover, the proposed model demonstrates a better overall performance compared to other alternative measures of training load exposures used in the literature. The widely known and used acute chronic workload ratio (ACWR) (Hulin et al., 2014), and its variants (Lolli et al., 2019; Wang et al., 2020), limit to summarize past observations into a predefined unweighted metric, through a ratio of two rolling averages –last 7 days (acute load) over last 28 days (chronic load). The same applies to the exponentially weighted moving averages (EWMA) (Williams et al., 2017) metric suggested as an alternative measure of rolling averages. While EWMA more accurately accounts for the decaying nature of fitness and fatigue effects over time than rolling averages, both may fail to accurately reflect several changes in past training exposures, as well as to consider predefined time windows that are either superfluous or insufficient. Instead, our method’s estimates of the cumulative effect and the relevant time window are based on the data.

Future research should evaluate negative binomial and zero-inflated models within the PAMM framework to address issues related to overdispersion and excess zeros (e.g. the low number of injuries). Although this is beyond the scope of the current study, comparing out-of-sample injury predictions from our model with other methods (e.g. regression models with ACWR-transformed training load variables and machine learning approaches like gradient boosting machines, neural networks, and random forests) would be a valuable area for further investigation. As this approach requires the estimation of many parameters and effects, which can be challenging when based on a limited number of events, another direction might involve pooling data from multiple teams or sources. This could be achieved using federated learning techniques (Archetti & Matteucci, 2023) to ensure data protection.

Simulation studies suggest that the model can recover a number of clinically plausible shapes for the true weight function under various levels of heterogeneity. Without prior knowledge about the form of association for time-varying exposures, the model proved to capture well a variety of shapes, estimating them from the data via P-splines. However, for nonsmooth effect shapes, like piece-wise constant or bi-linear, there exist other alternative methods such as adaptive splines (J. H. Friedman, 1991) or treed distributed lag nonlinear models (TDLNM) (Mork & Wilson, 2022), might be of interest.

Additionally, future research should explore the impact of the number of events per subject, held fixed in our simulation study, on model performance. It would also be worthwhile to investigate distributions other than Gaussian, such as Gamma-distributed random effects, which are commonly used in survival analysis (Balan & Putter, 2020).

By highlighting the potential value of PAMMs with WCE effects in assessing recurrent events in sports medicine, our study contributes to enriching the existing literature. We believe that this methodology would help in designing and comparing personalized training plans with insights regarding the risk of injury.

Acknowledgements

The authors thank the Medical Services of Athletic Club for data support and all the individuals who participated in the study.

Funding

This work was supported by the Basque Government through the BERC 2022–2025 program; of the Ministry of Science, Innovation and Universities through BCAM Severo Ochoa accreditation and through SEV-2017-0718 PRE2018-084007 funding; of AEI/FEDER, UE through the PID2020-115882RB-I00 with acronym ‘S3M1P4R’ and PID2023-153222OB-I00 with acronym ‘SPHERES’; and of Provincial Council of Bizkaia through the 6/12/TT/2022/00006 and acronym ‘MATH4SPORTS’.

Data availability

The application data used in this study is not publicly available due to privacy concerns. However, the simulation data and code used to generate the results are fully available and reproducible at: https://github.com/lzumeta/flex-mod-training-loads-recu-injuries.

Supplementary material

Supplementary material is available online at Journal of the Royal Statistical Society: Series C.

References

Author notes