Development and Validation of a Machine Learning Model for Predicting Osteoporosis in Patients with Parkinson’s Disease

Introduction

Osteoporosis is a systemic skeletal disorder characterized by reduced bone mass and impaired bone microarchitecture, which markedly increases the risk of fragility fractures.^1,2 In the general elderly population, osteoporosis represents a major public health burden, but in patients with Parkinson’s disease (PD) this risk is particularly pronounced.^3,4 Parkinson’s disease is a progressive neurodegenerative disorder characterized not only by motor dysfunction but also by non-motor comorbidities that increase skeletal fragility.^5,6 Factors such as impaired mobility, postural instability, frequent falls, reduced sunlight exposure, altered nutrition, and long-term medication use further contribute to the development of osteoporosis in these patients.^7–10 Several epidemiological studies have reported that the prevalence of osteoporosis and fracture incidence in PD patients is significantly higher than in age-matched controls, and fractures in this population often lead to prolonged hospitalization, loss of independence, and increased mortality.^4,11

Early identification of osteoporosis in PD patients is therefore essential for timely intervention and prevention of fractures. However, conventional clinical tools, such as the Fracture Risk Assessment Tool (FRAX), were designed for the general population and do not incorporate PD-specific risk factors.¹² Disease duration, Hoehn–Yahr stage, frequency of falls, use of walking aids, and exposure to medications such as corticosteroids, proton pump inhibitors, or selective serotonin reuptake inhibitors are relevant in the PD setting but are not adequately captured by traditional models.^13–16 As a result, clinicians often lack reliable methods to stratify osteoporosis risk in PD patients, which may delay the initiation of preventive therapy.

In recent years, machine learning (ML) has emerged as a promising approach for clinical risk prediction. Unlike conventional regression-based methods, ML can handle complex, nonlinear relationships and interactions among multiple variables, offering superior accuracy in high-dimensional datasets. ML algorithms have already demonstrated value in neurology, including prediction of cognitive decline, disease progression, and treatment outcomes in PD, as well as in metabolic bone diseases, where they have been used to predict bone mineral density, osteoporosis, and fracture risk.^17–19 Recent methodological work has also explored hybrid frameworks combining machine learning and fuzzy logic to improve osteoporosis risk prediction.²⁰ Despite these advances, few studies have specifically applied ML to predict osteoporosis in PD patients, and those that exist are limited by small sample sizes, single-center designs, or lack of external validation.

To address these gaps, we conducted a multicenter retrospective cohort study including 3935 PD patients across three tertiary hospitals in China. Our objective was to develop and externally validate machine learning models for predicting osteoporosis in PD, using a comprehensive set of demographics, lifestyle, clinical, and treatment-related variables. By incorporating PD-specific predictors and applying modern ML methods with external validation, we sought to establish a robust and interpretable prediction tool that could support individualized risk stratification and guide early preventive strategies in this vulnerable population.

Materials and Methods

Study Design

This was a multicenter, retrospective cohort study conducted at three tertiary teaching hospitals in China between January 2019 and December 2024. Two centers contributed data to the development cohort, and one center provided an independent external validation cohort. The study was designed and reported according to the TRIPOD (Transparent Reporting of a multivariable prediction model for Individual Prognosis or Diagnosis) statement.²¹ This study was approved by the institutional review board of the participating hospital, and the requirement for informed consent was waived due to the retrospective design and anonymization of patient data.

Study Population

We included patients with a confirmed diagnosis of PD and complete demographic, clinical, and treatment information. Patients with missing key variables, unclear PD diagnosis, or incomplete outcome data were excluded. After applying the inclusion and exclusion criteria, patients from two hospitals were used for model training and internal testing, while those from the third hospital were used for external validation. Figure 1 illustrates the flow of patient inclusion and exclusion across the participating centers.

Figure 1 Study profile and cohort construction.

Sample Size

The sample size was determined by the total number of eligible patients available in the participating centers during the study period. According to established recommendations for prediction model development, at least 10–20 outcome events per predictor variable (EPV) are required to ensure model stability.²² In our study, 3,935 patients with Parkinson’s disease were included, of whom 907 (23.1%) had osteoporosis. Fourteen predictors were retained after LASSO feature selection, corresponding to a minimum requirement of 140–280 events. The actual number of outcome events far exceeded this threshold, indicating that the study sample size was adequate for reliable model development and validation.

