Probabilistic graphical models, which integrate graphical structures into probabilistic inference, are widely used and effective frameworks for studying these complex systems. The concept is to factorize the joint probability distribution p for the variables X=(X1,⋯,Xm)⊤ concerning a graph G=(V,E), where V is the set of vertices representing the variables and E the set of edges encoding the independence relationships (31, 32). Bayesian networks are a special class of probabilistic graphical models, where G is a directed acyclic graph (DAG) or named DAG models (33, 34). The joint probability distribution p can be specified by a set of parameters θ and factorizes based on G as:

Where θ=∪i=1mθi, x is a realization of X, and we assume that the subsets {θi}i=1m are disjoint. Denote the parents of node i by pa(i)), where there is a directed edge from j to i if j∈pa(i). Thus, Equation (1) can also be interpreted as stating that a variable x_i is conditionally independent of its non-descendants, given its parents X_pa(i) in G. This is the Markov property (31).B=(G,θ) denotes a Bayesian network. Given a data sample X, learning a Bayesian network, therefore, involves estimating both the network structure G and θ. If the joint distribution of X is a Gaussian distribution, then

In the case of Gaussian data, to address the uncertainty regarding the graphical structure, Maathuis et al. (35) provided lower bounds for causal effects after identifying a Markov equivalence class that is consistent with the data. By using a Bayesian approach, one can combine structure learning and effect estimation into a process that produces the posterior distribution of causal effects. A significant advantage is that this method accounts for both graphical and parameter uncertainty, as first proposed and demonstrated in a psychology application by Moffa et al. (36) for binary data.

where the matrix Ω = Σ⁻¹ is symmetric, positive definite, and Markov relative to G. Under the assumption of normality, the Gaussian DAG model is almost always faithful to the DAG within the parameter space, which means that the conditional independence relationships implied by the distribution are precisely the same as those represented by the DAG through the Markov property (37). For a Gaussian DAG model, we can rewrite the factorization in Equation 1 as Equation (3): (ϕ denotes the normal univariate density function) (33).

The do-operator is a fundamental concept in causal inference, proposed by Pearl (37). Provides a theoretical framework for quantifying the effects of interventions in a system. By modifying the observed probability distribution, the do-operator helps distinguish between correlation and causation. It is a core tool in causal inference and decision theory, allowing researchers to systematically answer the question “What would happen if?”

DAG provides an alternative approach to causal inference. Using DAG to describe the data generation process is appealing because the edges of the graph naturally represent the causal relationships between variables. Under a known DAG, the do-operator determines the causal effect of one variable on another. We have discussed the DAG model obtained by the Gaussian model, but we are interested in the causal relationships between the nodes, as shown in Figure 1.

Our goal is to determine the causal effect of variable X_i on the variable X_o. We use Pearl's do-operator to describe the effects of intervention (37), where the distribution of X_o under an intervention in X_i is generally indicated as ℙ(X_o = k∣do(X_i = l)). Changes in the distribution or shifts in the distribution of the outcome variable across different levels of the intervention variable often serve as target estimands with practical significance (38). Evaluating and contrasting the change in the probability of X_o belonging to level k, when the intervention variable X_i is set to level l′ versus level l offers a measure of the distribution shift:

for each l≠l′ and l,l′∈{1,…,Li} and k∈{1, …, L_o}. We can evaluate ordinal causal effects (OCE) as represented by the target causal estimands in Equation (4).

In the data X=(X1,⋯,Xm)⊤, the random variable X_k = (x₁, x₂, ⋯ , x_n), where k = 1, 2, ⋯ , m, often includes both continuous data and ordinal data. Before causal analysis, data must be organized and optimized to ensure precision and reliability of the results. Inconsistent data types can have a significant impact on analysis results. For example, in regression analysis, it is necessary to standardize the data to ensure consistent units of measurement. The objective of this paper is to analyze ordinal data, so the first step is to properly transform unordered data into ordinal data, which supports the subsequent analysis. Therefore, we propose a quantile threshold function (QTF) that transforms continuous variables into ordinal variables and ensures the consistency of the classification results while effectively preserving the ordered nature of the data.

