Van Eetvelde et al.'s review (Van Eetvelde et al., 2021) presents mostly prospective case-control studies (McCullagh and Whitfort, 2013; Rossi et al., 2017; Thornton et al., 2017; Carey et al., 2018; López-Valenciano et al., 2018; Ruddy et al., 2018; Ayala et al., 2019; Oliver et al., 2020; Rommers et al., 2020), analyzing data retrospectively and predicting injuries based on the athletes' features value on a given day. In this prospective case-control study, we predicted the athletes' injury status in the upcoming 1 to 7 days. This will help the performance support staff anticipate athletes' upcoming injuries over the next training micro-cycle, which corresponds to 7 days within the Canadian National Female Short-Track Speed Skating Program.

We put a data mining system (Adriaans and Zantinge, 1997) in place within the Canadian National Female Short-Track Speed Skating Program to train and evaluate the random forest injury prediction models. Figure 1 provides a schematic of the different steps of data acquisition, preprocessing, feature engineering, and model training and evaluation. Each step presented in Figure 1 will be explained in the following subsections.

We collected data over the 2018–2019 and 2019–2020 seasons. Eleven women, members of the Canadian Short-Track Speed Skating National Team, consented to take part in the study. Every athlete signed an “athlete consent form” allowing the use of their data for research purpose and program evaluation. The consent form was approved according to provincial legislations and ethical recommendations. Sensitive information that can directly identify the athletes (e.g. names, addresses) were safeguarded and maintained under controlled conditions according to provincial laws. The athletes went through a familiarization and education process to ensure data were reported in a standard and appropriate way. Table 1 describes the characteristics of the athletes who participated in the study.

The conceptual model of injuries and sports performance, presented in Figure 2, was created to guide our interpretation of injuries and their underlying mechanisms and provide adequate perspective on the influence of the three categories of factors on injuries. It is inspired by Selye's general adaptation syndrome (Selye, 1950), Meeuwisse et al. (2007) dynamic model of sports injury etiology, and D'Isanto (2019) factors influencing sports performance. The model (Figure 2) guided the choice of variables to monitor. The amount of stress imposed on the athlete was quantified by monitoring the external (training stress imposed on the athlete) and internal training load (training stress perceived by the athlete) (Foster et al., 2017). Factors influencing stress response were monitored through psychological wellness questionnaire scores (Junge, 2000; Shrier and Hallé, 2011), tracking of athletes' injury type and location history, and monitoring of their history of training completion (ability to complete the prescribed training). The stress response was quantified through the tracking of the athletes' heart rate variability (HRV) (Goessl et al., 2017) and neuromuscular function (Gathercole et al., 2015).

We tracked 40 different variables through the monitoring process, grouped into eight subsets, which are described in Table 2. Combining the data from the 11 athletes over the two seasons provided a total cumulative dataset of 7,370 daily observations (a daily observation regroups the 40 monitored variable values with the associated injury outcome). The monitored data were stored in a database.

The monitoring process implies inevitable missing values (Benson et al., 2021). We observed three causes for missing values in the dataset: (1) the athletes did not train on a given day; (2) some variables were not measured daily; and (3) the athletes omitted to report data. Even if all athletes who participated in the study had participation rates ranging from 77 to 100% and averaging 90%, some daily instances were incomplete and displayed, for some variables, missing values. Imputation strategies had to be employed and are reported, for each monitored variable subset, in Table 2. Replacing missing values by the most recent measurement was the preferred imputation strategy for the psychological wellbeing metrics, neuromuscular function, and HRV variables, which were measured every 2 to 3 days. Although not perfect, the imputation strategy is easy to implement, avoids information leakage (Van Eetvelde et al., 2021), and provides indicators of the stress response variables over a short time frame. It is the most simplistic single imputation technique described by Benson et al. (2021).

