INTRODUCTION
Osteoporosis is a skeletal disorder that compromises bone strength, increasing the risk of fractures in the spine, hip, and other sites. Clinically, it is defined by low bone mineral density (BMD), with a T-score below −2.5 as a diagnostic threshold for postmenopausal women and men over 50. 1 In Europe alone, osteoporosis affects 22 million women and 5.5 million men, contributing to 3.5 million fragility fractures annually. 2
Osteoporotic vertebral compression fractures (OVCF) are among the most common and debilitating complications of osteoporosis, leading to pain, disability, and increased mortality. 3 4 5 6 7 OVCF are associated with an increased incidence of morbidity, including back pain, height loss, deformity, disability, and mortality. 6 7 Multiple thoracic fractures can cause restrictive lung disease, while lumbar fractures can alter abdominal anatomy, resulting in digestive issues. 8 The severe physical, functional, and psychological impact of OVCFs highlights the need for accurate and early prediction. 9 10
Despite the strong correlation between declining BMD and the risk of OVCF, there is no standardized predictive model for OVCF. Multimodal machine learning, which integrates data from multiple sources, has superior predictive accuracy compared with unimodal approaches. 11 Data fusion, which combines different data types, enhances model robustness by capturing a broader range of risk factors. Common fusion strategies include early-, intermediate-, and late-stage data integration. 12 13 14 15 16
In this work, “multimodal” refers specifically to structured clinical variables (including DEXA-derived body composition) and spinal imaging. We acknowledge that this scope is narrower than broader multimodal paradigms combining, for example, text or genomic data; our focus is a practical clinical–imaging fusion relevant to spine care.
Our aims were to explore whether combining body-composition signals and imaging representations improves risk stratification for a second OVCF after an index OVCF and to identify interpretable features that may guide preventive care pathways.
MATERIALS AND METHODS
Study participants
This study utilized retrospectively collected data from patients at a single institution who sustained their first OVCF between January 1, 2000, and December 31, 2019. The study was approved by the Institutional Review Board of Ajou University Hospital (IRB No. AJOUIRB-DB-2023-030). The inclusion criteria were as follows: 1) post- or perimenopausal women and men aged ≥50 years; 2) DEXA T-scores ≤−2.5 representing the BMD of the femoral neck (BMD-F) or the lumbar spine (BMD-L); 3) initial images did not show OVCF, whereas follow-up data showed OVCF; and 4) during the follow-up period, patients were diagnosed with OVCF at least once.
The exclusion criteria were as follows: 1) severe disorders that cause osteoporosis and limit the use of osteoporosis medication (multiple myeloma, renal failure, rheumatoid arthritis, and gastrointestinal disease); 2) spinal tumors, severe trauma occurring before OVCF diagnosis or after follow-up; and 3) no history of OVCF.
Participants were divided into two groups: those who experienced a second OVCF and those who did not. In this study, “second OVCF” was defined as a new OVCF occurring at a vertebral level that was not fractured at the time of the index OVCF. The event had to be confirmed on follow-up radiographs or MRI and associated with new or worsened back pain. Progressive collapse, recollapse, or peri-fracture extension at the index fractured level were not classified as second OVCF. Patients who underwent vertebral augmentation (vertebroplasty or kyphoplasty) at the index level remained eligible for inclusion; however, imaging or morphological changes confined to the treated vertebra were not considered second OVCF. Only new fractures at previously intact vertebral levels were counted as meeting the endpoint definition.
Model inputs
Terminology BMD-L and BMD-F denote lumbar-spine and femoral-neck T-scores, respectively. BMD-L was defined as the mean lumbar spine T-score across evaluable vertebrae from L1–L4. Vertebrae affected by degenerative change, fracture, or hardware were excluded, and the mean of the remaining evaluable levels (minimum of two vertebrae) was used. If fewer than two lumbar vertebrae were evaluable, BMD-L was not computed, and the case was excluded from analyses requiring BMD-L.
1) General conditions: age, sex, BMD-L, BMD-F, body mass index (BMI), diabetes mellitus, hypertension (HTN). Smoking history was coded positive if the individual had smoked ≥100 cigarettes in their lifetime; otherwise, negative. Alcohol consumption history was coded positive if the individual reported drinking on ≥5 days per week for >6 consecutive months; otherwise, negative.
2) Each medication group was further divided into two groups: injection and oral groups. Teriparatide (Forsteo), denosumab, and other bisphosphonates were the injectable agents administered. The oral medications administered were bazedoxifene, raloxifene, and bisphosphonates. This approach was chosen to reflect clinically relevant differences in adherence and dosing patterns between daily oral and intermittently administered injectable agents and to avoid excessively sparse categories, given the limited number of patients treated with each specific injectable drug (teriparatide, denosumab, injectable bisphosphonates).
3) Duration of drug use was defined as usage for >1 year [medication exposure was captured by route (oral vs. injection) and duration (>1 year)]. We recognize that drug-class–level categorization (e.g., bisphosphonate subtype, denosumab, teriparatide) would be clinically preferable and plan to refine this in future work.
4) Dual-energy X-ray absorptiometry (DEXA) parameters: regional tissue and fat measures (trunk, android, and gynoid regions) were obtained by DEXA using a GE HealthCare Lunar Prodigy densitometer (GE HealthCare, Wauwatosa, WI, USA) with enCORE software version [17.3]. The trunk region was defined as the area between the upper horizontal border just below the chin and the lower boundary that was formed by an imaginary line passing through the hip joints. For the android region, the upper boundary was placed above the pelvic line of demarcation at a position equivalent to 20% of the distance between the pelvis and the femoral neck, whereas the lower boundary was at the top of the pelvic line of demarcation. The lateral boundaries comprised the lines for the arms in the normal position for a whole-body scan. The gynoid region was defined as the upper boundary, which was positioned at 1.5 times the height of the android region below the pelvic line, and the lower boundary was situated at twice the height of the android region. The lateral boundaries were defined by the outer leg lines of the demarcation ( Fig. 1 ). The trunk, android, and gynoid tissue weights were determined by the total muscle and fat mass at each location.
5) Image data: the primary imaging input was sagittal T1-weighted spinal MRI. When MRI was unavailable at the baseline visit, a lateral spine radiograph from the same visit was used instead. In the internal cohort, MRI was available for 139 patients (80.3%), and lateral radiographs were used for 34 patients (19.7%). The distribution of imaging modalities across the 5-fold cross-validation splits and in the external cohort is provided in Supplementary Table 1 (only online). All imaging inputs used for model development were baseline studies obtained after the index OVCF and before the occurrence of any second OVCF. For patients who subsequently sustained a second OVCF, the imaging used for prediction corresponded to this baseline visit, and the mean interval between baseline imaging and the second OVCF was 21.3 months. Only imaging acquired before the occurrence of a second OVCF was included; any post-event imaging was excluded from model training and evaluation. All images were screened for motion, distortion, and noise artifacts before model ingestion; studies failing quality control were excluded. Grayscale images were replicated across three channels to match the network input format.
External validation
To assess the generalizability of the model, external validation was performed using an independent dataset collected from a separate institution, which included patients who met the same inclusion/exclusion criteria but were not involved in the model development cohort. The external cohort included 83 patients who experienced their first OVCF between January 1, 2017, and December 31, 2024, and the same outcome measures and preprocessing protocols were applied.
Model development
The cornerstone of our multimodal model development was the successful integration of image data and clinical data, including general conditions, medication history, and DEXA parameters. We tested two intermediate-level fusion methods. A schematic of the two models is shown in Fig. 2 .
Fig. 2A shows an intermediate-level fusion of image and clinical inputs using a dense layer. Clinical data were preprocessed using StandardScaler, which transforms the data such that each feature has zero mean and a standard deviation of one. Image data were preprocessed to ensure compatibility as input for the pre-trained EfficientNetB0 model (resizing, scaling, color channel configuration). Since EfficientNetB0 requires RGB-format images, grayscale MRI and X-ray matrices were replicated across the color channels.
Next, EfficientNetB0, pretrained on ImageNet weights, was employed for feature extraction from the image data. This architecture was selected for its lightweight design, which minimizes the floating-point operations per second (FLOPs) and reduces the number of parameters, ensuring that the model is compact, efficient for training, and capable of fast inference. 17 After feature extraction, a max pooling operation was applied to the imaging data, and the resulting features were concatenated with the preprocessed clinical data. The fused data were then processed using dense layers (fully connected layers).
Fig. 2B illustrates an alternative intermediate-level fusion method that combines images and clinical inputs using principal component analysis (PCA) and classifiers. The preprocessing step and usage of EfficientNetB0 are identical to those of the previous approach; however, PCA is applied to the extracted image features to reduce their dimensions and extract meaningful features. After fusion with the clinical data, a separate classifier was used for prediction. Random forest, 18 CatBoost, 19 linear regression (LR), and support vector machine (SVM) 20 were tested. The number of PCA components was systematically varied from 10 to 100 to identify the optimal configuration for ensuring predictive accuracy.
Training and validation
We used 5-fold patient-level cross-validation. Within each fold, standardization and PCA were fit on the training split only and applied to validation/test splits to minimize data leakage. Given class imbalance (second OVCF 38.2%), we prioritized area under the receiver operating characteristic curve (AUROC), area under the precision-recall curve (AUPRC), and F1-score; no re-sampling was applied. As a sensitivity analysis of modality contribution, we additionally trained a clinical/DEXA-only random forest model using the same 5-fold cross-validation framework but excluding all image-derived features. To avoid information leakage, all preprocessing steps—including feature standardization, convoluted neural network (CNN) feature extraction, PCA fitting, and classifier training—were performed within each training fold only and then applied to the corresponding held-out fold. Cross-validation was conducted at the patient level, ensuring no overlap of images between folds. Performance metrics were computed exclusively from out-of-fold predictions pooled across folds.
Model evaluation and Interpretation
The performance of each model was assessed using seven evaluation metrics: the AUROC, the AUPRC, F1-score, accuracy, recall, precision, and specificity.
Among the evaluated metrics, AUROC and F1-score were prioritized for model selection. Performance was summarized using 5-fold cross-validation and reported as the mean±standard deviation across folds. To quantify uncertainty, 95% confidence intervals (CIs) for each metric were obtained by nonparametric bootstrapping of the pooled out-of-fold predictions from cross-validation (1000 resamples).
Model calibration was assessed using reliability (calibration) curves generated with the calibration_curve function in scikit-learn. Predicted probabilities were grouped into deciles, and the observed event rate in each bin was plotted against the mean predicted risk. Overall probability accuracy was summarized using the Brier score.
To account for potential overfitting, we additionally performed bootstrap internal validation with 1000 resamples following Harrell’s approach. For each bootstrap sample, the full modeling pipeline was refit and evaluated on both the bootstrap sample (apparent performance) and the original dataset (test performance) to estimate optimism. Optimism-adjusted AUROC and F1-scores were then computed by subtracting the estimated optimism from the apparent performance.
In this study, the best-performing model was further interpreted using feature importance and Shapley Additive Explanations (SHAP) analysis. 21 Based on game theory, we quantified the positive and negative contributions of each variable to the model. Through SHAP analysis, we assessed both the global contribution of each feature to the target and the local contributions of individual features.
In addition, to better understand how the image components [principal components (PCs)] contribute to classification, a heatmap visualization was generated to highlight the areas represented by the PC components. We sorted significant image components according to feature importance. The feature map output of the CNN was then multiplied by the corresponding image component to generate a weighted importance map for each spatial location. The weighted feature map was converted into a heatmap and overlaid onto the original image. This method provided valuable insight into salient image features and the CNN decision-making process.
Statistical analysis
Python 3.12.3 was used for model development and statistical analysis. TensorFlow, Pandas, and scikit-learn were used for the overall coding process. For statistical analysis of clinical data, we employed the χ 2 -square test for calculating p -values for categorical variables and the independent two-sample t-test for continuous variables. Odds ratios (ORs) for continuous variables are reported per 1-unit increase in T-score for BMD-L and BMD-F, and per 100 g increase for regional fat measures, unless otherwise specified. Because BMD-L and BMD-F capture different skeletal compartments and showed only weak mutual correlation, opposite directions of association are statistically plausible in multivariable models. Statistical significance was set at p <0.05.
RESULTS
Patient demographics
We analyzed 178 patients with an index OVCF; 68 (38.2%) experienced a second OVCF during follow-up (mean interval of 21.3 months after the first OVCF). For patients who developed a second OVCF, this interval corresponds to the time between baseline imaging (index-fracture visit) and the subsequent fracture event. The proportion of women was higher in the second-OVCF group than in the no-second-OVCF group (94.1% vs. 78.2%). Between-group comparisons identified significant differences in medication exposure, route of administration (oral vs. injection), and DEXA-derived body-composition measures (trunk/android/gynoid tissue and fat) (all p <0.001). Age, BMI, smoking, alcohol use, diabetes, and HTN were not significantly different between groups ( Table 1 ). Given the sex imbalance, associations involving regional fat distribution should be interpreted cautiously.
To characterize collinearity among covariates, we computed a Pearson correlation matrix including sex, BMI, BMD, and DEXA-derived body-composition metrics ( Supplementary Table 2 , only online). Associations between sex (binary) and continuous variables were assessed using point-biserial correlations and are reported in the same table. BMI was positivey correlated with trunk fat (r≈0.62), trunk tissue (r≈0.57), and gynoid fat (r≈0.31), indicating that higher BMI tracked with greater regional adiposity. Trunk fat showed strong correlations with android fat (r≈0.82) and gynoid fat (r≈0.73), and gynoid fat was also closely correlated with android fat (r≈0.70), reflecting substantial collinearity among adiposity measures. In contrast, BMD-L and BMD-F showed only weak correlations with both BMI and regional fat/tissue measures (|r| generally <0.30), suggesting that BMD and body-composition metrics captured largely complementary information in this cohort.
Model results
Dense layer fusion using CNN
Table 2 shows the results of dense layer fusion using CNN (DLF model) for the prediction of a second OVCF. As 5-fold cross-validation was used for reliable results, the mean and standard deviation are both reported.
PCA and classifier-based fusion using CNN
Performance across PCA component counts (10/30/50/100) and classifiers (random forest, CatBoost, LR, SVM) is summarized in Tables 3 , 4 , 5 , 6 and Fig. 3 . Random forest and CatBoost generally outperformed other classifiers across PCA settings; SVM showed inconsistent performance. The final model selected for external testing was PCF with 50 PCA components+random forest, yielding mean internal accuracy 0.979 (0.022), precision 0.986 (0.029), recall 0.965 (0.058), F1-score 0.969 (0.029), AUROC 0.998 (0.002), and AUPRC 0.997 (0.003) ( Table 3 ). AUROC values approaching 1.0 in some configurations likely reflect optimistic estimates given sample size and should be interpreted with caution.
In a complementary sensitivity analysis, a model using only clinical and DEXA-derived variables (no imaging inputs) achieved good discrimination (mean accuracy 0.9052±0.0755, precision 0.9857±0.0286, recall 0.7692±0.2024, F1-score 0.8462±0.1395, AUROC 0.9736±0.0303, AUPRC 0.9490±0.0429) ( Supplementary Table 3 , only online). Compared with this tabular baseline, the full multimodal PCA-fusion random forest model with 50 components showed higher recall and F1-score as well as improved AUROC and AUPRC, indicating that image-derived features contributed additional predictive signal beyond clinical and DEXA variables.
To evaluate robustness to imaging-modality heterogeneity, we conducted sensitivity analyses restricted to MRI-only and X-ray–only subsets. In the MRI-only subset, the multimodal PCA-fusion model maintained good discriminative performance (accuracy 0.9571±0.0416, precision 0.9818±0.0364, recall 0.9167±0.0745, F1-score 0.9470±0.0513, AUROC 0.9948±0.0066, AUPRC 0.9928±0.0091). These results are summarized in Supplementary Table 4 (only online). The X-ray–only subset was too small to yield stable estimates. Near-ceiling AUROC/AUPRC values observed in some internal cross-validation settings likely reflect optimistic estimates related to limited sample size and feature selection within folds and should be interpreted cautiously.
Internal calibration
Calibration analysis using reliability curves ( Supplementary Fig. 1 , only online) showed generally good agreement between predicted probabilities and observed event frequencies. The model slightly overestimated risk in the low–intermediate probability range (0.1–0.5) and underestimated risk at the highest predicted probabilities, where observed event rates approached 1.0. The Brier score based on pooled out-of-fold predicted probabilities was 0.0599, suggesting reasonable probability calibration in internal validation. Bootstrap internal validation with 1000 resamples demonstrated modest optimism. The optimism-adjusted AUROC was 0.999 (95% bootstrap CI 0.982–0.999), the optimism-adjusted AUPRC was 0.999 (95% bootstrap CI 0.984–0.999), and the optimism-adjusted F1-score was 0.983 (95% CI 0.952–0.995) ( Supplementary Table 5 , only online).
Results of external validation
External testing at an independent institution (n=83; identical inclusion/exclusion criteria) used the final PCF (50-component) random forest model trained on the full internal development set. Discriminative performance was preserved but lower than internal estimates (accuracy 0.904, precision 0.902, recall 0.902, F1-score 0.902, AUROC 0.923, AUPRC 0.911), consistent with dataset shift and reduced overfitting relative to cross-validation ( Table 7 ).
Model interpretation using feature importance analysis
Feature importance analysis for the final model is shown in Fig. 4 . The most influential variables were body composition metrics such as gynoid fat, gynoid tissue, trunk fat, android fat, and android tissue. Notably, quantitative image features also contributed substantially to the predictive performance of the final model, underscoring the value of incorporating radiological parameters.
Model interpretation using SHAP analysis
Summary plots of SHAP values for the final model are shown in Fig. 5 . The overall pattern of influential features was consistent with the feature-importance analysis: gynoid fat, gynoid tissue, trunk fat, android fat, android tissue, and medication were among the variables with the largest model-based contributions. In the SHAP plots, higher gynoid, trunk, and android fat values tended to push predictions toward second OVCF, whereas higher gynoid and android tissue tended to shift predictions toward no second OVCF. These SHAP patterns reflect model behavior and do not by themselves establish independent associations. Consistent with this, in a multivariable logistic regression including sex, age, BMI, BMD-L, BMD-F, android fat, and gynoid fat, only gynoid fat remained significantly associated with second OVCF after adjustment, whereas android fat did not ( Supplementary Table 6 , only online). The apparently opposite directions of association for BMD-L and BMD-F reflect differences in skeletal compartments rather than a coding error; both variables were entered as continuous T-scores with identical scaling.
Multivariable analysis of clinical and DEXA variables
Because model-derived importance does not establish independent or causal associations, we additionally performed a multivariable logistic regression including sex, age, BMI, BMD-L, BMD-F, android fat, and gynoid fat to assess whether regional fat measures remained associated with second OVCF after adjustment.
In this adjusted model, gynoid fat mass remained independently associated with second OVCF (OR 1.49 per 100 g increase, 95% CI 1.22–1.82, p <0.001), corresponding to an OR of 1.87 per 1-SD increase, indicating a moderate but clinically interpretable effect size, whereas android fat was not significant. BMD-L (OR 1.74, 95% CI 1.06–2.86, p =0.029) and BMD-F (OR 0.37, 95% CI 0.15–0.92, p =0.033) also showed independent associations with risk, while BMI demonstrated only a borderline trend p ≈0.052). These results suggest that, within the constraints of this cohort, gynoid fat and BMD measures contribute a unique signal beyond basic demographics, whereas some of the apparent contribution of android fat in the machine learning model likely reflects overlap with correlated fat-distribution variables. Full regression results are provided in Supplementary Table 6 (only online).
Model interpretation using pca component visualization
To enhance interpretability, we selected a representative patient image and generated a heatmap of a weighted feature map of image components. These heatmaps were overlaid onto sagittal spinal MRI and lateral spine radiographs to visualize the spatial influence of key image features. As shown in Fig. 6 , the weighted components consistently focused on the lumbar vertebral bodies, adjacent anterior paravertebral soft tissues, and overall spinal alignment. Several components highlighted the index-level abnormality and neighboring vertebrae, providing qualitative support for the contribution of imaging features to second OVCF risk classification.
DISCUSSION
OVCF constitutes a prevalent complication of osteoporosis that affects approximately 25% of postmenopausal women in the United States. 22 OVCF significantly reduces the quality of life for up to a year post-fracture and increases morbidity and mortality. 23 Consequently, accurate predictive models are critical for identifying high-risk patients and enabling early intervention. Although traditional research has focused on statistical risk factor analysis, 24 25 our study developed a multimodal machine learning model to predict secondary OVCF with high precision.
Clinical machine learning models are generally considered reliable when their AUROC and F1-scores exceed 0.8. 26 Our model demonstrated exceptional performance, with F1-score, precision, and recall all above 0.96, supporting its robustness in identifying patients at risk of recurrence. Nonetheless, these performance estimates should be interpreted cautiously given the modest, single-region sample and retrospective design. Although external testing showed preserved discrimination, metrics were lower than in cross-validation, consistent with optimism from model selection and residual overfitting. Bootstrap internal validation with 1000 resamples demonstrated only modest optimism after correction; however, the near-ceiling internal performance observed in this relatively small, single-center cohort suggests a residual risk of overfitting. Accordingly, the present results should be interpreted as exploratory and hypothesis-generating, and confirmation in larger, multi-center cohorts with external calibration and decision-curve analyses will be required before clinical translation.
A pragmatic strength of our work is the clinical–imaging fusion itself: integrating DEXA-derived body composition with image-based representations. In internal cross-validation, the fusion models showed higher discrimination than single-modality baselines, suggesting complementary signal between clinical and imaging inputs. A clinical/DEXA-only model already achieved good performance, but the addition of imaging features further improved recall, F1-score, and discrimination metrics. These findings suggest that imaging contributes complementary information rather than replacing clinical predictors. However, given the single-center design and limited sample size, these estimates should be interpreted as exploratory, and future multi-site studies with architecture-matched modality-specific models will be required to more precisely quantify the relative contribution of each data modality.
We do not claim a head-to-head superiority over prior CT-based approaches; the figures from Hu, et al. 27 reflect different cohorts, imaging protocols, and endpoints and are not directly comparable. In our pipeline, sagittal T1-weighted MRI was the primary imaging source (with lateral radiographs used only when MRI was unavailable), and the observed gains were driven predominantly by MRI features combined with DEXA measures. To integrate high-dimensional image features with tabular variables while mitigating overfitting, we applied PCA to the CNN embeddings before fusion. This reduced dimensionality and variance, simplified the downstream classifier, and improved stability and interpretability by retaining components most aligned with the predictive signal.
Rather than claiming novelty, our work extends prior literature by evaluating DEXA-derived body-composition features within a multimodal model specifically for secondary OVCF risk. While DEXA is primarily used to quantify bone mineral content, we explored additional regional measures (android/gynoid/trunk fat and tissue) as complementary predictors. In our cohort, these features— especially gynoid fat—were frequently ranked as important in model-based analyses (feature importance and SHAP). In a multivariable logistic regression including sex, age, BMI, BMD-L, BMD-F, android fat, and gynoid fat, only gynoid fat remained significantly associated with second OVCF after adjustment, whereas android fat did not. BMD-L and BMD-F also showed independent associations in this model. These adjusted findings are still exploratory and may be influenced by sample size and collinearity, but they suggest that lower-body fat distribution and bone density may provide complementary information beyond traditional BMD measures.
We note that univariable group differences were large, but we did not adjust p -values for multiple comparisons and therefore emphasize model-based signals over isolated hypothesis tests. Importantly, these findings do not establish causality or clinical cut-offs, and they do not contradict reports highlighting visceral adiposity or android-to-gynoid ratio; rather, they may reflect cohort composition (including sex distribution), imaging/DEXA protocols, and outcome definition (second vs. incident OVCF). Prospective, sex-stratified, multi-site studies will be needed to determine whether and how DEXA body-composition metrics should be operationalized alongside traditional risk factors (e.g., BMD) for secondary fracture prevention.
DEXA is a widely used modality for estimating whole-body and regional body composition because it separates bone from soft tissue and provides standardized android/gynoid/trunk metrics with very low radiation exposure (on the order of <0.1 µSv). 28 29 Rather than asserting superiority over CT or MRI, we note that DEXA offers a pragmatic, low-dose option for older adults and routine follow-up, whereas CT/MRI can provide finer anatomic detail but are less feasible for repeated body-composition assessment.
In our context, incorporating DEXA-derived parameters into a predictive model is a plausible adjunct for risk stratification, particularly in settings where advanced imaging is limited. Any workflow impact remains exploratory: prospective calibration (e.g., thresholds, calibration curves, decision-curve analysis) and multi-site validation are needed before clinical adoption. If validated, such models could flag patients for earlier preventive measures (optimization of osteoporosis pharmacotherapy; targeted exercise and fall prevention; nutrition counseling) while avoiding overuse of high-cost imaging.
While DEXA-derived body-composition metrics are not traditional primary risk factors (e.g., BMD T-scores) for incident OVCF, our question concerns recurrence after an index fracture. In this context, regional fat and tissue distributions may capture post-fracture physiology and frailty not reflected by BMD alone. Multi-site cohorts will be necessary to determine whether these signals persist after adjustment and to define clinically usable thresholds.
Accordingly, DEXA features are best viewed as complementary inputs within a multimodal framework, not replacements for established risk factors. If prospectively validated—with calibration and decision-curve analyses—such models could help flag higher-risk patients earlier to prompt guideline-concordant osteoporosis therapy optimization and targeted nonpharmacologic measures (exercise, nutrition, fall-prevention).
Limitations
This exploratory study has several limitations. First, the retrospective, single-region design and modest sample size—particularly for patients with a second OVCF—raise risks of selection and spectrum bias and limit precision. Although we used patient-level 5-fold cross-validation and an independent site for testing, the very high internal discrimination suggests potential overfitting, and generalizability beyond similar institutions remains uncertain. Second, sex distribution differed between groups, and subgroup analyses (e.g., sex-stratified effects or interactions) were underpowered; residual confounding by sex and other factors is possible.
Third, osteoporosis medications were coded pragmatically by route of administration (oral vs. injectable) rather than by drug class, dose, or cumulative exposure. Because the numbers of patients treated with individual injectable agents (teriparatide, denosumab, injectable bisphosphonates) were small, we grouped them into a single “injection” category to improve model stability. This route-based grouping does not distinguish anabolic from antiresorptive mechanisms and may obscure drug-specific associations; therefore, any findings related to injectable therapy should be regarded as exploratory and interpreted with caution. Moreover, granular data on chronic glucocorticoid exposure, fall history, and formal sarcopenia measures were not available, and residual confounding by these factors is likely.
Fourth, imaging inputs were heterogeneous: MRI served as the primary modality, whereas lateral radiographs were used when MRI was not available. Although MRI-only analyses showed preserved high performance, the radiograph subset was small, limiting reliable modality-specific conclusions. Future multi-site datasets with standardized imaging protocols will be needed to evaluate modality-specific effects more rigorously.
Lastly, we did not have complete and standardized information on several important fracture risk determinants, including chronic glucocorticoid exposure (dose and duration), history of falls, and formal measures of sarcopenia (e.g., grip strength, gait speed, or validated muscle mass indices). Although DEXA-derived lean mass and comorbidity profiles may partially reflect these constructs, they are not adequate substitutes, and residual confounding by these unmeasured factors may have influenced both exposure and outcome. Our model should therefore be viewed as a predictive tool rather than a causal framework.
Conclusion
In summary, a pragmatic clinical–imaging fusion model integrating DEXA-derived body-composition measures with image-based features showed strong internal discrimination and preserved performance at an external site, suggesting signal for predicting second OVCF after an index fracture. Given sample size, potential confounding, and limited external heterogeneity, these findings should be viewed as hypothesis-generating rather than practice-changing. Next steps include multi-site prospective validation with protocolized imaging, sex-stratified and confounder-adjusted analyses, and reporting of calibration and decision-curve analyses to define risk thresholds and net benefit. If validated, such models could augment—not replace—established risk factors, supporting earlier, targeted secondary-prevention strategies in routine care.
The authors have no potential conflicts of interest to disclose.
AUTHOR CONTRIBUTIONS: Conceptualization: Hoyeon Cho and Sung Hyeon Noh. Data curation: Sung Hyeon Noh. Formal analysis: Hoyeon Cho, Min Kyoon Joo, and Sung Hyeon Noh. Investigation: Hoyeon Cho, Min Kyoon Joo, and Sung Hyeon Noh. Methodology Hoyeon Cho and Sung Hyeon Noh. Project administration: Hoyeon Cho, Sung Hyeon Noh, Pyung Goo Cho, and Sang Hyun Kim. Resources: Sung Hyeon Noh. Software: Hoyeon Cho. Supervision: Sung Hyeon Noh, Pyung Goo Cho, and Sang Hyun Kim. Validation: Hoyeon Cho and Sung Hyeon Noh. Visualization: Hoyeon Cho and Min Kyoon Joo. Writing—original draft: Hoyeon Cho and Min Kyoon Joo. Writing—review & editing: Hoyeon Cho, Min Kyoon Joo, and Sung Hyeon Noh. Approval of final manuscript: all authors.
1 Lane NE Epidemiology, etiology, and diagnosis of osteoporosis Am J Obstet Gynecol 2006 194 2 Suppl S3 S11 16448873 10.1016/j.ajog.2005.08.047
2 Hernlund E Svedbom A Ivergård M Compston J Cooper C Stenmark J Osteoporosis in the European Union: medical management, epidemiology and economic burden. A report prepared in collaboration with the International Osteoporosis Foundation (IOF) and the European Federation of Pharmaceutical Industry Associations (EFPIA) Arch Osteoporos 2013 8 136 24113837 10.1007/s11657-013-0136-1 PMC3880487
3 Marshall D Johnell O Wedel H Meta-analysis of how well measures of bone mineral density predict occurrence of osteoporotic fractures BMJ 1996 312 1254 1259 8634613 10.1136/bmj.312.7041.1254 PMC2351094
4 European Prospective Osteoporosis Study (EPOS) Group O’Neill TW The relationship between bone density and incident vertebral fracture in men and women J Bone Miner Res 2002 17 2214 2221 12469915 10.1359/jbmr.2002.17.12.2214
5 Cummings SR Black DM Nevitt MC Browner W Cauley J Ensrud K Bone density at various sites for prediction of hip fractures Lancet 1993 341 72 75 8093403 10.1016/0140-6736(93)92555-8
6 Guzon-Illescas O Perez Fernandez E Crespí Villarias N Quirós Donate FJ Peña M Alonso-Blas C Mortality after osteoporotic hip fracture: incidence, trends, and associated factors J Orthop Surg Res 2019 14 203 31272470 10.1186/s13018-019-1226-6 PMC6610901
7 Miyakoshi N Itoi E Kobayashi M Kodama H Impact of postural deformities and spinal mobility on quality of life in postmenopausal osteoporosis Osteoporos Int 2003 14 1007 1012 14557854 10.1007/s00198-003-1510-4
8 Johnell O Kanis J Epidemiology of osteoporotic fractures Osteoporos Int 2005 16 Suppl 2 S3 S7 15365697 10.1007/s00198-004-1702-6
9 Lips P Cooper C Agnusdei D Caulin F Egger P Johnell O Quality of life in patients with vertebral fractures: validation of the Quality of Life Questionnaire of the European Foundation for Osteoporosis (QUALEFFO) Osteoporos Int 1999 10 150 160 10501796 10.1007/s001980050210
10 Ao J Xu Z Gao Z Ge T Wu J Li J What factors contribute to the poor prognosis of conservative treatment for osteoporotic vertebral compression fracture (OVCF): a systematic review Clin Interv Aging 2025 20 2717 2736 41480122 10.2147/CIA.S556710 PMC12755102
11 Kline A Wang H Li Y Dennis S Hutch M Xu Z Multimodal machine learning in precision health: a scoping review NPJ Digit Med 2022 5 171 36344814 10.1038/s41746-022-00712-8 PMC9640667
12 Castanedo F A review of data fusion techniques ScientificWorldJournal 2013 2013 704504 24288502 10.1155/2013/704504 PMC3826336
13 Steinberg AN Bowman CL White FE Revisions to the JDL data fusion model accessed on 2025 February 3 10.1117/12.341367
14 Molino-Minero-Re E Aguileta AA Brena RF Garcia-Ceja E Improved accuracy in predicting the best sensor fusion architecture for multiple domains Sensors (Basel) 2021 21 7007 34770318 10.3390/s21217007 PMC8586946
15 Blasch E Pham T Chong CY Koch W Leung H Braines D Machine learning/artificial intelligence for sensor data fusion–opportunities and challenges IEEE Aerosp Electron Syst Mag 2021 36 80 93
16 Lahat D Adali T Jutten C Multimodal data fusion: an overview of methods, challenges, and prospects Proc IEEE 2015 103 1449 1477
17 Tan M Le Q EfficientNet: rethinking model scaling for convolutional neural networks Proceedings of the 36th International Conference on Machine Learning 2019 June 9–15 Long Beach, CA, USA Proceedings of Machine Learning Research 2019 6105 6114
18 Breiman L Random forests Mach Learn 2001 45 5 32
19 Prokhorenkova L Gusev G Vorobev A Dorogush AV Gulin A CatBoost: unbiased boosting with categorical features Bengio S Wallach H Larochelle H Grauman K Cesa-Bianchi N Advances in Neural Information Processing Systems 31 (NeurIPS2018) Red Hook Curran Associates, Inc 2018 6639 6649
20 Tang Y Deep learning using linear support vector machines arXiv [Preprint] 2013 accessed on 2025 February 3 10.48550/arXiv.1306.0239
21 Lundberg SM Lee SI A unified approach to interpreting model predictions Guyon I von Luxburg U Bengio S Wallach H Fergus R Vishwanathan S Advances in Neural Information Processing Systems 30 (NeurIPS 2017) Red Hook Curran Associates, Inc. 2017 4768 4777
22 Greenberg M Handbook of neurosurgery 6th ed New York Thieme 2005
23 Jung HJ Park YS Seo HY Lee JC An KC Kim JH Quality of life in patients with osteoporotic vertebral compression fractures J Bone Metab 2017 24 187 196 28955695 10.11005/jbm.2017.24.3.187 PMC5613024
24 Hwang SH Cho PG Kim KT Kim KN Kim SH Noh SH What are the risk factors for a second osteoporotic vertebral compression fracture? Spine J 2023 23 1586 1592 37479141 10.1016/j.spinee.2023.07.010
25 Lee BG Choi JH Kim DY Choi WR Lee SG Kang CN Risk factors for newly developed osteoporotic vertebral compression fractures following treatment for osteoporotic vertebral compression fractures Spine J 2019 19 301 305 29959099 10.1016/j.spinee.2018.06.347
26 Çorbacıoğlu S¸K Aksel G Receiver operating characteristic curve analysis in diagnostic accuracy studies: a guide to interpreting the area under the curve value Turk J Emerg Med 2023 23 195 198 38024184 10.4103/tjem.tjem_182_23 PMC10664195
27 Hu X Zhu Y Qian Y Huang R Yin S Zeng Z Prediction of subsequent osteoporotic vertebral compression fracture on CT radiography via deep learning View 2022 3 20220012
28 Toombs RJ Ducher G Shepherd JA De Souza MJ The impact of recent technological advances on the trueness and precision of DXA to assess body composition Obesity (Silver Spring) 2012 20 30 39 21760631 10.1038/oby.2011.211
29 Damilakis J Adams JE Guglielmi G Link TM Radiation exposure in X-ray-based imaging techniques used in osteoporosis Eur Radiol 2010 20 2707 2714 20559834 10.1007/s00330-010-1845-0 PMC2948153
SUPPLEMENTARY MATERIALS
Supplementary Table 1
Distribution of Imaging Modalities (MRI vs. Lateral Radiograph) in the Internal and External Cohorts
Supplementary Table 2
Pearson correlation matrix among sex, BMD, BMI, and DEXA-derived body-composition metrics
Supplementary Table 3
Sensitivity Analysis: Model Performance in the clinical data only Subset (5-Fold Cross-Validation)
Supplementary Table 4
Sensitivity Analysis: Model Performance in the MRI-Only Subset, including clinical variables (5-Fold Cross-Validation)
Supplementary Table 5
Optimism-Adjusted Internal Validation via Bootstrap Resampling (1,000 Iterations)
Supplementary Table 6
Multivariable Logistic Regression Assessing Independent Associations of Clinical and DEXA-Derived Body-Composition Variables with Second OVCF
Supplementary Figure 1
Calibration curve of the final RF + PCA fusion model in the internal cohort.