In this paper, we apply the kernel density estimation to fit the probability density of the continuous variable X_k in Equation 5. A Gaussian kernel function is chosen as the smoothing kernel in Equation 6, which has the advantage of not requiring a predefined data distribution shape. Through bandwidth adjustment, it effectively approximates unknown distributions. This method overcomes the dependence on distributional assumptions inherent in traditional parametric methods.

Where f^(x) represents the probability density function, K(u) is the Gaussian kernel function, and h is the bandwidth, which is a key smoothing factor in kernel density estimation. We adopt the bandwidth selection method proposed by Ripley (39) in his book Modern Applied Statistics with S as the criterion (page 127) (39). The definition of the quantile threshold function is as follows:

Definition 1 (quantile threshold function). The random variable X = (x₁, x₂, ⋯ , x_n) has a probability density function f(x) and a distribution function F(x). If there exists a non-negative real-valued function g(x) such that:

where Q_i satisfies F(Qi)=∫-∞Qif(x)dx=i4,i∈{1,2}, then, the function g(x) is called quantile threshold function

g(x) is also a random variable. The value of x ≤ Q_i in X is defined as the lower level (assigned a value of 0), and within the Q_i<x ≤ Q_i+1 interval is defined as a medium level (assigned a value of 1), and we classify values of x>Q_i+1 as a high level (assigned a value of 2). Based on the above definition, we can naturally deduce the following conclusion.

Proposition 1. The sum of the probabilities of the three categories equals 1, that is, ∑j=02ℙ[g(x)=j]=1

Proof. See the Supplementary material: Proof of Proposition 1.

Definition 1 extends to multiple classification scenarios, as shown in the following Definition 2.

Definition 2. The random variable X = (x₁, x₂, ⋯ , x_n) has a probability density function f(x) and a distribution function F(x). If there exists a non-negative real-valued function g(x) such that:

where Q_i satisfies F(Qi)=∫-∞Qif(x)dx=in+2,i∈{1,2,⋯,n}.

Proposition 2. The sum of the probabilities of the categories n+1 is equal to 1, that is, ∑j=0nℙ[g(x)=j]=1.

Proof. See the Supplementary material: Proof of Proposition 2.

Proposition 3. limn→∞ℙ[g(x)=n]=0.

Proof. See the Supplementary material: Proof of Proposition 3.

From Proposition 3, we know that when performing ordered classification on random variables, the classification must be finite. Therefore, the function g(x) achieves a smooth transition from continuous variables to ordinal variables. g(x) provides a structured data representation that maintains both information retention and interpretability for subsequent analysis.

In this paper, we use n = i = 2 as the classification criterion. After the classification is completed, we conduct hypothesis testing on the results. This study employs the Analysis of Variance (ANOVA) method to verify the significance of differences between groups for reconstructed ordinal variables (40). Specifically, our objective is to determine whether the differences among the three groups are significant, which can be achieved by testing whether the means of each group are the same. We set the significance level at α = 0.05. Let

Where μ₀, μ₁, and μ₂ are the means of the three groups. If p < 0.05, then H₁ holds, indicating that there are significant differences among the three groups, which suggests that the above classification is effective. Conversely, if H₀ holds, it indicates that the classification levels are not significant.

We obtained ordinal data X*=(X1,⋯,Xm)⊤ by the function g(x). What we are interested in is the DAG that represents the relationships between ordinal variables. We introduce the Gaussian DAG-models in Subsection 2.1, but this model is built under the assumption that the data follow a Gaussian distribution. Therefore, for the construction of the DAG model for ordered data, we have used the OSEM algorithm proposed by Luo et al. (29). By assuming that each ordinal variable is obtained by marginally discretising a set of Gaussian variables, we can get the ordinality amongst the categories. And Gaussian variables jointly follow a DAG structure G=(V,E). The OSEM algorithm provides a new framework that effectively learns Bayesian networks from ordinal data and captures the orderliness among categories.