Outcome Definition

The primary outcome was osteoporosis, defined as a T-score ≤ −2.5 at the lumbar spine, femoral neck, or total hip, based on dual-energy X-ray absorptiometry (DXA). When multiple skeletal sites were available, the lowest T-score, defined as the most negative value across the lumbar spine, femoral neck, and total hip, was used for outcome determination, consistent with the diagnostic principle that osteoporosis can be diagnosed if any major site meets the threshold. All clinical and densitometric variables were collected at the time of hospital admission or outpatient evaluation. As the outcome was defined by baseline DXA T-scores, no longitudinal follow-up was required for outcome ascertainment.

Candidate Predictors

Based on prior literature, clinical expertise, and data availability across centers, we included a comprehensive set of candidate predictors for osteoporosis in patients with Parkinson’s disease. The variables comprised demographic and anthropometric factors including age, sex, and body mass index (BMI); lifestyle factors including smoking status, alcohol use, and regular physical activity; comorbidities such as hypertension, diabetes mellitus, dyslipidemia, chronic kidney disease, coronary artery disease, stroke, dementia, and depression; Parkinson’s disease–related characteristics including disease duration, Hoehn–Yahr stage, use of anti-Parkinson medications, and duration of PD treatment; medication exposures including use of anti-osteoporosis therapy, corticosteroids, proton pump inhibitors (PPI), and selective serotonin reuptake inhibitors (SSRI); and fracture risk–related variables including history of previous fracture, family history of osteoporosis, fall frequency, walking aid use, recent weight loss, and sarcopenia. The index date was defined as the date of the DXA scan. All candidate predictors were extracted from patients’ medical records within 6 months prior to the index date. Falls, recent weight loss, and medication exposures were operationalized using prespecified definitions based on medical records and prescription records within the six months look back window prior to the index date. Comorbidities, including sarcopenia, were identified based on clinician documented diagnoses and corresponding diagnostic codes recorded in the medical records.

Data Preprocessing and Missing Data

Raw data from each center were harmonized into a unified format. Continuous variables were converted to numeric and categorical variables were encoded with consistent factor levels. Implausible values were cross-checked against source records. Missing data were handled with multiple imputation by chained equations instead of case deletion. After imputation, the binary outcome was recoded as a two-level factor with “Yes” defined as the event.

Feature Selection with LASSO

We fitted binomial LASSO using glmnet with ten-fold cross-validation to identify the optimal penalty at the minimum cross-validated deviance (lambda.min). For stability assessment, we calculated the selection frequency of each predictor as the proportion of imputed datasets with a nonzero coefficient. A representative completed dataset was chosen whose lambda.min was closest to the median across imputations and was used to produce the coefficient path and cross-validation plots, and to enumerate predictors with nonzero coefficients at lambda.min.

Model Development and Internal Validation

Using the imputed development data, we split observations into a training set (70%) and a hold-out internal validation set (30%) by stratified sampling to preserve the event rate. We trained nine supervised learning algorithms under repeated cross-validation (ten folds, five repeats) with the area under the receiver operating characteristic curve (AUC) as the optimization metric and class probabilities enabled. The algorithms included support vector machine (SVM), gradient boosting machine (GBM), neural network, random forest, extreme gradient boosting (Xgboost), k-nearest neighbors (KNN), AdaBoost, LightGBM, and CatBoost.

External Validation

An independent dataset from a separate center was processed using the same numeric coercion and multiple-imputation strategy, with the outcome not imputed and any identifier excluded from the predictor matrix. To prevent leakage, we aligned the external features to the training schema. Absent one-hot encoded columns were added as zeros, and missing values were imputed using the same strategy as in the development cohort. We then generated model-specific predicted probabilities for the external cohort using the final models derived from the training set.

Performance Assessment

Discrimination was quantified by AUC with 95% confidence intervals from DeLong’s method and visualized with multi-model ROC curves. For each model we identified the Youden-optimal threshold on the relevant evaluation set and reported accuracy, sensitivity, specificity, precision, and F1 score. Calibration was assessed using nonparametric calibration curves constructed from predicted probabilities were grouped into deciles to construct calibration plots. Clinical utility was appraised with decision curve analysis under a cohort design with 50 bootstrap resamples, presented for the development, internal validation, and external validation datasets; when a center identifier was available in the external data, we additionally drew center-stratified decision curves to evaluate robustness across sites. The clinical decision was defined as whether to recommend DXA evaluation for osteoporosis assessment in patients with Parkinson’s disease. Treat all indicates recommending DXA to all patients and treat none indicates recommending DXA to no patients.

