This retrospective cohort study enrolled 560 patients diagnosed with OVCF who underwent PVP between August 2015 and August 2024 at a tertiary medical center. Re-fractures were defined as radiologically confirmed new vertebral compression fractures at any level (adjacent or non-adjacent) occurring post-PVP, excluding fractures at the initially treated level. In this analysis, mental disorders included documented diagnoses such as depression, anxiety disorders, bipolar disorder, and schizophrenia, based on ICD-10 codes from electronic medical records. Inclusion criteria were: (1) confirmed fresh vertebral fracture via MRI (low T1 and high T2 signals) with localized tenderness and low back pain; (2) osteoporosis diagnosis based on dual-energy X-ray absorptiometry (DXA) or quantitative computed tomography (QCT) (T-score ≤ −2.5 or bone mineral density < 80 mg/cm³); (3) absence of pathological fractures (e.g., spinal tumors or infections) or high-energy trauma. Patients with prolonged bedridden status or incomplete follow-up data were excluded. The cohort was stratified into re-fracture and non-re-fracture groups based on postoperative imaging and clinical evaluations.

Data were extracted from electronic health records and standardized to ensure consistency. Missing values were addressed using predictive mean matching, a multiple imputation method preserving data distribution integrity (11). Continuous variables (e.g., age, bone density) were normalized via z-score transformation. Feature selection was performed using logistic regression (LR), which identified 14 non-zero coefficients by minimizing the binomial deviance through 10-fold cross-validation (12)

The dataset was randomly split into a training set (70%) and a testing set (30%). Six ML models were developed and evaluated: RF, LR, XGBoost, SVM, GBM, and MLP (13, 14). Model training and evaluation were conducted using Python (version 3.9) with libraries such as Scikit-learn, XGBoost, and SHAP.

Model performance was assessed using the AUC-ROC, accuracy, sensitivity, specificity, F1 score, and PR curves. Calibration curves evaluated prediction reliability, while DCA quantified clinical utility by calculating net benefit across threshold probabilities (0.0–1.0) (15). To enhance interpretability, SHAP values were computed to rank feature importance and visualize directional impacts on predictions (16).

The flow chart for the study was shown in Figure 1.

The study cohort comprised 560 patients with OVCF, divided into a training set (n = 392, 70%) and a testing set (n = 168, 30%). Baseline characteristics, including demographic and clinical variables, were well-balanced between the two sets, with no significant differences observed in most indicators (all p > 0.05, Table 1). The mean age of the cohort was 69.91 ± 6.77 years, with a slightly higher proportion of females (54.5%) compared to males (45.5%). Key comorbidities included hypertension (62.1%), diabetes mellitus (DM) (42.3%), chronic obstructive pulmonary disease (COPD, 54.8%), and CKD, 29.5%. Notably, the prevalence of re-fractures was consistent across the training (33.4%) and testing sets (33.9%, p = 0.984), ensuring comparable risk profiles for model development and validation. The balanced data partitioning, supported by p-values > 0.05 for all indicators, confirms the absence of significant bias between the two groups, validating the rationality and robustness of the dataset for ML analysis.

The LR analysis identified 14 characteristic variables of 22 variables. As illustrated in the coefficient trajectory plot (Figure 2A), which visualizes the regularization paths of variables during parameter tuning, increasing penalty parameters (λ) progressively eliminated non-contributory variables, retaining 14 features with non-zero coefficients. Cross-validation (Figure 2B) determined the optimal λ by minimizing binomial deviance, balancing model simplicity and predictive accuracy. Characteristic variables included scoliosis, mental status, CKD, trauma history, number of treated vertebrae ≤ 1 in the initial surgery (OP_lte_1), coronary heart disease (CHD), hypertension, DM, alcohol consumption ≥ 10 year (Alcohol_gte_10a), COPD, osteoarthritis (ost), coronary stent implantation, gout, tumor. These selected features were subsequently utilized for training and validating ML models, ensuring robust risk stratification while avoiding over fitting.