Let X^* be a set of m ordinal variables, where X_k takes values in the collection of {τ(k, 1), τ(k, 2), …, τ(k, L_k)} and τ(k, 1) <τ(k, 2) <⋯ <τ(k, L_k), k = 1, …, m. We assume that the number of levels L_k≥2, therefore, each variable should at least be binary. It is typical to set τ(k, l) = l−1 for all 1 ≤ l ≤ L_k, i.e. τ(k, 1) = 0, τ(k, 2) = 1, and so on. Further, we assume that each X_k is obtained by discretising an underlying Gaussian variable Y_k using the thresholds −∞ =:α(k, 0) <α(k, 1) <⋯ <α(k, L_k−1) <α(k, L_k): = ∞. Let αi=(α(k,0),…,α(k,Lk))⊤ and α={αk}k=1m. Thus X_k is defined by the following rule:

Of course, Y=(Y1,⋯,Ym)⊤ are unobservable, and we observe X_k obtained from the continuous variables by discretisation. The diagram in Figure 2 offers a visual depiction of the setup in an example case with a few variables.

Formally, Luo et al. (29) proposed a DAG model for ordinal variables based on latent Gaussian variables. The model makes the following optimization definition based on Equation 2. The detailed derivation process can be found in the paper (29).

The OSEM combines the multinomial probit model (41) and the structural EM algorithm of (42) to solve the problem of learning Bayesian networks from ordinal data. Specifically, the method proposes an iterative scoring and search strategy - the Ordinal Structural EM (OSEM) algorithm for learning Bayesian networks from ordinal data.

Consider the general latent Gaussian DAG-model of Section 3.2. We are interested in computing the target causal estimand in Equation 4. For example, in Figure 2, the intervention variable X_i is X₁, and the outcome variable X_o is X₀. The Equation 4 can be written as Equation 11, representing the OCE on X_o of an intervention on X_i. When the intervention variable X_i is set to level l′ versus level l offers a measure of the distribution shift (37):

The direct computation of Equation 11 would result in 0 for each level of the intervention and outcome variables because there is no causal path between the ordinal X_i and X_o in the DAG. But it is evident that they are causally related to each other by Y. Therefore, we can consider that if we intervene on the latent variable Y_i in a way that changes the level of its ordinal child variable X_i, and then compute the effect of this intervention on the latent parent Y_o of X_o, it is possible that the level of X_o could also change as a result. The potential change in the level of X_o resulting from an intervention on the latent parent of X_i is the OCE studied in article (30). Using the α={αk}k=1m, the target causal estimand in Equation 11 on the ordinal variables can be equivalently computed as the following [(30); Definition 1,page9]:

for each 1≤k≤Lo,1≤l,l′≤Li, with l≠l′. The definition of OCE is anti-symmetric for the initial and end level of the intervention variable, implying that

Based on the above Equation 12, the Latent Gaussian DAG-model establishes a relationship for calculating the intervention effect between ordinal variables. This allows for the computation of the OCE between variables X^*, which is equivalent to the intervention effect between variables Y. The (30) provides a detailed proof and derivation of Equation 12 in both the main text and the Appendix. Building on this result, Scauda et al. (30) proposed Proposition 5 - a method for computing OCE. The specific details can be found in Proposition 5 (Computation of the Ordinal Causal Effect) on page 12 of the (30). According to Proposition 5, we can calculate the OCE efficiently. Figure 3 illustrates the flowchart of the proposed algorithm.

Maintaining an injury-free condition is a crucial factor for success in sports. Although injuries are difficult to predict, the application of emerging technologies and data science can offer valuable insights. Even with a well-specified model, inaccurate data can compromise causal analysis. Data quality is often more critical than sample size (43), particularly in sports science, where physiological measures (e.g., maximal oxygen uptake, gene transcription activity) are inherently noisy due to biological and technical variability. Additionally, exercise intervention studies face challenges such as participant dropout, missing data, measurement errors, and inconsistencies in data processing, all of which hinder reliable causal interpretation (44).