Model Explainability

For post-hoc interpretation, we applied KernelSHAP with shapviz to the best-performing models identified in the internal validation and external test cohorts. Background samples were drawn from the training data to approximate the underlying distribution, and a random subset of cases from the respective evaluation set was explained. We generated global beeswarm plots to illustrate the relative contribution and directionality of each predictor to the outcome.

Statistical Analysis

All analyses were performed using R (version 4.5.1). Continuous variables were summarized as medians with interquartile ranges and compared using the Mann–Whitney U-test, whereas categorical variables were presented as counts and percentages and compared using the χ²-test or Fisher’s exact test as appropriate. Missing data were handled using multiple imputation by chained equations (MICE) implemented in the mice package (version 3.19.0). Predictor selection was performed using LASSO penalized logistic regression implemented in the glmnet package (version 4.1.10). P values < 0.05 were considered statistically significant.

Results

Patient Characteristics

Among the 3,935 patients with Parkinson’s disease, 907 (23.1%) had osteoporosis. Compared with those without osteoporosis, these patients were older, had lower BMI, longer disease duration, and significantly reduced T-scores (all P<0.001). Women were more common in the osteoporosis group, and they were also more likely to have a history of fracture, family history of osteoporosis, and prior steroid use (all P<0.001). In terms of disease severity, advanced Hoehn–Yahr stage and higher fall frequency were markedly more prevalent among patients with osteoporosis. Other comorbidities and lifestyle factors showed no significant between-group differences. (Table 1)

Table 1 Baseline Characteristics of Patients with and without Osteoporosis

Feature Selection

LASSO regression with ten-fold cross-validation identified fourteen predictors with nonzero coefficients at the optimal penalty. The selected variables included Age, Sex, BMI, Smoking status, PD duration, Hoehn–Yahr stage, Anti-Parkinson medication use, Steroid use, PPI use, Previous fracture, Family history of osteoporosis, Fall frequency, Recent weight loss, and Sarcopenia. These variables were carried forward into model development (Figure 2).

Figure 2 LASSO regression for feature selection. (A) Coefficient profiles of the candidate variables plotted against the log(λ) sequence. Each colored line represents the trajectory of a variable coefficient as the penalty parameter λ varies. (B) Ten-fold cross-validation for tuning parameter selection in the LASSO model. The mean area under the receiver operating characteristic curve (AUC) and standard error are plotted versus log(λ). The dotted vertical lines indicate the optimal values of λ corresponding to the minimum mean cross-validated error (lambda.min) and the largest value within one standard error of the minimum (lambda.1se).

Model Development and Internal Validation

In the development cohort, all nine algorithms showed good discrimination with AUCs ranging from 0.881 to 0.935. LightGBM achieved the highest performance (AUC = 0.935, 95% CI: 0.923–0.947), followed by GBM (AUC = 0.922, 95% CI: 0.909–0.935). Random forest and AdaBoost performed similarly (both AUC = 0.910; 95% CI: 0.896–0.924 and 0.896–0.923, respectively). In the internal validation cohort, discrimination remained robust with AUCs between 0.809 and 0.869; the neural network performed best (AUC = 0.869, 95% CI: 0.843–0.895), followed by SVM (AUC = 0.866, 95% CI: 0.839–0.893) and GBM/AdaBoost (both AUC = 0.862, 95% CI: 0.834–0.889 and 0.835–0.889, respectively) (Figure 3).

Figure 3 Receiver operating characteristic (ROC) curves of nine machine-learning models for predicting osteoporosis in patients with Parkinson’s disease. (A) Training cohort. (B) Internal testing cohort. (C) External validation cohort.

External Validation

