Multimodel Interpretable Machine Learning for Osteoporosis Prediction in Diabetes Using a Large Critical Care Electronic Health Record Database

Summary

This protocol describes the development and evaluation of an interpretable ensemble tree-based machine learning model using routinely collected clinical data to predict osteoporosis risk in patients with diabetes.

Abstract

Diabetes mellitus is associated with increased skeletal fragility and elevated fracture risk; however, early identification of osteoporosis (OP) in diabetic populations remains challenging. This study aimed to develop and compare multiple machine-learning (ML) models for predicting OP among patients with diabetes, while evaluating their clinical utility and enhancing model transparency using SHapley Additive exPlanations (SHAP). Using routinely collected demographic and laboratory data from the MIMIC-IV database, multiple ML algorithms, including tree ensembles, gradient boosting, and linear baselines, were trained and validated. Model performance was comprehensively evaluated using discrimination metrics, including the area under the receiver operating characteristic curve (AUC) and precision–recall analysis, as well as probability calibration curves and decision curve analysis (DCA) to assess net clinical benefit. The best-performing model was further interpreted using SHAP to quantify and visualize feature contributions at both global and individual levels. Among all evaluated models, ensemble tree-based methods showed improved performance. The ExtraTrees classifier achieved the highest validation performance (AUC = 0.862; average precision = 0.866). Furthermore, the selected model exhibited excellent probability calibration (Brier score = 0.154) and demonstrated substantial net clinical benefit across a wide range of risk thresholds in the DCA. SHAP analysis identified gender and age as the most influential predictors, followed by routine hematologic and metabolic laboratory indicators. Local explanations provided clinically interpretable insights into individual predictions. An interpretable ensemble tree-based model effectively predicts OP risk in patients with diabetes using routinely available clinical variables. The integration of SHAP improves clinical transparency, and strong calibration coupled with positive net clinical benefit supports its potential implementation as a reliable decision-support tool for early OP risk stratification in diabetic populations.

Introduction

Diabetes mellitus and osteoporosis (OP) are highly prevalent chronic conditions that increasingly coexist, creating a substantial burden through fragility fractures, disability, and healthcare costs. In patients with type 2 diabetes (T2D), fracture risk is paradoxically elevated even when areal bone mineral density (BMD) measured by dual-energy X-ray absorptiometry (DXA) is preserved or higher than expected, highlighting that bone “quality” and systemic metabolic factors contribute importantly to skeletal fragility^1,2. Despite DXA being the clinical standard for diagnosing OP, its accessibility and uptake remain limited in many settings, and DXA-based assessment may undercapture diabetes-related skeletal risk. At the population level, meta-analytic evidence suggests that OP is common among individuals with diabetes, though prevalence estimates vary across cohorts and clinical contexts³. Therefore, there is a critical clinical need and strong motivation to develop pragmatic, scalable risk-prediction tools that leverage routinely collected clinical data to identify high-risk patients earlier, thereby guiding targeted screening and individualized prevention strategies.

In recent years, machine learning (ML) methods have been increasingly used to model complex, non-linear relationships in clinical data and to build OP risk predictors for patients with diabetes. Prior studies have demonstrated promising discrimination using algorithms such as gradient boosting, support vector machines, neural networks, and nomogram-style tools. Specifically, recent advancements have explored innovative modeling strategies, such as the fusion of ML with fuzzy logic, to predict OP and diabetes-related risks based on modifiable factors^4,5,6. However, a significant research gap remains: many existing models were developed in restricted subgroups (e.g., exclusively postmenopausal women or elderly T2D patients) and evaluated within a limited set of candidate algorithms, leaving uncertainty about which modeling choices are optimal for broader diabetic populations and heterogeneous electronic health record (EHR) data. Furthermore, a critical barrier to the clinical adoption of these high-performing ML models is their “black box” nature, which limits clinician trust and hinders translation into actionable decision support. SHapley Additive exPlanations (SHAP) provides a solution by offering case-level explanations, attributing each prediction to individual features using Shapley values, thereby enabling clinicians to verify whether the model’s reasoning aligns with clinical plausibility.

