Exploring the effect of the triglyceride-glucose index on bone metabolism in prepubertal children, a retrospective study: insights from traditional methods and machine-learning-based bone remodeling prediction

Data selection and study design

A total of 356 prepubertal children who visited the Department of Pediatric Endocrinology at the Second Affiliated Hospital of Wenzhou Medical University between June 2021 and October 2023 were retrospectively analyzed. Of these, 332 participants met the inclusion criteria, comprising 163 boys and 169 girls, aged 7–11 years. To avoid the influence of puberty on the results of the study, we chose children in the Tanner stage 1 as the study participants (Eckert-Lind et al., 2020). This is because the rapid skeletal development of children during puberty is accompanied by an increased secretion of growth hormones and sex hormones (e.g., testosterone and estrogen), which play key roles in bone formation, remodeling, and metabolism. Written informed consent was obtained from both participants and their guardians. This study adhered to the principles of the Declaration of Helsinki and was approved by the Medical Ethics Committee of the Second Affiliated Hospital of Wenzhou Medical University (Approval No. 2024-K-031-02).

Inclusion criteria were as follows: (1) written informed consent was obtained from both participants and guardians, (2) participants cooperated with medical history taking, anthropometric measurements, and fasting blood sample collection, and (3) Tanner stage 1 children aged 7–11 years. Exclusion criteria included: (1) missing key data, (2) presence of metabolic bone disease, (3) chromosomal or genetic abnormalities affecting endocrine or bone metabolism, (4) history of bone fracture within the past year, (5) use of medications influencing bone metabolism, (6) Tanner stages 2–5, and (7) chronic systemic diseases.

Grouping and variable definitions

According to the standardized growth charts from the Centers for Disease Control and Prevention, adjusted for age and gender, we defined a healthy weight as a body mass index (BMI) between the 5th and 84th percentiles, overweight as a BMI between the 85th and 94th percentiles, and obesity as a BMI at or above the 95th percentile, as described by Peinado Fabregat, Saynina & Sanders (2023). We categorized the participants into three groups: 135 in the healthy weight group (64 boys and 71 girls), 86 in the overweight group (41 boys and 45 girls), and 111 in the obese group (58 boys and 53 girls).

Since there was no clear cut-off value for BTMs, we classified β-CTx (threshold: 488.345 pg/ml), T-P1NP (threshold: 374.000 ng/ml), and N-MID (threshold: 29.280 ng/ml) according to the median level of BTMs in the participants to suggest enhanced or inhibited bone metabolism and remodeling, respectively. The TyG index is calculated using the formula: ln[fasting triglycerides (mg/dL) × fasting glucose (mg/dL)/2] (Tao et al., 2022). Similar to the TyG index, the TyG-BMI index is a novel biomarker for IR, calculated by multiplying the TyG index by BMI (Huang et al., 2023). Tang et al. (2022) demonstrated that systemic immune inflammation (SII) is a reliable and stable hematological indicator of the body’s systemic immune and inflammatory status, calculated using the formula: platelet count × neutrophil count/lymphocyte count. The HOMA-IR was calculated using the formula: HOMA-IR = Fasting Insulin (μU/mL) × Fasting Glucose (mmol/L)/22.5 (Tahapary et al., 2022). Spexin, a recently identified adipokine belonging to the galanin/kisspeptin/spexin peptide family, plays a significant role in the pathophysiology of obesity-induced IR and T2DM (Fang et al., 2022). Fibroblast growth factor 23 (FGF23) is an osteogenic hormone produced and secreted by the skeletal system. Alongside parathyroid hormone (PTH), 1,25-dihydroxyvitamin D (1,25(OH)₂D), and calcitonin, FGF23 plays a critical role in regulating phosphate and calcium homeostasis across the bone, kidneys, and gastrointestinal tract to support optimal bone mineralization (Cipriani et al., 2022).

Missing value handles

In medical research, there are often missing values in medical data, which can result in wasted sample size if only complete variables are considered for analysis, and statistical analyses that ignore missing values have the potential to introduce bias in parameter estimation (Cao & Hu, 2024). Imputation of missing values plays a vital role in maintaining the integrity of the data and improving the validity of the study results. In this study, we used the R package “mice” to multiply impute variables with less than 15% missing data (Blazek et al., 2021). Variables with more than 15% missing data were excluded from the analysis.

Feature engineering and modeling strategies