Six ML models were evaluated on an independent testing set (30% of the cohort), with the RF algorithm demonstrating superior performance across both training and testing phases. On the training set, RF achieved near-perfect discrimination (AUC = 0.99, 95% CI: 0.96–1.03; Figure 3A) and exceptional PR performance (AUC = 0.99; Figure 3E1), indicating robust learning without over fitting. Other models, including XGBoost (training set PR AUC = 0.96), GBM (0.95), and MLP (0.93), also showed strong performance, while SVM (0.90) and LR (0.82) lagged behind. This superiority extended to the test set, where RF maintained an AUC of 0.88 (95% CI: 0.83–0.93; Figure 3B), significantly outperforming LR (0.87), SVM (0.86), and XGBoost (0.87). Detailed performance metrics (Table 2) further validated RF’s dominance: on the training set, RF achieved the highest accuracy (0.95), sensitivity (0.98), and F1 score (0.93), while on the test set, it retained robust performance with accuracy (0.84), sensitivity (0.77), and specificity (0.87), surpassing SVM (accuracy = 0.78, F1 = 0.70) and LR (accuracy = 0.79, F1 = 0.74). DCA (Figure 3C) highlighted RF’s clinical utility, yielding the highest net benefit across threshold probabilities (0%–100%), while calibration curves (Figure 3D) confirmed strong alignment between predicted and observed outcomes (Brier score = 0.12). PR analysis on the test set (Figure 3E2) further reinforced RF’s reliability (AUC = 0.78), exceeding SVM (0.81) and GBM (0.82). Collectively, these results underscore RF’s consistency, generalizability, and clinical applicability in post-PVP re-fracture risk prediction.

The RF model demonstrated robust performance across multiple evaluation metrics, as detailed in Figure 4. On the training set, RF exhibited near-perfect discriminative ability with an AUC of 0.99 (95% CI: 0.96–1.03; Figure 4A), while maintaining strong generalizability to the test set (AUC = 0.88, 95% CI: 0.83–0.93; Figure 4B). PR analysis on the test set (Figure 4C) further validated RF’s reliability, achieving an AUC of 0.78, which, although slightly lower than SVM (0.81) and GBM (0.82), reflected its balanced performance in clinical risk stratification. DCA (Figure 4D) highlighted RF’s superior clinical utility, yielding the highest net benefit across threshold probabilities (0%–100%), particularly at the clinically relevant 20% threshold (net benefit = 0.6), outperforming alternative strategies. The Kolmogorov-Smirnov curve (Figure 4E) underscored RF’s discriminative power, with a KS statistic of 0.646, indicating clear separation between high- and low-risk patients. Finally, the confusion matrix (Figure 4F) quantified RF’s classification performance on the test set: 100 true negatives (specificity = 0.83), 36 true positives (sensitivity = 0.78), and an overall accuracy of 0.81, aligning with its robust calibration (Brier score = 0.12). Collectively, these results position RF as a comprehensive tool for post-PVP re-fracture risk prediction, excelling in both statistical rigor and clinical applicability 0.5.

SHAP analysis elucidated feature variables contributions to RF predictions (Figures 5A, B). The SHAP scatter plot (Figure 5A) further elucidated the relationship between feature values and their impact on the model’s output. Directional impacts were visualized through a SHAP scatter plot (Figure 5A), where binary features (e.g., “scoliosis_Yes”) exhibited strong positive associations with elevated risk (red dots, SHAP > 0). Conversely, absence of these features (blue dots, SHAP < 0) correlated with reduced risk. SHAP analysis revealed scoliosis (mean SHAP value = 0.14), impaired mental status (0.12), and CKD (0.10) as the top three predictors of re-fracture risk (Figure 5B). Secondary contributors included trauma history (SHAP = 0.08), severe osteoporosis (0.07), and coronary heart disease (0.06). Positive SHAP values for these features emphasizing the clinical relevance of these predictors (17).

While SVM and LR demonstrated comparable AUC performance (0.86–0.88) on the test set, the RF model exhibited significantly greater clinical net benefit (DCA curve area, Figure 3C), underscoring its superior decision-making utility for high-risk thresholds. This finding contrasts with Cai et al. (10), potentially due to their limited sample size (n = 385) compromising model stability. The RF model’s outperformance is consistent with established literature highlighting its efficacy in modeling complex, non-linear clinical data (18), particularly in capturing interactions among heterogeneous variables (e.g., comorbidities and biomechanical factors) often overlooked by linear models (19, 20). Notably, our results align with Xu et al. (21), who demonstrated RF’s robustness with imbalanced medical data—a frequent challenge in fracture risk prediction. The observed advantage of ensemble methods (RF/XGBoost) over simpler models (e.g., LR) reinforces the hypothesis that re-fracture risk is driven by multifactorial, non-additive interactions. Nevertheless, LR’s competitive AUC (0.88) implies that linear relationships may still dominate certain risk pathways, suggesting the need for future research into hybrid modeling strategies.