To address these gaps, we developed and validated an OP risk prediction framework for patients with diabetes using the Medical Information Mart for Intensive Care IV (MIMIC-IV) critical care EHR database⁷. The novelty of this study lies in the comprehensive benchmarking across multiple ML model families combined with transparent SHAP-based interpretability within a large, clinically complex diabetic cohort. Unlike studies centered on a single algorithm, our end-to-end approach systematically selects the best-performing model based on rigorous validation and provides transparent, patient-level explanations. This framework aims to deliver both robust predictive performance and clinically actionable insights, supporting earlier risk stratification for OP among patients with diabetes.

Protocol

Data source and ethics statement
This retrospective study used data from the MIMIC-IV (version 3.1), a publicly available, de-identified critical care EHR database hosted on PhysioNet. The database includes demographic information, vital signs, laboratory tests, diagnoses, procedures, and medications.

Access to MIMIC-IV requires credentialed authorization and adherence to the PhysioNet data use agreement and training requirements. The investigators completed the required training and were granted access (certificate number: 68351548). Because the database is de-identified, this study analyzed anonymized data and did not involve direct contact with patients; therefore, informed consent was not required. The complete, step-by-step executable modeling workflow is illustrated in Figure 1.

Figure 1. Overall modeling workflow for osteoporosis (OP) risk prediction in patients with diabetes. Schematic representation of the study pipeline, including cohort identification, data extraction, feature engineering, multi-model benchmarking, model selection based on validation performance, and SHAP-based interpretability of the final model. Please click here to view a larger version of this figure.

Study population and data extraction
Data were extracted from the Medical Information Mart for Intensive Care IV (MIMIC-IV) database using a database management tool. To ensure reproducibility, the exact SQL queries, including table names and filtering logic used for cohort retrieval, are provided in Supplementary File 1. Adult patients (≥18 years old) diagnosed with diabetes mellitus were identified using ICD-9 codes (249, 250) and ICD-10 codes (E08–E13). The primary outcome, OP, was defined using ICD-9 codes (733.x) and ICD-10 codes (M80, M81, M82). Within the initial eligible diabetic cohort (n = 46,226), 3,672 positive events (patients with OP) and 42,554 negative events (patients without OP) were identified. A detailed patient selection and attrition flowchart is presented in Figure 2.

Figure 2. Flowchart of patient selection and cohort construction. The flowchart illustrates the stepwise cohort derivation from the MIMIC-IV (v3.1) database. Adult patients with diabetes mellitus were identified using ICD codes, followed by exclusion of patients with age <18 years, missing critical data >30%, or invalid records. The final cohort was divided into OP (positive events) and non-OP (negative events). Class imbalance was addressed using random undersampling and SMOTE applied to the training data prior to stratified dataset splitting. Please click here to view a larger version of this figure.

Feature engineering
To ensure fair benchmarking across multiple ML algorithms and precise reproducibility, an identical sequence of preprocessing steps was applied: feature selection, missing value imputation, continuous feature scaling, class balancing, and dataset splitting. The full list of input features, including variable names, units, and data sources, is detailed in Table 1.

