Osteoporotic vertebral fractures (OVFs) represent one of the most common fragility fractures in older adults and are associated with substantial morbidity, functional decline, and loss of independence. The affected population is typically characterized by advanced age, a female predominance, and multimorbidity, reflecting the epidemiology of vertebral fragility fractures in aging societies [1,2]. Beyond acute pain, OVFs frequently result in impaired mobility and reduced activities of daily living, with a considerable impact on discharge disposition and short-term recovery. Reliable prediction of discharge mobility from admission data would enable earlier targeting of rehabilitation resources, yet functional trajectories remain difficult to capture: individual clinical predictors show only a modest predictive signal, and recovery differs substantially between patients with ostensibly similar admission profiles.

Functional trajectories after hospitalization for OVF are heterogeneous. While fracture severity and treatment strategy guide therapeutic decision-making, growing evidence suggests that baseline mobility and frailty-related factors substantially influence postfracture functional outcomes [3-5]. However, functional reserve does not divide neatly into 2 states—gradations between mobility levels carry distinct implications for discharge planning, and modeling these distinctions as a multiclass outcome raises questions about class imbalance, threshold selection, and which prediction errors matter most clinically.

Prefracture mobility, comorbidity burden, and frailty have been associated with postfracture functional status [6]. Frailty in particular has been linked to increased disability and reduced health-related quality of life in patients with vertebral compression fractures. Nutritional status, rehabilitation timing, and mobilization strategies may further influence discharge outcomes in older fracture populations [7,8].

Machine learning models have been applied to predict discharge outcomes in spine surgery and fracture populations using routine clinical variables [9,10], including functional impairment at discharge in OVF [4], and comparable approaches have predicted residual symptoms following OVF treatment [11]. These models share a common design: a binary end point derived by dichotomizing functional scales, evaluated by aggregate metrics such as the area under the curve (AUC) [12]. This framing obscures a clinically relevant distinction—confusing a patient with high mobility with one who is nearly nonambulatory is a qualitatively different error than misclassifying adjacent mobility categories, yet existing models do not differentiate between these outcomes. Using the same data for both threshold selection and performance estimation is a recognized form of data leakage that leads to overoptimistic results [13,14], and, to our knowledge, no published OVF prediction model has explicitly incorporated a predefined decision policy for reducing clinically relevant high-consequence prediction errors.

The present study therefore modeled discharge mobility as a 3-level ordinal outcome (LOW, MID, and HIGH) and developed a 2-stage machine learning architecture in which a predefined decision policy governs the final 3-class assignment while reducing clinically relevant extreme misclassifications. Ordinal regression was used in parallel for clinical characterization of predictor associations. The primary methodological aim was to demonstrate that explicit error-type control is achievable in a real-world multiclass clinical prediction setting without sacrificing interpretability.

This prospective, single-center cohort study was conducted at a university hospital between June 2023 and December 2024. Consecutive patients aged 65 years or older who required inpatient treatment for an acute symptomatic thoracic or lumbar OVF were eligible. An osteoporotic fracture (OF) was defined as a fracture resulting from low-energy trauma or occurring without an adequate trauma mechanism.

Exclusion criteria comprised high-energy trauma, pathological fractures, spinal infections, or incomplete datasets regarding the primary outcome or core predictors. Treatment was performed according to routine clinical decision-making by the treating team [15,16].

The study protocol was approved by the institutional ethics committee from the Faculty of Medicine, Leipzig University (approval 214/23-ek), and the study was conducted in accordance with the Declaration of Helsinki. Written informed consent was obtained from all participants prior to enrollment.

The primary outcome was functional mobility at hospital discharge, categorized into 3 clinically meaningful groups: LOW (walking distance <50 m, transfer only, or bedridden), MID (walking distance 50‐500 m), and HIGH (walking distance >500 m or ability to climb stairs). These category boundaries were defined a priori by clinical consensus and were not derived from the study data. They were based on functionally meaningful walking distances described in previous mobility and rehabilitation research (approximately 50 m and 500 m, with stair climbing indicating higher functional mobility) [17,18]. No validated discharge-mobility classification specific to this population has been established. This 3-level ordinal classification served as the outcome variable for all analyses. Discharge mobility was assessed by the treating physiotherapists through direct observation and was routinely documented in standardized physiotherapy progress notes as part of clinical care.