Based on previous studies (Schini et al., 2023; Vasikaran et al., 2024) and expert opinions from pediatric and orthopedic specialists at the Second Affiliated Hospital of Wenzhou Medical University, we selected demographics, laboratory parameters, and anthropometric measurements as the initial model features for this study, comprising a total of 26 variables (Data S1). We screened the model variables through correlation analysis and near-zero variance tests to mitigate multicollinearity and overfitting. Based on the results of the correlation analysis (correlation coefficient threshold > 0.85) and the comparison of the explanatory power of the variables on the target variables, we removed four redundant variables: waist circumference, TyG-BMI index, insulin, and triglycerides, and ended up with a dataset containing 22 features and three target variables (β-CTx, T-P1NP, and N-MID). The dataset includes one categorical feature and 21 continuous features. To more robustly evaluate the performance of the model in unknown data, we randomly split the original predictive BTMs dataset into a training set and a test set based on 80:20. We used a five-fold cross-validation and grid search to determine the optimal combination of hyperparameters for the model, trained the model using the training set data, and evaluated the model’s performance in predicting BTMs in the test set data (Chen et al., 2023). The data were normalized using Scikit-Learn’s StandardScaler, and categorical variables were transformed into numerical form through one-hot encoding (Singh et al., 2020, 2023). Different ML algorithms have different applicability due to their ability to handle category variables, data nonlinearities, and computational efficacy. We employed five ML algorithms—support vector machine (SVM), Random Forest (RF), extreme gradient boosting (XGBoost), light gradient boosting machine (LGBM), and categorical boosting (CatBoost)—to investigate the relationship between TyG indexes and BTMs, including β-CTx, T-P1NP, and N-MID. After evaluating seven performance metrics for each ML model, the model that best captured the relationship between the TyG index and BTMs was selected and interpreted using SHapley Additive exPlanation (SHAP) analysis.

SHAP-based model interpretability approach

With the extensive application of ML and artificial intelligence in healthcare, model interpretability has emerged as a critical topic in both research and practical applications. The high complexity of these models, along with the nonlinearity of their features, often complicates the understanding of their internal decision logic and influencing factors. SHAP-based interpretation methods offer a novel approach for visualizing and quantifying model behavior, enhancing the transparency and usability of complex predictive models in healthcare settings (Fan et al., 2023). The SHAP value quantifies the importance of each feature by calculating its marginal contribution across various feature combinations, ensuring fair consideration of each feature’s role in different permutations (Wang et al., 2021). This approach enhances our understanding of how individual features influence model decisions. In the healthcare domain, SHAP values facilitate insights into how specific features impact a given patient’s diagnostic outcome, thereby supporting personalized diagnosis and treatment.

Statistical analysis

The study sample size was determined using the pmsampsize R package (Riley et al., 2020). Based on data distribution, normally distributed data were represented as mean ± standard deviation (SD), while non-normally distributed data were expressed as median and interquartile range (25%, 75%). The Kruskal-Wallis H test was used for comparisons across three groups, with categorical variables presented as frequencies (percentages) and analyzed using the Pearson chi-square test for group comparisons. Multivariate linear regression assessed the correlation between the TyG index and BTMs. The performance of the ML model was evaluated using metrics including the area under the receiver operating characteristic curve (AUROC), accuracy, precision, recall, F1 score, brier score, and the area under the P–R curve (AP). Statistical significance was defined as a two-sided P-value < 0.05. Data analysis was performed using R software version 4.3.1, with the following packages: tidyverse, dplyr, tsummary, ggplot2, tableone, glmnet, car, MASS, psych, rms, ggExtra, mice, pmsampsize, and ggpubr. ML analyses were conducted in Python version 3.11.5, utilizing the scikit-learn library (version 1.2.2).

Baseline characteristics of the study population

Figure 1 shows the process for inclusion and exclusion of participants, as well as the flow of the study design, which resulted in the inclusion of 332 prepubertal children as the study sample. Table 1 summarizes the baseline characteristics of participants grouped according to BMI. Gender, age, and 1,25(OH)₂D levels were not statistically different when compared between the healthy weight, overweight, and obese groups, respectively (P > 0.05). In the obese group, height, high-density lipoprotein (HDL), Spexin, β-CTx, T-P1NP, and N-MID levels were lower than in the healthy weight and overweight groups, while the rest of the variables were higher than in the healthy weight and overweight groups, and the difference was statistically significant (P < 0.05).