Injuries were forecasted through a feature-based approach. New variables were computed from the athlete's monitored variables history, stored in the database, to reflect injuries' time dependence (Clarsen et al., 2013; Bittencourt et al., 2016) and improve the models' forecasting performances. We calculated the exponentially weighted moving average (EWMA) over the last 7 and 28 days (referred to as the acute and chronic indices) for each monitored variable. These time frame periods were chosen based on previous studies exploring the association between these indices and injuries (Gabbett and Jenkins, 2011; Hulin et al., 2014; Gabbett et al., 2016). It also corresponds to the microcycles (7 days) and mesocycles (28 days) used for the athletes' global periodization. The acute to chronic work ratio (ACR) (Gabbett and Jenkins, 2011; Hulin et al., 2014; Gabbett et al., 2016) is an additional feature engineered for all the variables. The EWMA ratios are a standard forecasting practice to improve the forecasting models' performances (Holt, 2004). In addition, we computed the acute to chronic difference (ACD), inspired by Banister's training impulse concept (Banister et al., 1975; Jovanovic, 2018).

For each variable, to quantify the influence of the past mesocycles and microcycles, we computed rolling means, standard deviations, and maximums over time windows of 7, 14, 21, and 28 days. At the end of the feature engineering process, a single daily observation is constituted of 1014 features associated with the athletes' injury status (injured or not injured).

To make predictions, the random forest algorithms were “trained” on available data and “learned” a set of rules that allow appropriate classification (injured or not injured), given the set of features provided as input. This set of features was composed of: (1) monitored variables on any given day; (2) imputed variables (if some monitored variables were missing); (3) engineered features based on the athletes' monitored variables history stored in the database.

To evaluate the classifiers' ability to replicate their predicting performances, they were evaluated on a test set, i.e., a set of unseen data (not used in the models' training process). Two sets of forecasting models were developed. One including only training load variables (TL), and a second (ALL) combining a wide array of monitored variables (neuromuscular function, heart rate variability, training load, psychological wellbeing, past injury type, and location). From the 7,370 available daily instances, we created two test sets of 1,100 instances, composed of each of the 11 athletes' last (chronologically) 100 daily instances. The ALL-Test-Set was composed of each of the 40 monitored variables and their respective engineered features, based on each athlete's complete data history, for a total of 1,014 features, while the TL-Test-Set was composed only of internal and external training load monitored variables and their corresponding engineered features, for a total of 182 features. With the remaining 6,270 daily instances, we created two training sets: the ALL-Training-Set composed of the same features as the ALL-Test-Set and the TL-Training set composed of the same features as the TL-Test-Set. As depicted in Figure 1, the test sets feature engineering was performed on the athletes' complete data history. The training sets feature engineering was only performed on the remaining 6,270 daily observations. This process replicates the way data will be handled to make predictions once the models are implemented in the training environment and avoid data leakage, where information of the test set would leak in the training set, leading to overfitting and biased interpretation of the models' evaluation (Van Eetvelde et al., 2021). Two distinct data sets (ALL and TL) were created to assess the impact of building forecasting models composed of variables from each of the three components of the injury conceptual model (ALL) in comparison to models composed of variables only reflecting the stress applied to the athlete (TL).

Two sets of models were trained using the Waikaito Environment for Knowledge Analysis (WEKA) (Hall et al., 2009). WEKA provides built-in functions to facilitate classifications of imbalanced datasets, such as the one in this study, where non-injured outweighs injured daily observations. To improve the random forest injury classifier, the injured and non-injured instances of the test set were balanced using WEKA's class balancer filter, which applies a weight function to the instances such that injured and non-injured instances display a similar total sum of weight (Hall et al., 2009). A cost function was also introduced using WEKA's cost-sensitive classifier function. The function penalizes injury misclassifications through the random forests training process, which improves injury forecasting performances (Hall et al., 2009). Default setups were used for all of WEKA's built-in functions. The ALL forecasting models were trained on the ALL Training Set, and the TL forecasting models were trained on the TL Training Set. For each set, seven classifiers were trained to predict whether the athletes would be injured in an upcoming time frame of 1 to 7 days, respectively. The random forest algorithm introduces randomness, notably through the bootstrap aggregating process (Breiman, 2001). To account for the classifiers' random factors, the training of each classifier was performed 30 times for each of the seven forecasting time frames. At the end of the training process, we obtained, for both sets of models (ALL and TL), 30 trained random forest classifiers for each forecasting time window (1 to 7 days).