Only variables available at admission or within −2 to +3 days relative to hospital admission were considered. Postoperative or in-hospital course variables were excluded.

Feature inclusion was determined a priori based on clinical reasoning and predefined availability criteria. The final feature set was fixed before model training and stored as a frozen feature matrix. No data-driven variable selection procedures were applied.

Continuous variables are reported as median with IQR. Categorical variables are presented as counts and percentages. Comparisons across discharge mobility categories (LOW, MID, and HIGH) were performed using Kruskal-Wallis tests for continuous variables and χ² tests or Fisher exact tests for categorical variables, as appropriate.

Regression analyses were conducted post hoc for clinical characterization and were independent of machine learning model development. Associations between selected clinical variables (age, CCI, OF grade, treatment modality, and dichotomized prefracture mobility) and discharge mobility were analyzed using proportional odds ordinal regression (OrderedModel). The proportional odds assumption was formally assessed prior to final model specification.

Additional binary logistic regression models (HIGH vs non-HIGH and LOW vs non-LOW) were estimated as supplementary analyses. Complete-case analysis was applied; as the regression model variables had no missing values, no cases were excluded.

A predefined 2-stage machine learning architecture was used to predict discharge mobility (Figure 1).

Stage A estimated the probability of HIGH vs non-HIGH mobility using a gradient-boosted decision tree model (HistGradientBoostingClassifier). Stage L estimated the probability of LOW vs non-LOW mobility using a random forest classifier (RandomForestClassifier), selected for its robustness to class imbalance and heterogeneous predictor effects. Missing values were imputed within the training folds of each cross-validation split, and the fitted imputers were applied to the held-out fold. Stage A used median imputation for continuous variables and most-frequent imputation with one-hot encoding for categorical variables; stage L used multivariable iterative imputation (multiple imputation by chained equations).

Structural hyperparameters for both stages were determined during a preceding model-development phase and were then applied unchanged across all 5 cross-validation folds. Stage A used fixed gradient-boosting hyperparameters without class weighting, and stage L used fixed random-forest hyperparameters with balanced class weights. No oversampling or undersampling procedures were applied to address class imbalance. No additional L1/L2 penalization was applied.

The 2 models output only class probabilities (P_high and P_low); the final 3-class label is assigned by a predefined decision policy.

Final 3-class predictions were derived from the 2 probability estimates (P_high and P_low). HIGH was assigned for P_high values greater than or equal to 0.63 or P_high values greater than or equal to 0.58 in patients with OF grade 4 or higher, operative treatment, and high prefracture mobility. LOW was assigned only for P_low values greater than or equal to 0.46 and P_high values less than or equal to 0.15; all remaining cases were assigned MID.

The decision-policy thresholds were determined during model development from the pooled out-of-fold predictions. Candidate configurations were restricted to those with a HIGH recall of at least 0.40, a LOW recall of at least 0.18, and a predicted MID proportion of at most 0.75. Among the remaining configurations, the one with the fewest major errors (LOW-HIGH) was selected, with macro-F₁, HIGH recall, and LOW recall as successive tie-breakers [20]. If no configuration met these constraints, the same ranking was applied without them. The absence of HIGH-to-LOW errors was not imposed as a constraint. The resulting thresholds were fixed and were not reestimated within each fold. The policy structure and this selection procedure were specified in advance, whereas the numerical thresholds were derived from the data as described.

Model performance was estimated from 5-fold stratified cross-validation using the out-of-fold predictions. Imputation and model fitting were performed within the training folds; the structural hyperparameters and decision-policy thresholds had been fixed during the preceding model-development phase. Major errors were defined a priori as misclassification between LOW and HIGH. Secondary metrics included class-wise recall and macro-F₁. All 31 candidate predictors were specified a priori on clinical grounds. The number of events was 70 (HIGH), 59 (MID), and 29 (LOW), corresponding to 2.26 and 0.94 events per candidate predictor for the HIGH and LOW stages, respectively [21].