Applied to the independent external cohort, the models preserved strong discrimination with AUCs ranging from 0.855 to 0.937. The neural network achieved the highest AUC (0.937, 95% CI: 0.920–0.954), closely followed by SVM (0.935, 95% CI: 0.918–0.953) and random forest (0.935, 95% CI: 0.916–0.953). GBM and AdaBoost also performed well (AUC = 0.919, 95% CI: 0.898–0.941; and AUC = 0.922, 95% CI: 0.902–0.943, respectively), while CatBoost, XGBoost, KNN, and LightGBM maintained acceptable discrimination (AUCs = 0.916, 0.881, 0.858, and 0.855, respectively) consistent with the ROC curves and the threshold-based metrics summarized in Table 2.

Table 2 Threshold-Based Classification Performance in the External Validation Cohort

Model Calibration and Clinical Utility

Calibration analyses using decile-based plots showed close agreement between predicted risks and observed event rates in the development, internal testing, and external cohorts, with curves tracking the 45° line and only modest departures at the extremes (Figure 4).

Figure 4 Calibration and decision curve analyses across cohorts. (A–C) Decision curve analysis in the development cohort (A), internal testing cohort (B), and external validation cohort (C). Decision curves show net benefit for the strategy of recommending DXA evaluation at different risk thresholds. Treat all means recommending DXA for all patients, and treat none means recommending DXA for no patients. (D–F) Decile-based calibration plots in the development cohort (D), internal testing cohort (E), and external validation cohort (F).

Decision-curve analysis further demonstrated that random forest, LightGBM, and GBM provided higher net clinical benefit than “treat all” or “treat none” across clinically meaningful threshold probabilities, with consistent patterns replicated in both the testing and external validation cohorts (Figure 4).

Model Explainability

SSHAP analysis of the neural network model indicated that higher BMI was most strongly associated with a lower predicted risk of osteoporosis, whereas higher Hoehn–Yahr stage and longer PD duration were associated with a higher predicted risk. Sarcopenia and anti-Parkinson medication use also contributed notably, and other factors such as PPI use, smoking, recent weight loss, steroid use, previous fracture, and family history of osteoporosis showed consistent but smaller effects (Figure 5).

Figure 5 Global SHAP importance for the neural network model in the internal and external cohorts. (A) Beeswarm plot depicting SHAP value distributions for each feature in the internal cohort. (B) Corresponding SHAP beeswarm plot in the external validation cohort. Points are colored by feature value, and horizontal spread reflects the magnitude and direction of each feature’s contribution to the predicted probability.

Discussion

In this multicenter retrospective cohort study, we developed and externally validated a machine learning model to predict the risk of osteoporosis in patients with PD. By integrating routinely available clinical variables with PD-specific factors, the model demonstrated stable discrimination, good calibration, and favorable net benefit across clinically relevant thresholds. External validation in an independent center confirmed robust generalizability, while SHAP-based explainability analysis highlighted the major contributions of age, BMI, disease duration, Hoehn–Yahr stage, fall frequency, and medication exposure to individual predictions, thereby enhancing clinical interpretability.

The novelty of our study lies in being the first large-scale, multicenter effort to systematically incorporate PD-specific determinants into an ML framework for osteoporosis risk prediction. Unlike conventional tools developed for the general population, our model accounts for disease stage, functional impairment, fall frequency, use of walking aids, sarcopenia, recent weight loss, and exposure to bone-affecting medications such as corticosteroids, proton pump inhibitors, and SSRIs. Feature selection through LASSO followed by training across multiple algorithms ensured model parsimony and robustness. The resulting models achieved favorable AUCs in both development and validation cohorts, and calibration and decision curve analyses further supported their potential clinical value.

Previous studies in PD populations have mainly focused on epidemiological associations, reporting higher prevalence of osteoporosis and fracture compared with controls, or on identifying risk factors using traditional regression models.^23–26 However, these approaches have typically been limited by small sample sizes, single-center settings, or lack of independent validation. Compared with traditional risk assessment tools, our ML approach incorporates disease-specific variables and demonstrates improved discrimination. Although ML has been applied in osteoporosis and fracture risk prediction in general populations, studies specifically targeting PD patients, with external validation and model interpretability, remain scarce.^19,27 Our work fills this gap and shows that integrating PD-related predictors yields tangible improvements in risk stratification.