This study utilizes a comprehensive training log dataset from Kaggle (45, 46) collected by a Dutch team in 2012–2019. This dataset, developed by Lövdal et al., employs machine learning to predict injuries based on data that focuses on middle- and long-distance events (800m to the marathon), with detailed performance records for 74 athletes (27 females and 47 males). The dataset adhered to the ethical principles of the Declaration of Helsinki and received formal approval from the ethics committee (45). This study follows the data structure established in (45), utilizing the Day approach [(45), page 1523, Table 1], with this research focusing exclusively on the data from a single day to analyze the causal effects. This dataset contains a total of 42,766 samples. Table 1 presents detailed characteristics of the variables. Variables 1, 6, and 11 are directly considered ordinal, while the remaining variables require ordinal classification using QTF.

We simulated the probability density function for continuous variables and classified the variables based on QTF. The graph of the probability density function and its corresponding quantile ranges are presented in the Supplementary Figure 1. Taking variable 2 as an example, from Supplementary Figure 1 we can observe that Q₂ and Q₃ divide the range of the values of the variable into three parts. The portion less than Q₂ is defined as low level and assigned a value of 0. The values between Q₂ and Q₃ are categorized as medium level and assigned a value of 1. Finally, values greater than Q₃ are defined as high level and assigned a value of 2. From Supplementary Figure 1, we can see that the probability of values at the low level (green) is 0.5, the probability at the medium level (orange) is 0.25, and the probability at the high level (purple) is also 0.25. Through this method, we transform the variable from a continuous variable to an ordinal variable. Subsequently, we perform ANOVA on the three groups levels and conduct a significance test (Table 2) (40).

The mean values of the variable increase monotonically from group “0” to group “2”, and the standard deviations also progressively increase. For variable 2, the means and standard deviations at different levels are as follows: In group “0”, the mean is 3.186, the standard deviation is 3.991, and the sample size is 30, 437; in group “1”, the mean is 12.882, the standard deviation is 0.98, and the sample size is 5, 996; in group “2”, the mean is 20.016, the standard deviation is 5.2547, and the sample size is 6, 333. Other variables exhibit similar patterns, where the mean in group “0” is relatively small, while the means in groups “1” and “2” are significantly higher (e.g., for variable 5, the mean increases from 0.017 to 2.536). Several variables have extremely high F-values (e.g., 156, 092 for the variable), showing that the between-group variance is far greater than the within-group variance. The ANOVA results indicate that the differences in means across groups are statistically significant (p < 0.05), suggesting that the distribution of the variable differs significantly across levels and shows a clear gradient pattern (group “0” <“1” <“2”). This result implies that the method used can effectively distinguish differences between levels, leading to the rejection of the null hypothesis H₀ and the acceptance of the alternative hypothesis H₁.

In the Supplementary material, we present the ANOVA results for the three discretization methods: equal-width, equal-frequency, and k-means. Supplementary Table 1 shows the ANOVA results for equal-width. This method uses intervals of the same width, making it simple and intuitive. However, the ANOVA results indicate that its F-value is lower, suggesting the between-group differences are less significant than with the QTF method. Moreover, for outlier variables (such as Variable 4 or Variable 5), equal-width fails to reflect the data distribution accurately. Supplementary Table 2 shows the ANOVA results for equal-frequency. This method ensures that each interval contains roughly the same number of samples. When sample values are highly concentrated and repeated, the resulting zero mean and standard deviation for some groups reduces the interpretability of the data. Supplementary Table 3 presents the results of k-means clustering. The ANOVA shows large F-values for the three clusters, suggesting significant between-group differences. However, k-means is a typical data-driven method with unfixed boundaries. If clusters are ordered by their mean values to define “low,” “medium,” and “high” levels, the boundaries will change with the data because cluster centers are randomly initialized. This randomness weakens theoretical interpretability.