To assess the added value of the 2-stage framework, we specified a baseline multinomial logistic regression model using the 3-level prefracture mobility category as the sole predictor (MID as the reference), with balanced class weights and L2 regularization (C=1.0), evaluated using 5-fold stratified cross-validation. Baseline predictions were assigned based on the maximum predicted probability; no threshold-based decision policy equivalent to that of the 2-stage framework was applied.

Calibration of both stages was assessed using the logistic calibration slope and intercept with 95% CIs, based on out-of-fold predictions, and displayed as calibration plots [22]. No post hoc probability recalibration was applied. Analyses were conducted in Python 3.12.7 using scikit-learn 1.8.0, pandas 3.0.0, NumPy 2.4.1, and Matplotlib 3.10.8, with fixed random seeds for reproducibility.

Because the thresholds were selected on the same out-of-fold predictions used for reporting, the decision policy was additionally evaluated in a nested design. Within each fold of an outer 5-fold cross-validation, stratified by the 3-level outcome so that both stages shared the same partition, both stages were refitted in an inner 5-fold cross-validation of the training partition. Thresholds were then selected from the inner out-of-fold predictions using the procedure described above and applied to the held-out outer fold. The entire procedure was repeated across ten 5-fold cross-validation resampling partitions, each covering all 158 patients. On the same partitions, the fixed thresholds of the final policy and the baseline model were evaluated for comparison. Structural hyperparameters remained fixed throughout and were not selected within the nested loops.

Feature importance was assessed for both stages using permutation importance computed on the held-out test fold of each cross-validation split and summarized across folds (30 repeats, receiver operating characteristic–AUC scoring; median with IQR). Permutation importance was preferred over impurity-based measures, which are biased toward continuous and high-cardinality variables. Importance measures were used for interpretation only and did not influence model specification, feature selection, or threshold determination.

The study included 158 patients aged 65 years or older with acute osteoporotic thoracolumbar fractures (Table 1). Overall missingness was low (44/4898, 0.9% feature cells) and confined to 6 laboratory variables, with vitamin D showing the highest variable-specific proportion of missing values (29/158, 18.4%). The proportion of operatively treated patients increased from 0% (0/5) in OF grade 1 to 96.2% (25/26) in OF grade 5 (P<.001).

Prefracture HIGH mobility was independently associated with higher discharge mobility (odds ratio [OR] 3.38, 95% CI 1.68-6.80; P<.001). In contrast, age, CCI, OF grade, and operative treatment were not significantly associated with discharge mobility (age per 10 y; OR 0.855, 95% CI 0.512-1.430; CCI OR 0.896, 95% CI 0.787-1.021 per point; OF grade OR 0.921, 95% CI 0.649-1.308; treatment OR 1.995, 95% CI 0.953-4.174).

Overall model fit was modest (Akaike information criterion [AIC] 317.2; McFadden R² 0.077), and the proportional odds assumption was met (Brant test P=.81; Figure 2). A null model restricted to prefracture mobility alone yielded McFadden R² 0.049 and AIC 318.5, indicating only marginal incremental fit from additional clinical predictors (ΔAIC=1.3).

Overall accuracy of the 2-stage model was 0.494, with balanced accuracy of 0.450 and macro-F₁ of 0.444. Detailed performance, discrimination, and calibration metrics are summarized in Table 2, and calibration plots are shown in Figure 3. Calibration slopes were below 1 for both stages, indicating imperfect calibration of the predicted probabilities. These values refer to the prespecified analysis using fixed thresholds; leakage-free estimates for the complete threshold-selection procedure are reported below.

Compared with the prefracture mobility–only baseline, the 2-stage framework showed higher macro-F₁ (0.444 vs 0.349) but only modest differences in discrimination, whereas calibration slopes were closer to 1 for the baseline (Table 3). MID recall was 0.73 vs 0.02 for the baseline, and major LOW-HIGH errors were 5 vs 22, including HIGH-to-LOW errors of 0 vs 7.

LOW mobility was correctly identified in 6 of 29 (20.7%) cases, with most remaining cases classified as MID (18/29, 62.1%). MID mobility was correctly identified in 43 of 59 (72.9%) cases, and HIGH in 29 of 70 (41.4%) cases. No HIGH cases were misclassified as LOW (Figure 4).