Clinical implementation of a PD-specific osteoporosis risk model requires attention to care setting, threshold selection, workflow integration, and the behavioral realities of shared decision-making. In ambulatory neurology and peri-hospital follow-up, individualized probabilities can function as low-cost pre-screening to prioritize DXA or earlier reassessment, addressing the consistent under-recognition of bone disease in PD despite a clearly elevated fracture burden. A recent review emphasized that standard assessment tools such as FRAX may underestimate the probability of major osteoporotic fractures in patients with Parkinson’s disease, partly because these tools fail to adequately account for PD-specific risk factors.²⁸ Emerging variants such as FRAXplus incorporate additional risk factors like recent falls, and QFracture includes PD as a named factor, but direct evaluations in PD populations remain limited.^12,29 A transportable, PD-tuned model can therefore complement existing calculators and help operationalize proactive screening pathways in movement-disorder clinics.

Decision-curve analysis offers a principled framework to identify threshold ranges where net benefit exceeds “treat all” or “treat none”, and to compare candidate algorithms on the dimension of clinical utility. Within those ranges, institutions can align thresholds with local DXA capacity, medication access, and fracture-liaison resources, while individual encounters can tailor choices to patient preferences regarding testing and pharmacotherapy burden. Clear reporting and interpretation of decision curves, as recommended by methodological primers and critiques, will be essential to avoid misapplication.^30–32

Embedding risk estimation into electronic health records may provide passive safety nets by generating contextual risk flags during neurology or geriatrics visits, issuing automated DXA orders for high-risk phenotypes, and prompting recalls after hospitalization for falls. Together, these functions create a closed loop of identification, testing, intervention, and re-evaluation. At the bedside, case-level explanations produced by post-hoc methods such as SHAP can make risk drivers legible to clinicians, patients, and caregivers, supporting concordant decisions around calcium and vitamin D optimization, initiation of anti-osteoporosis pharmacotherapy, and individualized fall-prevention programs. Recent syntheses of explainable AI in healthcare underscore both the promise of such explanations for trust and auditability and the need for careful use to avoid over-interpretation; pairing global importance summaries with simple, patient-facing narratives appears to be a practical compromise.^33–36

The clinical significance of a PD-specific model lies in its capacity to aggregate the cumulative effects of motor phenotypes, functional decline, body-composition changes, and medication exposures on skeletal fragility, translating diffuse risk into timely prevention. Implementation at scale could plausibly reduce preventable fractures, shorten hospital stays, and lower disability and mortality while optimizing the use of DXA scanners and fracture-liaison services. Recent neurology and movement-disorders reviews have also called for explicit incorporation of bone-health pathways into PD care guidelines, highlighting a system-level gap that predictive tools can help close.^28,37

Future work should extend beyond static baselines. Longitudinal updating that ingests routinely collected data, repeated weights and T-scores, medication changes, intercurrent falls—would allow dynamic risk trajectories and timely recalibration across centers. Integration of digital mobility measures from wearables and smartphones is particularly promising: gait speed, variability, turning metrics, and trunk motion derived from short in-clinic sensors or brief home monitoring have repeatedly associated with future falls and may refine fracture risk stratification when combined with traditional covariates. Parallel incorporation of biochemical markers, including vitamin D and bone-turnover indicators, and systematic assessment of sarcopenia and low lean mass could further tighten the mechanistic link between motor impairment, muscle-bone crosstalk, and skeletal fragility in PD.

Several limitations should be acknowledged. The retrospective design may leave residual confounding and selection bias. Despite harmonization efforts, measurement heterogeneity across centers cannot be fully eliminated. Osteoporosis was defined solely by T-scores, without systematically incorporating vertebral morphometric fractures or detailed fall severity indices, potentially underestimating clinically relevant outcomes. Moreover, the study population was restricted to three tertiary hospitals in China, limiting generalizability to other ethnic groups and healthcare systems. Prospective studies are required to assess transportability, recalibration needs, and temporal stability. Finally, although SHAP enhanced interpretability, consistency of explanations across algorithms and causal generalizability warrant further investigation.

Conclusion

In conclusion, we developed and externally validated a machine learning model for predicting osteoporosis in PD patients, incorporating disease stage, functional status, and medication exposures to achieve strong performance and clinical interpretability. Future work should focus on prospective, multi-regional validation and electronic health record integration to translate predictive accuracy into measurable improvements in screening efficiency, fracture prevention, and patient outcomes.