We followed the approach described in Luo et al. to derive DAG estimates from 500 bootstrap samples of the data (29, 30), utilizing the OSEM algorithm with a Monte Carlo sample size of K = 5 and a penalty coefficient of λ = 6. And the resulting CPDAG is shown in Figure 4.

Figure 4 illustrates the causal relationships between different variables, highlighting both direct and indirect variables surrounding the Injury variable. The direct influencing variable for Injury is Kmsprinting, which has a significant impact on the occurrence of Injury. Kmsprinting training typically places high stress on muscles, joints, and ligaments, especially when the intensity is excessive or recovery is insufficient, making it prone to cause tissue damage or excessive fatigue. Kmsprinting may lead to rapid muscle contraction and high load in a short period, thereby increasing the risk of exercise-related injuries. Perceivedrecovery and Perceivedtrainingsuccess are also direct variables influencing Injury. Perceivedrecovery plays a critical role in injury risk. When perceived recovery is poor, muscles and joints may not withstand higher loads, potentially leading to improper movement and decreased endurance, thus increasing the risk of injury. Perceivedtrainingsuccess directly affects injury risk. If an athlete perceives training success as high, it may indicate good physical condition and effective recovery, thereby reducing the risk of injury. Conversely, a lower perception of training success may indirectly reflect accumulated fatigue and insufficient recovery, increasing the probability of injury.

Injury's indirect influencing variables, Pathway 1 is Totalkm→Sessions→Perceivedrecovery→Injury: Increasing total running distance can lead to higher training intensity and frequency. A rise in total distance is often accompanied by increased training frequency, which may result in insufficient recovery time. Excessive training frequency and load can compromise recovery quality, leading to poor perceived recovery and indirectly increasing the risk of injury. If increases in running distance and frequency are not balanced with adequate recovery and proper load management, the perceived recovery level may decline, significantly elevating the risk of injury.

Pathway 2 is Strengthtraining→Perceivedtrainingsuccess→Injury: Strengthtraining enhances muscle strength, joint stability, and exercise efficiency, thereby improving athletic performance and training outcomes. It also increases athletes confidence and perception of training success, which may indirectly indicate improved physical adaptation and recovery levels. Perceivedtrainingsuccess can reduce the risk of injury, as the body is in better physical condition, movements are more precise, and energy distribution is more efficient. Conversely, a lower perception of training success may have the opposite effect. Based on the above analysis, it is recommended to control the intensity and frequency of Strengthtraining to prevent muscle injuries caused by overtraining. Attention should be placed on athletes Perceivedrecovery and sense of Perceivedtrainingsuccess, and by adjusting the training plan, recovery outcomes can be effectively improved.

To visually represent the bootstrapped estimates, we present the adjacency matrices of the DAG derived using OSEM, converted to CPDAG, as a heatmap in Supplementary Figure 2 of the Supplementary material. The intensity of each cell corresponds to the frequency with which each edge appears in the bootstrapped samples. The shade in the grid indicates the proportion of times a directed edge occurs in the 500 bootstrapped CPDAG, with an undirected edge being split equally between both directions. Darker shading corresponds to a higher frequency of the respective directed edge. Additionally, we examine the causal relationship along the most frequently observed directed edge in the 500 bootstrapped CPDAG by estimating the ordinal causal direct effects of Kmsprinting (variable 5) on Injury (variable 11) within the sample's DAG. Raincloud plots, which incorporate histograms and boxplots of the estimated effects, are displayed in Figure 5.

Figure 5 shows the ordinal causal direct effects of Kmsprinting (variable 5) on Injury (variable 11). “Injury=1” indicates no injury, while the level of “Injury=2” indicates injury. We use the do-operator to describe this result.