Each instance of the test-set was provided as input to the forecasting models, and the trained classifier generated predictions on the upcoming injury status (on a time frame of 1 to 7 days depending on the classifier). To obtain a confusion matrix, each prediction was compared to the true outcome, which is already known. Confusion matrices compile True Positives (TP), False Positives (FP), False Negatives (FN), and True Negatives (TN) (Visa et al., 2011). We computed the following three metrics from confusion matrices that summarized each of the 30 models' predictive performances for every forecasting time window (1 to 7 days).

The sensitivity.

The specificity

the Matthews Correlation Coefficient (MCC) (Matthews, 1975).

The sensitivity and specificity are commonly used evaluation metrics and were computed to compare our models with existing sports injury prediction models (Van Eetvelde et al., 2021). The MCC (Matthews, 1975) is considered one of the best methods to evaluate binary classifiers (Powers, 2020). It is a wellbalanced metric that can be used even if the two classes are unbalanced (Boughorbel et al., 2017). An MCC of 1 indicates complete agreement between the predictions and the observation, and inversely, an MCC of −1 denotes complete disagreement (Boughorbel et al., 2017). An MCC of 0 would be considered a model where the predictions were made randomly (Boughorbel et al., 2017). The unbalanced nature of the data set is the reason why accuracy (Yin et al., 2019) was not used as an evaluation metric. In our case, out of the total 7,370 daily observations, only 12% constitute injuries. Therefore, even if the models were always to predict non-injured athletes, the accuracy would remain relatively high (~0.88) and would not provide valuable information to the practitioners. We computed the mean and standard deviation of the 30 sensitivity scores, 30 specificity scores, and 30 MCC of every forecasting time window. At the end of the evaluation process, we obtained three evaluation metrics (mean ± standard deviation) for each model and each forecast window (1 to 7 days) for both the ALL and TL forecasting models.

To analyse the effect of the forecasting model types (ALL or TL) and forecasting time windows (1 to 7 days) on each evaluation metric (sensitivity, specificity, MCC), two-way ANOVAs were performed on each metric. In addition to ANOVAs, paired T-tests (Meek et al., 1987) were performed on each separate forecasting window to evaluate the effect of the model type on the evaluation metrics, with the null hypothesis that both forecasting models (ALL and TL) evaluation metrics were equal at a significance level of p < 0.01. In addition, we computed Cohen effect sizes (Cohen, 1994; Fritz et al., 2012) to quantify the difference between the two types of model evaluation metrics.

Random forest classifiers enable highlighting variables of importance, which have a practical application for the performance support staff, trying to understand better the causes and mechanisms of injuries (Meeuwisse et al., 2007; Bittencourt et al., 2016; Van Eetvelde et al., 2021). The higher the random forests' minimum GINI decrease coefficient, the greater the importance of the variable in the model (Archer and Kimes, 2008).

The conceptual model of injuries presented in Figure 2 is essential for developing an injury prediction machine learning framework. Practitioners who want to elaborate a data mining system need to have in mind the three conceptual components of injuries, i.e., the stress applied to the athlete, the factors influencing the athletes' stress response, and the stress response itself, which can lead to positive adaptations, or negative outcomes like injuries. They should define key monitoring variables reflecting each component of the conceptual model of injuries, which could vary depending on the context (exercise modality, athletes' level, age, financial resources, etc.) and should be adapted to match the practitioners' needs (Halson, 2014; Starling et al., 2020; Clarsen et al., 2021). Monitored data should be stored in a database, as engineered features based on the data history are essential to reflect the time-dependent nature of injuries (Bittencourt et al., 2016). Handling the database is easier when a problem conceptualization has been made (Figure 2), and a data-mining plan (Figure 1) has been established. This planning provides perspective on the problem and ensures data are monitored and stored purposefully.