The model’s ability to provide interpretable predictions through SHAP values is a significant advancement over traditional “black-box” ML models. SHAP analysis not only identified the most important predictors but also quantified their impact on the model’s output, offering clinicians actionable insights (22). This level of interpretability is crucial for clinical adoption, as it allows healthcare providers to understand the rationale behind each prediction and tailor interventions accordingly. Similar approaches have been successfully applied in other medical domains, such as cardiovascular risk prediction (23) and cancer prognosis (24), further validating the utility of explainable ML models in healthcare.

Scoliosis has been identified as the strongest predictor of re-fractures in elderly patients with underlying diseases and OVCF (SHAP = 0.14). This result is supported by the previous studies (25–27) which demonstrating its significant biomechanical and metabolic impact. The abnormal spinal curvature disrupts load distribution, creating asymmetric stress concentrations that increase fracture susceptibility—evidenced by a postoperative Cobb angle ≥ 20° doubling the hazard ratio (HR = 6.243, p < 0.001) and finite element analyses highlighting uneven stress patterns in the “vertebral fractured arc” (T10–L4), where 93.6% of re-fractures occur. Additionally, scoliosis exacerbates osteoporosis progression through a bidirectional relationship: spinal malalignment accelerates bone loss via mechanical strain-induced osteoclast activation and impaired nutrient diffusion, leading to lower BMD at fracture sites (−3.7 vs. −3.2, p = 0.014) and further weakening structural integrity. This vicious cycle of biomechanical stress and bone fragility underscores scoliosis as a critical risk factor for post-surgical re-fractures, particularly after procedures like percutaneous kyphoplasty, where cement augmentation intensifies adjacent-segment stress. The identification of scoliosis as a top predictive factor in machine learning models highlights its critical role in vertebral re-fracture risk stratification. Clinically, this enables early intervention for high-risk patients—particularly elderly individuals with degenerative scoliosis—through multimodal approaches: biomechanical stabilization (bracing to correct load imbalance), osteoporosis management (antiresorptives, calcium/vitamin D supplementation), and nutritional optimization (protein-calorie support). Future interventions should combine these strategies with close monitoring of Cobb angle progression and BMD changes to disrupt the vicious cycle of spinal deformity and bone fragility.

Mental disorders (SHAP = 0.12) increased fracture risk may through neuroendocrine dysregulation, chronic inflammation, and oxidative stress, which impair bone remodeling and reduce bone density. Depression disrupts bone homeostasis through chronic inflammation, hypothalamic-pituitary-adrenal axis dysregulation, and oxidative stress, impairing osteoblast function and accelerating bone loss, as evidenced by a pooled hazard ratio of 1.24 for fractures in depressed individuals (28). Pharmacologically, selective serotonin reuptake inhibitors exacerbate fracture risk by inhibiting serotonin transporters in bone cells, suppressing osteoblast activity and enhancing osteoclastogenesis, with cohort studies showing adjusted HRs of 1.43 and 1.48 for major osteoporotic and hip fractures, respectively, even after adjusting for depression severity (29, 30). Notably, conventional risk assessment tools like FRAX underestimate fracture risk by 29%–36% in these populations due to the exclusion of mental health and psychotropic medication parameters, delaying critical interventions such as bone mineral density (BMD) monitoring or anti-resorptive therapies (30). To mitigate re-fracture risk, a multidisciplinary approach is essential: integrating mental health history into risk models, optimizing psychotropic prescriptions (e.g., favoring serotonin-norepinephrine reuptake inhibitors), addressing modifiable lifestyle risks, and prioritizing anti-osteoporotic therapies. Future studies should evaluate SSRIs as a potential mediator of fracture risk given their possible contribution to bone fragility, while also working to disentangle the independent contributions of mental disorders versus psychotropic medications, explore serotonin’s role in bone metabolism, and validate mental health-inclusive prediction models to refine preventive strategies for this vulnerable cohort.