Observed and predicted discharge mobility distributions were examined across CCI strata, treatment categories, and predefined risk profiles. Across CCI strata, observed HIGH mobility decreased with increasing comorbidity (32/57, 56.1%; 29/61, 47.5%; and 9/40, 22.5% for CCI ≤6, 7‐9, and ≥10, respectively), whereas predicted HIGH proportions were consistently lower (26/57, 45.6%; 13/61, 21.3%; and 4/40, 10.0%). Similar patterns were observed across treatment groups and risk profiles, with predicted MID proportions exceeding observed MID proportions across all strata. Major errors (LOW-HIGH) remained limited across strata (CCI: 1/3/1; conservative: 1; operative: 4; high-risk: 1; low-risk: 0; Figure 5). Two predefined risk profiles were compared: a high-risk profile (high prefracture mobility, CCI ≥9, operative treatment; n=14) and a low-risk profile (high prefracture mobility, CCI≤6, conservative treatment; n=25).

All 5 major misclassifications represented LOW-to-HIGH transitions. Predicted probabilities for HIGH mobility were consistently elevated (0.721-0.829), placing these cases within the upper 16% (26/158) of the cohort’s P_high distribution. The median P_high among misclassified patients (0.782, IQR 0.741-0.794) exceeded that of all true HIGH patients (0.516, IQR 0.325-0.728), and 74% (52/70) true HIGH cases had a lower predicted probability than the lowest misclassified case. No HIGH-to-LOW misclassifications occurred.

Clinically, 3 patients had high prefracture mobility and underwent operative treatment for OF grade 4 or higher fractures. Two showed markedly elevated inflammatory markers with reduced hemoglobin, and 1, managed conservatively, was of advanced age with a high comorbidity burden. No single structural factor alone explained the discrepancy between predicted and observed mobility (Table 4).

Among true HIGH patients misclassified as MID, patients were older and had a higher comorbidity burden, lower hemoglobin, reduced estimated glomerular filtration rate (eGFR), and lower P_high probabilities compared with correctly classified true HIGH patients (81.9 vs 76.6 y; CCI 8.0 vs 5.0; hemoglobin 7.70 vs 8.50 mmol/L; eGFR 64.0 vs 84.0 mL/min; P_high 0.361 vs 0.766). Among true MID patients misclassified as LOW, elevated C-reactive protein (CRP), lower hemoglobin, reduced renal function, and higher P_low probabilities were observed (CRP 52.9 vs 16.3 mg/L; hemoglobin 6.50 vs 7.80 mmol/L; eGFR 44.0 vs 55.0 mL/min; P_low 0.580 vs 0.382). No consistent pattern was observed for LOW-to-MID or MID-to-HIGH misclassifications (Figure 6).

In stage A (HIGH vs non-HIGH), test-fold permutation importance was highest for aspartate aminotransferase/glutamic-oxaloacetic transaminase and eGFR, followed by age and prefracture mobility. Comorbidity burden (CCI) and CRP contributed moderately, while the remaining variables showed low and increasingly uncertain importance across folds (Figure 7).

In stage L (LOW vs non-LOW), hematocrit showed the highest importance, followed by several prefracture mobility variables and the metabolic markers vitamin D and calcium. Beyond these highest-ranked features, importance was low and showed wide fold-to-fold variability, reflecting limited and heterogeneous predictive signal (Figure 7).

Across 10 partitions, nested threshold selection yielded a macro-F₁ of 0.388 (SD 0.031), an accuracy of 0.457 (SD 0.038), and 8.3 (SD 3.4) major errors. These comprised 2.6 (SD 2.2) HIGH-to-LOW errors and 5.7 (SD 1.9) LOW-to-HIGH errors. No HIGH-to-LOW error occurred in 1 of 10 partitions.

Applying the fixed thresholds to the same partitions yielded a macro-F₁ of 0.390 (SD 0.031) and 8.4 (SD 3.1) major errors. The baseline model yielded a macro-F₁ of 0.350 (SD 0.009) and 22.5 (SD 1.1) major errors. The 2-stage framework produced fewer major errors in all 10 partitions (paired difference of −14.2, SD 3.1; Figure 8).