Missing data have to be handled (Benson et al., 2021). In this study, we used the simplest imputation form, which replaces the missing data with the most recent measurement. The imputation strategy does not minimize the errors and introduces information loss that could potentially be detrimental to the models' predictive performances. Other imputation strategies that use mean values or regression techniques minimize the models' error prediction (Benson et al., 2021). However, these strategies are less optimal in the perspective of the models' implementation and could introduce higher risks of information leakage and overfitting (Van Eetvelde et al., 2021). To minimize the missing values' imputation detrimental effects on future models' performance, we will be looking at minimizing missing values, by measuring certain variables more frequently. It is a challenge to maximize measurement frequency without overwhelming the athletes (Starling and Lambert, 2018). New monitoring technology to track the athletes' HRV and neuromuscular function may represent an interesting avenue for that purpose (de Zambotti et al., 2019; Roberts et al., 2020). Practitioners are free to use the imputation strategies that better suit their needs. However, they should always keep in mind the models practical value and overfitting risks in their decision-making.

The ALL forecasting models show promising results, especially over short time windows of 1 and 2 days, where more than half the future injury days and more than 70% of future injury-free days can be predicted correctly by the models. Injuries can be compared to extreme meteorological events, which are particularly difficult to forecast (Ghil et al., 2011; Nayak and Ghosh, 2013; Boers et al., 2014). Although modest, the results suggest that, in short time frames, the models could bring added information and help the performance support staff, providing better reaction time and facilitating their injury management decisions. In previous studies, reported sensitivity scores usually ranged between 0.65 and 0.85, while reported specificity scores ranged between 75 and 95% (McCullagh and Whitfort, 2013; Whiteside et al., 2016; Rossi et al., 2017; López-Valenciano et al., 2018; Ayala et al., 2019; Rommers et al., 2020); every study differs in context and methodology.

In this study, methodological choices were made regarding to the models' implementation in the speed skating environment, which could explain the lower sensitivities and specificities. We brought a novel approach that predicts the athletes' health status in upcoming 1 to 7-day time windows. Inevitably, forecasting implies more uncertainty and diminished performance. Injury forecasting was preferred as it has more potential to help and assist the performance support staff. In the experimental setup, we replicated the models' implementation process in the training environment, where models will be trained with all the available information and data history at hand and where new data entries, to predict future outcomes will be provided by the same athletes who participated in the models' training process. Therefore, the test sets data was collected on the same athletes that contributed to train the models. We did not test the models' ability to predict injuries from athletes that didn't contribute to train the models. The training sets' engineered features were computed separately and are not influenced by the test set data. These choices introduce some form of overfitting and information leakage (Van Eetvelde et al., 2021), because new data entries for injury predictions will inevitably be influenced by the athletes' data history used to train the models. However, this process must be seen as a calibration process rather than a major hurdle or methodological issue. The models are calibrated to predict the upcoming injuries of the athletes within the Canadian National Female Short-Track Speed Skating Program and could only work in that particular context, with the specific athletes who took part in the data-mining process. Practitioners should adapt the framework to work in their specific environment and context, and with the athletes participating in their data-mining framework.

This study appreciates the significant impact of introducing variables from all three components of the injury conceptual model on the injury forecasting performances. Keeping the methodological aspect constant, the forecasting models, composed of variables from the three components of the injury conceptual model (ALL forecasting models), have displayed significantly higher MCC (p < 0.01), with large Cohen effect sizes (d >1) on each of the seven forecasting windows, compared to models built only using training load (TL forecasting models). In fact, the TL forecasting models all had an MCC below zero. Therefore, they behave similarly to a model making predictions randomly and would not be helpful for the performance support staff.