Dose-response effects of TyG index and BTMs

Table 2 summarizes the results of the multivariate linear regression analysis of the relationship between the TyG index and BTMs. Initially, the unadjusted model showed that the TyG index was associated with the response variables β-CTx (β = −75.98, 95% CI [−90.85 to −61.10], P < 0.001), T-P1NP (β = −91.99, 95% CI [−115.24 to −68.73], P < 0.001), and N-MID (β = −6.62, 95% CI [−8.54 to −4.70], P < 0.001), respectively, and was negatively correlated. In model 2, this association remained statistically significant even after adjusting for gender and age factors, β-CTx (β = −76.00, 95% CI [−90.90 to −61.11], P < 0.001), T-P1NP (β = −92.03, 95% CI [−115.23 to −68.83], P < 0.001), N-MID (β = −6.61, 95% CI [−8.52 to −4.70], P < 0.001). In Model 3, after adjusting for more covariates, for each unit increase in the TyG index, β-CTx decreased by 54.66 pg/ml (β = −54.66, 95% CI [−75.91 to −33.42], P < 0.001), T-P1NP decreased by 72.38 ng/ml (β = −72.38, 95% CI [−105.60 to −39.17], P < 0.001), and N-MID decreased by 4.21 ng/ml (β = −4.21, 95% CI [−6.97 to −1.46], P = 0.003), and the correlation remained significant. In Model 3, we observed that individuals in the highest quartile of the TyG index exhibited reductions in β-CTx, T-P1NP, and N-MID levels compared to those in the lowest quartile of the TyG index. Specifically, β-CTx decreased by 62.18 pg/mL (β = −62.18, 95% CI [−89.22 to −35.14], P < 0.001), T-P1NP decreased by 90.35 ng/mL (β = −90.35, 95% CI [−132.47 to −48.23], P < 0.001), and N-MID decreased by 5.07 ng/mL (β = −5.07, 95% CI [−8.56 to −1.58], P = 0.005). In model 3, there was a significant dose-response effect of the TyG index with β-CTx, T-P1NP, and N-MID when the TyG index was ≥8.74, ≥8.44, and ≥8.74, respectively (P < 0.05). Normality, independence, homoscedasticity, and linearity assumptions of the multivariate linear regression equations were evaluated using Q-Q plots, the Durbin-Watson test, the non-constant variance score test, and the variance inflation factor (VIF), respectively. Results indicate that all model assumptions are met, supporting the model’s robustness.

Comparison of ML model performance metrics

Table 3 summarizes fifteen ML models’ performance metrics for predicting BTMs. We trained fifteen ML models predicting β-CTx, T-P1NP, and N-MID based on the training set data and tested each model on the test set data for performance evaluation metrics including AUROC, accuracy, precision, recall, F1-score, Brier score, and AP, and 1,000 times resampling using Bootstrap to improve the robustness of the model predictions. Based on the confusion matrix and ROC curves, the CatBoost model has higher accuracy and robustness in predicting β-CTx, T-P1NP, and N-MID compared to other models (Fig. 2, Figs. S1–S4). The AUROC (95% CI) for predicting β- CTx, T-P1NP and N-MID were 0.782 [0.683–0.885], 0.789 [0.691–0.874], and 0.727 [0.619–0.827], respectively.

SHAP features importance and modeling decisions

The SHAP global summary plot shows the positive and negative impact of all features in the CatBoost model on predicting β-CTx, T-P1NP, and N-MID, sorted by the importance of each feature’s contribution to the model output. The model SHAP values suggested a negative correlation between TyG index and BMI for the prediction of BTMs (β-CTx, T-P1NP, N-MID) (Figs. 3A, 4A, and 5A). The SHAP summary plot suggests that the most important features for predicting β-CTx, T-P1NP, and N-MID in the CatBoost model are BMI, BMI, and HbA1c, respectively (Figs. 3B, 4B, and 5B). The SHAP individual decision plot reveals the individual prediction of BTMs (β-CTx, T-P1NP, N-MID) for each feature in the CatBoost model (Figs. 3C, 4C, and 5C). In the decision plot, the gray vertical line represents the base value of the CatBoost model’s prediction of BTMs, the red curve indicates a positive prediction of BTMs, and the blue curve indicates a negative prediction of BTMs. The decision curves reflect the specific contribution of each feature to the CatBoost model’s BTMs prediction and label the SHAP value of each feature.