Features	Total (n = 14,688)	Training set (n = 9,546)	Validation set (n = 2,204)	Test set (n = 2,938)	P value
Gender					0.345
Male	6222 (42.36%)	4045 (42.37%)	935 (42.42%)	1242 (42.27%)
Female	8466 (57.64%)	5501 (57.63%)	1269 (57.58%)	1696 (57.73%)
Age	0.28 (−0.40, 0.94)	0.28 (−0.43, 0.94)	0.23 (−0.46, 0.89)	0.28 (−0.39, 0.95)	0.1655
Red Blood Cells	−0.14 (−0.79, 0.49)	−0.13 (−0.79, 0.49)	−0.17 (−0.83, 0.48)	−0.14 (−0.76, 0.49)	0.3641
Hematocrit	−0.14 (−0.81, 0.52)	−0.13 (−0.80, 0.52)	−0.14 (−0.87, 0.51)	−0.16 (−0.81, 0.52)	0.3709
Hemoglobin	−0.12 (−0.79, 0.52)	−0.12 (−0.79, 0.51)	−0.15 (−0.86, 0.53)	−0.13 (−0.78, 0.52)	0.3616
RDW	−0.16 (−0.59, 0.45)	−0.17 (−0.59, 0.46)	−0.14 (−0.59, 0.45)	−0.17 (−0.60, 0.42)	0.5803
Phosphate	−0.14 (−0.60, 0.38)	−0.15 (−0.61, 0.38)	−0.11 (−0.57, 0.38)	−0.13 (−0.56, 0.34)	0.415
MCV	0.05 (−0.50, 0.63)	0.06 (−0.50, 0.65)	0.03 (−0.49, 0.65)	0.02 (−0.50, 0.59)	0.5408
Magnesium	−0.10 (−0.52, 0.37)	−0.10 (−0.52, 0.37)	−0.10 (−0.52, 0.37)	−0.09 (−0.50, 0.37)	0.7797
White Blood Cells	−0.16 (−0.47, 0.21)	−0.16 (−0.47, 0.20)	−0.15 (−0.48, 0.25)	−0.16 (−0.46, 0.20)	0.6856
Anion Gap	−0.12 (−0.64, 0.45)	−0.13 (−0.64, 0.45)	−0.07 (−0.61, 0.48)	−0.13 (−0.64, 0.45)	0.1217
Creatinine	−0.30 (−0.48, 0.01)	−0.30 (−0.48, 0.01)	−0.30 (−0.46, 0.01)	−0.30 (−0.48, −0.00)	0.3323
MCH	0.11 (−0.48, 0.62)	0.12 (−0.48, 0.62)	0.11 (−0.47, 0.62)	0.09 (−0.48, 0.61)	0.5195
Platelet Count	−0.09 (−0.61, 0.49)	−0.09 (−0.61, 0.49)	−0.08 (−0.62, 0.47)	−0.10 (−0.62, 0.48)	0.9709
MCHC	−0.00 (−0.59, 0.59)	−0.00 (−0.60, 0.59)	−0.00 (−0.57, 0.58)	0.00 (−0.57, 0.60)	0.9263
Chloride	0.05 (−0.54, 0.64)	0.05 (−0.52, 0.65)	0.03 (−0.59, 0.62)	0.07 (−0.52, 0.62)	0.4675
Potassium	−0.07 (−0.67, 0.53)	−0.09 (−0.67, 0.53)	−0.07 (−0.67, 0.53)	−0.07 (−0.62, 0.53)	0.3071
Sodium	0.05 (−0.56, 0.60)	0.05 (−0.55, 0.60)	−0.00 (−0.60, 0.60)	0.07 (−0.57, 0.61)	0.3492
Urea Nitrogen	−0.28 (−0.60, 0.29)	−0.28 (−0.60, 0.29)	−0.26 (−0.59, 0.31)	−0.28 (−0.59, 0.25)	0.6826
Total Calcium	0.02 (−0.61, 0.59)	0.02 (−0.61, 0.59)	−0.01 (−0.59, 0.56)	0.01 (−0.61, 0.60)	0.8203

Table 1. Baseline characteristics and engineered features of the study cohort. Categorical variables are presented as n (%). Continuous variables are presented as median (interquartile range) of standardized values. P values were calculated to compare distributions across the training, validation, and test sets using appropriate statistical tests.

Continuous variables with repeated measurements were aggregated using the arithmetic mean within the entire ICU stay observation window to yield a single feature vector per patient. Variables with a missingness rate exceeding 30% were excluded. For the remaining numerical variables, missing values were imputed using the K-Nearest Neighbors (KNN) algorithm (n_neighbors = 5, weights = 'uniform'). Categorical variables were imputed using the most frequent category and subsequently transformed using binary mapping, with any unseen categories assigned a vector of zeros.

Continuous variables were standardized using the z-score scaling formula (). Features exhibiting zero variance were explicitly removed prior to scaling. All preprocessing parameters (µ and σ.) were estimated strictly on the training set and applied consistently to the validation and test sets.

Given the severe class imbalance (3,672 vs. 42,554), a hybrid balancing strategy was applied exclusively to the training data to avoid inflating performance estimates. First, random undersampling was used to reduce the majority negative class to achieve a 1:2 (positive:negative) ratio (yielding 7,344 negative samples). Next, the Synthetic Minority Over-sampling Technique (SMOTE) was applied to the positive class to achieve a final 1:1 balanced ratio (7,344 vs. 7,344)⁸.