The ALL forecasting models all have an MCC above zero. The MCC values gradually decrease to near zero (0.00 ± 0.02) on 7-day forecasting windows. As the forecasting window gets larger, the ALL forecasting models' ability to predict injury occurrences from the test set decreases. This is expected, especially for complex problems such as injuries, which are influenced by many interrelated variables and some random factors. For instance, injuries can occur accidentally by falling on the ice. These types of injuries are unfortunately impossible to predict. As the forecasting windows get larger, the influence of the different variables and random factors are multiplied, thus explaining the decrease observed in ALL forecasting models' sensitivity and MCC. As for the ALL forecasting models' specificity increase with larger forecasting windows, there are two main reasons. The first one is that patterns leading to an athlete staying healthy are probably better defined and predictable by the models. There are more cases where athletes are healthy, which facilitates pattern recognition. While injuries may be compared to extreme meteorological events, injury-free days resemble day-to-day weather, which is relatively predictable. The second reason is the cost function for injury misclassification, which ties the sensitivity and specificity together. As injuries get harder to predict, the only way to minimize the cost function with larger forecasting windows is to maximize the models' specificity. We applied the same cost function for misclassification on every forecasting window. It would be interesting to test different cost functions depending on the forecasting windows.

It is crucial for the performance support staff to understand the variables that have the strongest influence on athletes' injuries. These variables can provide valuable information on injury risk factors and guide their interventions (Meeuwisse et al., 2007; Clarsen et al., 2021; Van Eetvelde et al., 2021). The random forest algorithms were useful to provide the most important variables [maximizing the minimum GINI decrease (Archer and Kimes, 2008)] of the ALL forecasting models to forecast injuries over the seven time windows. The most important variables are derived from three main features: athletes' ability to complete their training (training completion); lower limb and trunk/lumbar injuries; and sFatigue, an internal training load marker. The training completion variable is intimately tied to our definition of injury (Meeuwisse et al., 2007; Clarsen et al., 2021). Injuries impede the athletes' ability to complete their training because of physical or psychological complaints. Tracking the athletes' ability to complete their training sessions over time is crucial to identify future injuries. Lower limb injuries are another variable worth tracking in short-track speed skating. Lower limb asymmetries have already been pointed out as potential injury risk factors in this sport (Konieczny et al., 2020). In addition, the athletes skate in a crouched posture, which could be responsible for lower back and lumbar discomfort (Hesford et al., 2012). Both variables (training completion, as well as lower limbs and trunk/lumbar injuries) provide indices of the athletes' health status, a factor influencing the athletes' future stress response. The other monitored variables may not be as important for the ALL forecasting models' predictions, but they certainly play a role in the forecasting process. As a matter of fact, sFatigue, an internal training load marker, was among the most important variables for five-day forecasting windows. In addition, as the forecasting windows got larger, the minimum GINI decrease coefficients of the important variables got smaller, potentially linked to the ALL forecasting models' MCC and sensitivity decrease. Practitioners building models in different sports and contexts will necessarily observe different variables of importance.

The injury conceptual model presented in Figure 2, although useful to understand the interactions of different categories of variables and their influence, remains simplistic. The etiology of injuries remains much more complex (Meeuwisse et al., 2007; Bittencourt et al., 2016). We will keep improving the conceptual model of injuries, making it more robust and accurate, based on the continuous observations provided by the performance support staff. Improvements to the conceptual model will help orient the data mining framework and future variable monitoring process. Thus far, we have focused on the negative stress aspect of the conceptual model presented in Figure 2 (injuries). However, in the future, basing our work on the same conceptual model, we could explore the positive stress adaptations to predict performance through a similar framework.

The random forest forecasting models act as black boxes. Even if they provide some indicators of the variables of importance, they cannot assist the performance support staff in understanding how each variable influences the injury conceptual model. Nonetheless, the ALL forecasting windows on short time windows (1 or 2 days) are promising, and we believe they could add value to the performance support staff. We plan to implement the models similarly to meteorological models, assisting the performance support staff and providing injury warnings in the upcoming 1 to 7 days. Doing so could facilitate the performance support staff interventions. Much like meteorological models are used in real life, they do not override human judgment. Ultimately, decisions are based on all available information, some of which is not necessarily taken into consideration by the forecasting models. The implementation of the models will be closely documented to help validate and improve our framework.