*2.5. Analyses*

For descriptive statistics and summary tables, chronic pain (CP) was defined as pain severity (NPS) of greater than or equal to 4 at the second time point. The continuous underlying indices of NPS and pain interference were used in other analyses.

Demographic and clinical characteristics were summarized using counts with percentages for categorical variables and mean with standard deviation and range for continuous variables. These were compared between groups using the chi-squared test and Mann-Whitney test, respectively. The comparisons of CP between individual injury mechanisms used Fisher's exact tests with permutation-based adjustment over the possible mechanisms to control the overall type I error rate.

Based on initial bivariate exploratory analyses, the biomarkers (2-AG, AEA, cortisol) were log-transformed to improve the linearity of the relationships and reduce skewness. Pearson's correlation coefficient was used to quantify the strength of association between different biomarkers and between biomarkers and pain indices as continuous variables. The *p*-values were adjusted for multiple comparisons using the Benjamini–Hochberg method that controls the false discovery rate. For these exploratory analyses, FDR < 0.1 was considered a significant correlation.

A path model was developed to analyze the relationships among the circulating eCB and cortisol concentrations and pain measures at both time points, adjusting for covariates. The initial model structure was constructed based on biological plausibility, measurement timing, and the results of our exploratory correlational data. Specifically, the following variable groups were considered: sex and injury severity score (ISS); hospital and follow-up 2-AG and cortisol concentrations; self-reported pain score at hospitalization; and pain severity and interference scores obtained from the BPI at follow-up. In the initial model, sex and ISS were allowed to influence both hospitalization and follow-up pain measures. In addition, concentrations of 2-AG and cortisol were assumed to be correlated and could affect pain and biomarker measures at either time point. No direct effect of pain on 2-AG or cortisol was included. Hospitalization pain scores were included as a predictor of the follow-up pain measures. The effects of sex, cortisol, and 2-AG on pain were constrained to have equal strength at both hospitalization and follow-up, and the correlation between 2-AG and cortisol was also constrained to be the same at both time points. The model was fitted using the full information maximum likelihood method, which is a maximumlikelihood-based method that can incorporate missing-at-random observations [39].

The initial model was then simplified to find a more parsimonious description. Paths with non-significant effects with standardized coefficients under 0.1 in absolute value were removed, monitoring that goodness of fit indices continue to fall in their acceptable ranges and prioritizing models with lower Bayesian Information Criterion (BIC). Table 1 shows the goodness of fit indices of the initial and reduced models, indicating an excellent fit for both.

All analyses were performed using SAS 9.4 (SAS Institute, Cary, NC, USA) using the CALIS procedure for the path analysis. Unless otherwise noted, a two-sided 5% significance level was used.


**Table 1.** Goodness of fit indices of the initial and reduced models.