Class-wise recall under nested evaluation was 0.117 for LOW (SD 0.037), 0.669 for MID (SD 0.068), and 0.419 for HIGH (SD 0.045). Thresholds selected within folds were close to the fixed values (with a median of 0.63 for P_high and 0.48 for P_low). Results for each partition are provided in Multimedia Appendix 1.

Prefracture functional capacity was the dominant determinant of discharge mobility, irrespective of fracture grade and treatment strategy. Consistent with this, prefracture mobility was the only significant independent predictor in ordinal regression. In the machine-learning analyses, permutation importance ranked aspartate aminotransferase/glutamic-oxaloacetic transaminase and eGFR highest in stage A, followed by age and prefracture mobility, whereas hematocrit ranked highest in stage L, followed by several prefracture mobility variables. Two considerations qualify this apparent divergence between the regression and machine-learning rankings: the mobility signal entered the models through 4 correlated encodings (the ordinal scale and 3 binary indicators), a constellation known to dilute the permutation importance attributed to each individual feature [23,24], and importance estimates, particularly in stage L, showed wide fold-to-fold variability. Prefracture functional status therefore remained a consistent contributor across both analytic approaches, whereas the prominence of individual laboratory markers should not be read as evidence of biological relevance.

Using the predefined decision policy, no HIGH-to-LOW misclassification occurred in the prespecified analysis. Under nested cross-validation, such errors occurred in 9 of 10 partitions and averaged 2.6. The decision policy therefore reduces extreme misclassification rather than eliminating it, and the reduction relative to the baseline was present in all 10 partitions. This profile reflects the predefined decision policy rather than the probabilistic models alone. The 2 directions differ in clinical consequence: LOW-to-HIGH errors overestimate mobility and may result in insufficient rehabilitation and support after discharge, whereas HIGH-to-LOW errors overestimate the need for support. Both were less frequent than with the baseline (5.7 vs 15.0 and 2.6 vs 7.5). Five of 158 (3.2%) LOW-to-HIGH errors were observed; all showed elevated P_high values substantially exceeding the threshold, with predicted HIGH probabilities higher than those of most correctly classified HIGH patients. This pattern suggests that these cases presented with clinical profiles resembling higher-mobility patients at admission and that the discrepancy between predicted and actual outcomes was driven by in-hospital processes not captured by admission variables.

The favorable extreme-error profile was achieved in part through conservative assignment to the intermediate category: among the 158 patients, MID was predicted in 102 (64.6%) patients, compared with an observed prevalence of 59 (37.3%), with correspondingly low MID precision (0.42; Figure 4). This behavior is an intended property of the decision policy, which defaults to MID whenever neither extreme probability exceeds its threshold. The higher MID recall relative to the baseline should accordingly be understood as trading overassignment to the middle category for protection against extreme errors, not as improved discrimination of the MID class per se. The same behavior limits detection of LOW mobility: LOW recall was 0.207 in the prespecified analysis and 0.117 under nested evaluation, compared with 0.407 for the baseline. The baseline reached this recall by assigning LOW liberally, at the cost of 7.5 HIGH-to-LOW errors vs 2.6 for the 2-stage framework.

Our cohort corresponds closely to contemporary descriptions of hospitalized patients with OVF—advanced age, female predominance, and substantial comorbidity burden align with established epidemiological data [1] and prior OVF inpatient series [3,6,25]. The distribution of fracture severity and operative treatment reflects recent OF classification-based series [26,27].

The finding that prefracture functional capacity rather than fracture grade or treatment predicts discharge mobility is consistent with prior OVF data [3,4] and broader vertebral fragility fracture cohorts demonstrating that functional recovery is largely driven by preexisting mobility. In contrast, OVF-specific data linking individual laboratory parameters to early functional recovery remain limited [25], though evidence from orthogeriatric fragility fracture cohorts indicates that anemia, systemic illness, and impaired metabolic reserve are associated with poorer short-term functional trajectories [28-31].