When the “Injury=1”, and the level of variable 5 changes from 1 → 2, the OCE <0. That means ℙ[injury = 1∣do(kmsprinting = 2)]−ℙ[injury = 1∣do(kmsprinting = 1)] <0. This result suggests that increasing the intensity of sprint training significantly increases the injury risk for uninjured athletes. High-intensity sprint training may lead to excessive load on muscles and joints, exceeding the body's ability to adapt, thereby increasing the likelihood of injury. Therefore, for uninjured athletes, maintaining a lower level of sprint training helps to minimize the injury risk. If the sprint training level changes from 1 → 3, or from 2 → 3, the OCE value remains negative. This finding indicates that increasing the intensity of sprint training has a significant impact on the injury probability for uninjured athletes, especially during the transitions from 1 → 2 and from 1 → 3. Among these three intervention levels, the absolute value of the OCE mean is largest for the intervention from 1 → 3, suggesting that transitioning directly from low to high intensity training is particularly dangerous for uninjured athletes. Therefore, training plans need to be carefully designed.

When the “Injury=2”, and the level of variable 5 changes from 1 → 2, the OCE > 0. That means ℙ[injury = 2∣do(kmsprinting = 2)]−ℙ[injury = 2∣do(kmsprinting = 1)]>0. When the sprint training level increases from 1 → 2 or 1 → 3, the probability of injury significantly increases. For injured athletes, continuing to increase the sprint training intensity before full recovery can worsen the existing injury or lead to incomplete recovery, which significantly increases the risk of injury. During this stage, athletes should avoid increasing sprint training intensity and prioritize basic training and recovery. When the sprint training level changes from 2 → 3, the OCE for both injured groups is more concentrated, with the mean close to 0, indicating that the impact of increasing from medium to high-intensity sprint training on injury is relatively small. This limited effect may be due to the athletes' adaptation to medium-intensity training and their more stable physical condition. Therefore, after sufficient medium-intensity training, gradually increasing high-intensity sprint training does not significantly increase the injury risk.

In summary, the transition from low to medium intensity (1 → 2) is a critical period for injury risk, requiring close monitoring of an athlete's physical condition and recovery. After reaching medium intensity, transitioning to high intensity is relatively safe, as the body has developed some level of adaptation, resulting in a lower injury risk. For uninjured athletes, it is recommended to prioritize a progressive increase in training load and avoid a direct jump from low to high intensity (1 → 3). For injured athletes, high-intensity sprint training should be strictly limited during the recovery period, and priority should be given to recovery training. High-intensity sprint training requires a focus on recovery quality after training to ensure athletes maintain a good recovery state during the gradual increase in intensity, minimizing fatigue accumulation and injury risk. The Kmsprinting has a significant causal effect on Injury. The risk from medium to high intensity (2 → 3) is relatively small, suggesting that training upgrades should be based on the athlete's adaptation. Scientifically planning sprint training intensity and pace, combined with recovery training, helps reduce injury probability while improving training effectiveness and safety.

Figure 6 shows the ordinal causal direct effects of Kmsprinting (variable 5) on Kmmiddle (variable 3) and Kmhigh (variable 4). The change of variable 5 from low level to medium level (1 → 2) shows that for low-level Kmmiddle athletes, an increase in sprinting results in a positive OCE value, indicating that improvements in sprinting may indirectly help these athletes enhance their middle-distance running ability. However, for medium and high-level Kmmiddle athletes, the OCE value is negative, suggesting that an increase in sprinting may decrease their performance in middle-distance running. A similar pattern applies to Kmhigh athletes. For low-level Kmhigh athletes, an increase in sprinting may improve their high-intensity running ability. In contrast, for medium and high-level athletes, an increase in sprinting may limit their high-intensity running ability.