Chronic kidney disease (SHAP value = 0.10) has been increasingly recognized as a significant risk factor for re-fracture following osteoporotic fractures, supported by converging clinical and epidemiological evidence. The association is multifaceted, involving direct skeletal alterations and systemic complications. First, CKD-induced mineral and bone disorders impair bone quality by disrupting calcium-phosphate homeostasis, leading to secondary hyperparathyroidism, abnormal bone turnover, and adynamic bone disease, which collectively reduce mechanical integrity. Additionally, uremia-induced oxidative stress and chronic inflammation accelerate bone resorption while suppressing osteoblast activity, as noted in Shimizu et al.’s (31) machine learning study, which ranked CKD as a top predictor of re-fracture due to dysregulated bone remodeling. Second, CKD patients often have comorbidities—such as cardiovascular disease, neuropathy, and muscle wasting—that synergistically increase fracture risk. Lourenço et al. (32) observed that CKD shortened the time to contralateral hip re-fracture, partly due to heightened fall propensity from frailty and uremia-related cognitive impairment. Third, suboptimal bone health management exacerbates risk; despite guidelines, Lin et al. (33) found only 6.8% of dialysis patients received anti-osteoporotic therapy, while Lourenço et al. (32) noted 78.9% of hip fracture patients (including CKD cases) lacked postoperative bone protection. This gap is critical, as CKD accelerates bone loss and complicates treatment (e.g., bisphosphonate contraindications in severe renal impairment). In conclusion, CKD drives re-fracture through bone metabolism disruption, fall-related risks, and systemic under treatment, necessitating integrated strategies targeting mineral homeostasis, fall prevention, and renal-adjusted osteoporosis therapy.

In addition to the top three risk factors, the model identified six other significant predictors: trauma history, OP_lte_1, CHD, hypertension, DM, and an Alcohol_gte_10a. Our SHAP analysis identified OP_lte_1 (≤ 1 surgical vertebra) as a significant predictor of re-fracture (SHAP = 0.07). This finding aligns with biomechanical studies demonstrating that single-level vertebroplasty increases adjacent-segment stress, whereas multi-level augmentation distributes forces more evenly (34, 35). Additionally, patients with OP_lte_1 may have untreated weak vertebrae, leading to subsequent fractures (27). Notably, multi-level cases (> 1 vertebra) often involve more aggressive surgical management, potentially masking their inherent risk (25). Thus, limited surgical intervention (≤ 1 vertebra) may serve as a marker for incomplete stabilization, warranting closer postoperative monitoring. Other factors are also supported by existing literature. For example, trauma history has been linked to increased fracture risk due to weakened bone structure (36). Hypertension and diabetes have been associated with bone loss and increased fracture risk due to their impact on bone metabolism (37). Alcohol consumption, particularly in long-term excessive alcohol intake, can cause alcohol-induced osteoporosis and increase fracture risk (38, 39).

Our RF model translates theoretical risk stratification into actionable perioperative pathways. Pre-operatively, it generates individualized risk scores using routine EHR data, enabling targeted interventions for high-risk patients (probability ≥ 30%): scoliosis stabilization consultations, psychiatric evaluations, and renal optimization prior to PVP. Postoperatively, risk-stratified monitoring tailors follow-up: high-risk patients receive 3 months clinical/imaging surveillance with osteoporosis/fall management, while low-risk patients (< 15%) follow standard 6 months schedules. SHAP interpretability (Figure 5) visualizes dominant risk drivers (e.g., scoliosis/CKD interaction), enabling personalized interventions beyond binary classification. We emphasize that this approach translates predictive net benefit (DCA) into modified care pathways—a pivotal advance toward real-world deployment.

Strengths: one of the key strengths of this study is the comprehensive evaluation of multiple ML models, which allowed for a robust comparison of their predictive performance. The use of SHAP values for model interpretability is another notable strength, as it provides clinicians with a clear understanding of the factors driving re-fracture risk. Additionally, the inclusion of a relatively large cohort of 560 patients, with a balanced distribution of re-fracture and non-re-fracture cases, enhances the generalizability of our findings.

However, limitations must be acknowledged. First, the retrospective design introduces potential selection bias, and unmeasured confounders (e.g., genetic predispositions, detailed lifestyle factors) were not included. Second, external validation in diverse populations is needed to confirm generalizability, as the cohort was derived from a single center. Third, while we adjusted for key comorbidities, unmeasured confounders—such as nutritional status (e.g., vitamin D/calcium levels), physical activity levels and genetic predispositions to osteoporosis—may influence re-fracture risk.