From a modeling perspective, the moderate discrimination observed (AUCs of 0.693 and 0.625 for the binary components) is consistent with comparable spine surgery prediction studies [32,33] and with systematic evidence showing that machine learning does not consistently outperform regression when dominant baseline predictors capture the principal prognostic signal [34]. Major misclassifications occurred exclusively in the LOW-to-HIGH direction in the prespecified analysis, whereas HIGH-to-LOW errors also occurred under nested evaluation; minor misclassifications were confined to adjacent mobility states, consistent with the recognized difficulty of separating neighboring functional categories in geriatric fracture populations and with evidence that in-hospital processes substantially influence short-term recovery [35-37]. Taken together, the added complexity was not justified by gains in discrimination or calibration, which were modest or absent; its value lay in 3-class differentiation and the reduction of extreme misclassifications. Because the baseline assigned classes by the maximum probability, this reduction is attributable to the combination of the stage-wise probabilities with the predefined decision policy rather than to the machine-learning components alone. Whether this trade-off warrants the additional complexity depends on the clinical weight assigned to avoiding extreme errors. The extreme-error profile depended on the stability of the stage L probabilities rather than on the thresholds: with the out-of-fold probabilities held fixed, no HIGH-to-LOW error occurred, whereas refitting the models for each partition produced such errors irrespective of whether the thresholds were reselected.

This was a single-center cohort from a certified level 1 spine center, which may limit generalizability. The local care pathway typically involved inpatient geriatric rehabilitation after operative treatment, whereas conservatively treated patients were more often discharged home. Together with the center’s referral and admission profile of predominantly older patients with osteoporotic thoracolumbar fractures requiring inpatient treatment, these characteristics may have shaped the observed distribution of discharge mobility. Transportability therefore depends on comparable case mixes and rehabilitation pathways, and external validation in other settings is required before clinical implementation. Because the decision-policy thresholds were derived from the same pooled out-of-fold predictions that underlie the reported performance estimates, these estimates should be regarded as best-case internal results. The nested reevaluation quantified this optimism. On identical partitions, nested and fixed-threshold evaluations differed by 0.002 in macro-F₁, suggesting that the numerical threshold values were relatively stable across these partitions. Only the nested evaluation, however, provides a leakage-free estimate of the complete decision pipeline and is therefore the more realistic estimate of performance: macro-F₁ was 0.388 (SD 0.031) across 10 partitions, and HIGH-to-LOW errors were no longer absent. Discrimination showed the same pattern, with an AUC of 0.663 for stage A and 0.594 for stage L across partitions. Structural hyperparameters were selected on this cohort and were not nested, which remains an unquantified source of optimism. Prediction was restricted to admission variables and did not account for dynamic in-hospital processes, including acute medical complications, delirium, and rehabilitation response, that may substantially influence short-term recovery. This restriction was a deliberate design choice to enable prediction at admission; incorporating the in-hospital course might improve predictive performance but would address a different clinical question, placing an inherent ceiling on the performance achievable with admission-only data. The 3-level outcome definition increases modeling complexity compared with binary end points and may contribute to moderate discrimination. Unmeasured factors such as pain burden, rehabilitation intensity, and social support were not captured. The stage A calibration slope below 1 indicates some overfitting of the predicted probabilities; the stage L slope was estimated too imprecisely to be informative (95% CI 0.156‐1.600, including 1), given the limited number of LOW outcomes (n=29). These calibration properties suggest that the model is better suited to risk stratification than to precise individual probability estimation; external validation and, if necessary, recalibration would be required before probability-based clinical use. The moderate sample size relative to the 31 candidate predictors may have contributed to this instability [38], consistent with the uncertainty seen in the calibration and feature-importance estimates.

Strengths include the prospective design, standardized outcome assessment, clinically well-characterized real-world cohort, and the use of an ordinal mobility end point reflecting clinically relevant gradations. Model evaluation included detailed misclassification analysis distinguishing major from adjacent errors, enabling structured assessment beyond global performance metrics.

In hospitalized patients with osteoporotic thoracolumbar fractures, discharge mobility is primarily determined by prefracture functional reserve. The 2-stage framework achieved moderate multiclass discrimination and reduced extreme LOW-HIGH misclassification relative to a mobility-only baseline; this reduction persisted under nested cross-validation. Using routinely available variables at admission, the model may support early risk stratification, with potential applications in geriatric comanagement and physiotherapy resource planning. External validation in multicenter cohorts is a prerequisite before clinical implementation.