It is worth noting that the OCE for Kmhigh is more scattered, indicating significant variation in the impact of sprinting on high-intensity running. The change of variable 5 from medium level to high level (2 → 3) shows that for both Kmmiddle and Kmhigh, the OCE are close to 0. This result suggests that once athletes have adapted to a certain level of sprinting load, further increases in sprint intensity have minimal effect on middle-distance running and high-intensity running performance. For athletes with weaker middle-distance and high-intensity running abilities, increasing sprint training can indirectly enhance their running abilities by improving neuromuscular adaptation and rapid power output. For athletes who already possess strong middle-distance and high-intensity running capabilities, increasing sprint training may lead to overloading or accumulation of fatigue, thereby affecting other running abilities. Once sprinting has reached a high level, athletes' bodies gradually adapt to the high-intensity load, and further increases in sprinting have limited intervention effects on other running abilities.

Therefore, we can conclude that for low-level athletes, appropriately increasing sprint training can help improve their middle-distance and high-intensity running abilities. For medium and high-level athletes, the ratio of sprinting to other running training must be balanced to avoid excessive sprint loads that may affect overall running performance. Regularly monitoring athletes' performance in different running abilities and adjusting the sprinting load based on data feedback will ensure maximized training effects while preventing fatigue accumulation. As sprint load increases from medium to high intensity (2 → 3), gradual adaptation should be emphasized to avoid sudden increases in training intensity that could negatively impact middle-distance or high-intensity running abilities. Through reasonable design and adjustment of sprint training, coaches and athletes can enhance specific abilities while minimizing the risk of injuries caused by overtraining, ultimately achieving comprehensive optimization of athletic performance.

Figure 7 shows the ordinal causal direct effects of other variables on Injury (variable 11). From Figure 7, we can observe that the OCE of Kmmiddle(variable 3) and Hoursalternative(variable 8) on injury are close to zero across various levels of change. This finding indicates that these two variables have a minimal impact on injury risk. Kmmiddle typically involves moderate-intensity training, and athletes' bodies are generally more adaptable, so it does not significantly increase or decrease injury risk. Hoursalternative mainly refers to low-intensity recovery exercises, which place minimal stress on the body, thus limiting their impact on injury risk. These two variables can be part of a stable training load to alleviate the physical stress from high-intensity training, thereby improving the overall safety of the training program.

In contrast, the changes in Strengthtraining(variable 6) show the most significant fluctuations in the OCE for injury, indicating a significant causal intervention effect. This finding suggests that strength training can enhance athletic performance but may also increase injury risk if the intensity is too high or the progression is too rapid. High-intensity strength training places a heavy load on muscles and joints, and without adequate recovery or adaptive training, it can lead to muscle strains, ligament injuries, and other issues. Additionally, Strengthtraining's effects vary significantly among individuals, as different athletes have differing capacities to tolerate intensity and recover. This further amplifies its impact on injury risk. Therefore, strength training programs must strictly control intensity and progression, prioritizing the development of foundational strength and movement stability. This approach helps mitigate injury risk caused by excessive or improperly planned training.

The Pathway 2 analysis of the CPDAG in Figure 4 shows that Strengthtraining(variable 6) indirectly influences injury risk through Perceived Perceivedtrainingsuccess(variable 7), which supports the conclusion that their impact on injury risk is relatively limited. In practical sports, this result suggests that strength training should pay particular attention to intensity management. Proper planning of load progression is key to avoiding sports injuries, especially during high-intensity strength training, where it is important to incorporate restorative training [such as extending Hoursalternative(variable 8) ]. kmmiddle(variable 3) and Hoursalternative(variable 8) can serve as transition phases for high-intensity training, helping athletes gradually adapt to higher loads and reduce fatigue accumulation. Combining feedback from Perceivedrecovery(variable 10) and Perceivedtrainingsuccess(variable 7) allows for the real-time optimization of training plans, ensuring a balance between performance enhancement and injury risk. Through careful design and dynamic adjustment of training plans, athletes can not only reduce the risk of training-induced injuries but also improve overall training quality. Additionally, we provide the causal direct effects of other variables in the Supplementary Figure 3 for further reference.