Finally, the cohort was split at the patient level into training (65%), validation (15%), and test (20%) sets using stratified sampling. To strictly control all random processes across imputation, balancing, and splitting, a global random seed (random_state = 42) was enforced.

Model architecture
To identify the most suitable predictive model for OP risk in patients with diabetes, a large-scale benchmarking framework was used to compare 70 ML algorithm variants under an identical data representation and evaluation protocol. Candidate model families included regularized linear classifiers, support vector machines (SVMs), kNN, tree ensembles, gradient boosting architectures, and multilayer perceptrons^{9,10,11,12,13,14,15,16,17}.

Rather than using a traditional grid search, hyperparameter optimization was performed by evaluating predefined, robust configurations across these 70 model variants on the training split. Selection was strictly based on maximizing the area under the receiver operating characteristic (ROC) curve (AUC) on the validation set. The final best-performing model was configured with n_estimators = 200, with other parameters remaining at default settings of the software package listed in the Table of Materials. Evaluation metrics, including AUC, accuracy, F1-score, precision, and recall, were calculated using a default 0.5 probability threshold.

To support clinical transparency, SHAP was applied to the selected model using the TreeExplainer method. Given the tree-based nature of the final model, exact tree path-dependent expected values were used as the background configuration.

Results

Model performance comparison
To identify the optimal predictive approach for OP in patients with diabetes, multiple ML models were systematically compared. A regularized linear classifier (logistic regression) was included as a baseline control model. The overall performance comparison is summarized in Figures 3–6. Across all evaluated models, ensemble tree-based methods outperformed the linear baseline in terms of discriminative performance.

Discussion

In this study, multiple ML algorithms were systematically compared to identify an optimal predictive model for OP among patients with diabetes. Ensemble tree–based methods, particularly the ExtraTrees model, showed improved discrimination, reliable probability calibration, and greater net clinical benefit compared to linear baselines and other algorithm families. These findings indicate that nonlinear interactions and complex feature relationships are important for accurately capturing OP risk in diabetic populatio...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
CatBoost	Yandex	v1.2 (or applicable)	Gradient boosting algorithm optimized for categorical features
imbalanced-learn (SMOTE)	Scikit-learn Contrib	v0.11 (or applicable)	Library used to perform Synthetic Minority Over-sampling Technique for class balancing
LightGBM	Microsoft Corporation	v4.0 (or applicable)	Gradient boosting framework based on decision trees
Matplotlib	Matplotlib Development Team	v3.7 (or applicable)	Visualization library used for plotting figures and performance curves
MIMIC-IV Clinical Database	PhysioNet	v3.1	Publicly available critical care electronic health record dataset used as the data source
Navicat Premium	PremiumSoft CyberTech Ltd.	v17.0	Database management tool used for SQL-based data extraction
NumPy	NumPy Developers	v1.24 (or applicable)	Numerical computing library used for array operations and data processing
OpenPyXL	Eric Gazoni, Charlie Clark	v3.1 (or applicable)	Library used for reading and writing Excel files
Pandas	Pandas Development Team	v2.0 (or applicable)	Data manipulation and analysis library for structured datasets
PostgreSQL	PostgreSQL Global Development Group	v14 (or applicable)	Relational database system used to query and manage MIMIC-IV data
Python	Python Software Foundation	v3.8+ (or applicable)	Programming language used for data processing, modeling, and analysis
Scikit-learn	Scikit-learn Developers	v1.3 (or applicable)	Machine learning library used for model development, preprocessing, and evaluation
SciPy	SciPy Developers	v1.10 (or applicable)	Scientific computing library used for statistical analysis
Seaborn	Michael Waskom	v0.12 (or applicable)	Statistical data visualization library for enhanced graphical outputs
SHAP (SHapley Additive exPlanations)	Scott Lundberg	v0.42 (or applicable)	Framework used to interpret model predictions using feature attribution values
XGBoost	DMLC (Distributed Machine Learning Community)	v2.0 (or applicable)	Gradient boosting algorithm widely used for structured data modeling