This study utilized two cross-sectional data sets, namely the National Health and Nutrition Examination Survey data sets of the United States (NHANES) and South Korea (KNHANES). NHANES is a cross-sectional study conducted by the National Center for Health Statistics to assess the overall health and nutritional status of the population in the United States [18].

Four cycles of NHANES from 2005 to 2006, 2007 to 2008, 2009 to 2010, and 2013 to 2014 were used, and we incorporated only respondents over the age of 50 years or those who had gone through menopause. Two cycles of NHANES (from 2011 to 2012 and from 2015 to 2016) were excluded because femoral neck BMD and total femur BMD data were not available for those two cycles. In NHANES, the total femur BMD and femoral neck BMD were measured using DXA, which was performed using a Hologic QDR-4500A fan-beam densitometer (Hologic, Inc; Bedford, MA, USA).

KNHANES is a cross-sectional survey of South Korean health and nutrition statuses with a representative population set derived through sampling that has been conducted by the Korean Centers for Disease Control and Prevention (KCDC) since 1998 [19]. KNHANES data from 2008 to 2011 were utilized in this study. The data sets from 2008 to 2009 consisted of 9200 households and the data sets from 2010 to 2011 consisted of 7680 households. In KNHANES, femoral neck BMD and total femur BMD were investigated in respondents aged over 50 years or among those who had gone through menopause, which were measured using DXA performed with a Hologic DISCOVERY QDR-4500W device (Hologic, Inc).

The methods were performed in accordance with relevant guidelines and regulations and approved by the Research Ethics Review Board of the National Center for Health Statistics (Protocol #2005-06, Protocol #2011-17). NHANES has obtained written informed consent from all respondents. KNHANES was approved by the KCDC Institutional Review Board (2008-04EXP-01-C, 2009-01CON-03-2C, 2010-02CON-21-C, 2011-02CON-06-C) and written consent was obtained from all participants.

For both NHANES and KNHANES, two different criteria based on the femoral neck BMD and total femur BMD were used to classify respondents into osteoporosis, osteopenia, and normal groups. Respondents with a T-score of femoral neck BMD or total femur BMD under −2.5 were defined as having osteoporosis. Those with T-scores between −1.0 and −2.5 were classified into the osteopenia group. The group with T-scores above −1.0 was defined as normal. For KNHANES, the T-scores were calculated using a reference group of Japanese individuals aged over 20 years [20] separately for the femoral neck and total femur. In the case of NHANES, we calculated the T-scores of respondents using a reference group consisting of non-Hispanic white women aged 20 to 29 years from NHANES III according to World Health Organization recommendations [21]. Using the calculated T-scores, we classified respondents into three groups: osteoporosis, osteopenia, and normal.

We only included respondents who had femoral neck and total femoral BMD records. Additionally, multiple variables with relevant information were merged into single variables or all but one variable were deleted. For example, nine variables indicating signs and symptoms of depression used in the Patient Health Questionnaire-9 (PHQ-9) [22] in NHANES were combined into one variable according to the PHQ-9 criteria. Multiple variables related to smoking, drinking, and blood pressure were eliminated, and we included only representative variables, as discussed previously [23,24]. Additionally, several variables with the same measurements but different units were removed, leaving only one variable. For KNHANES, several variables related to smoking, income, education level, drinking, and diabetes were merged or eliminated, following the data preprocessing methods described previously [25].

After merging variables, we excluded variables with nonnumeric values and considered responses of “refused” or “do not know” as missing values. Respondents who had missing values for more than 10% of all variables or variables that had missing values for more than 10% of all respondents were excluded. Missing values were filled in using k-nearest neighbors (KNN) imputation. At the end of preprocessing, a total of four types of data sets were created: two from NHANES and two from KNHANES, each using femoral neck BMD or total femur BMD when categorizing respondents into osteoporosis, osteopenia, or normal groups. The two data sets from NHANES comprised 8274 respondents with 89 variables and those from KNHANES comprised 8680 respondents with 162 variables. A flow diagram of the overall data processing is presented in Figure 1.

We implemented a DL algorithm and compared the results to those of ML algorithms. The DL model we used consisted of three layers: two dense layers with rectified linear unit activation functions and a final layer with a softmax activation function. Batch normalization and dropout layers were included after the two dense layers. Through five-fold cross-validation, the hyperparameters of the DL models were optimized based on five subsets of the created data sets. The optimized models for NHANES had two dense layers (each with 128 and 16 nodes for the femoral neck, and 128 and 64 nodes for the total femur) with a dropout rate of 0.2 using the Adam optimizer with a 0.005 learning rate. The models for KNHANES had two dense layers (each with 128 and 16 nodes for the femoral neck, and 128 and 32 nodes for the total femur) and used dropout rates of 0.2 (for femoral neck) and 0.4 (for total femur) with the Adam optimizer and a 0.005 learning rate. The models were trained and optimized in Keras using a TensorFlow backend (Google Inc) [26]. For the ML methods, seven models were considered for comparison: a nonlinear support vector machine, decision trees, extra trees, light gradient boosting machine (LGBM) classifier, logistic regression, KNN, and multilayer perceptron (MLP). We set the hyperparameters of the ML models to the default values within the scikit-learn package [27].

After training the DL and ML models, we applied an algorithm to explain the DL model’s behavior and rank the features in order of importance. For the DL models, local interpretable model-agnostic explanations (LIME) [28] was adopted as an XAI technique. LIME is an algorithm that explains classifier predictions by perturbing inputs and performing searches using local interpretability. The output of LIME represents the relative importance of the variables for classification. Therefore, we analyzed the results in two manners. First, we added the resulting values from LIME for each variable and ranked contributions according to their values. This rank represents the order of the variables based on their impact on the DL model’s classification of osteoporosis. Second, we observed the positive or negative effects of each feature on the probability of osteoporosis with respect to the numerical or categorical values of the features.

For a similar reason, the Boruta [29] and least absolute shrinkage and selection operator (LASSO) [30] methods were applied to the ML model with the best performance to calculate the importance values of features. We utilized the measured relative feature importance from Boruta and coefficients from the output of LASSO to rank the contributions of features. Boruta was used with the 80th percentile of shadow feature importance and LASSO was used with an α value of 0.001, which is the weight given to the sum of the squares of coefficients.

Widely known clinical assessment tools for diagnosing osteoporosis, such as OST, ORAI, and OSIRIS, were used for comparisons to evaluate the performance of the DL models. We compared the performance of our DL model with that of the clinical osteoporosis assessment tools on both the NHANES and KNHANES data sets. When classifying osteoporosis, OST was used for both men and women, whereas ORAI and OSIRIS were applied only to the female data set.

Comparisons of the continuous variables from the NHANES and KNHANES data sets were performed using the Student t-test and are presented as mean (SD) values. Categorical variables were compared using the χ² test and are presented as counts and percentages. A P value less than or equal to .001 was considered statistically significant. CIs were computed using the Student t-distribution.

The baseline characteristics of the respondents from the preprocessed data sets are presented in Table 1. Among the 8274 respondents in NHANES, the average age was approximately 65 years and 52% were men. The KNHANES respondents exhibited a similar distribution with a male ratio of 45% and average age of approximately 64 years. For both NHANES and KNHANES, more respondents were classified into the osteoporosis or osteopenia groups when femoral neck BMD was used instead of total femur BMD.

The receiver operating characteristic curves and corresponding AUC values of the optimized DL and ML models are presented in Figure 2 and Figure 3, respectively. The results were averaged across five-fold cross-validation. The DL model using NHANES exhibited a lower microaveraged AUC for the femoral neck BMD data set than for the total femur BMD data set (Table 2). Regarding individual classes, the classification of osteoporosis exhibited a higher AUC compared to the osteopenia and normal classes in the NHANES data set (Table 2). Similarly, for KNHANES, we achieved a higher AUC for the total femur data set than for the femoral neck BMD data set (Table 2). Additionally, the AUC for classifying osteoporosis was the highest in the KNHANES data set. Some of the ML models exhibited slightly lower performance than the DL models (Figure 3). Logistic regression exhibited the highest AUC values for all four data sets. The LGBM classifier and MLP exhibited the second-highest AUC values for the NHANES femoral neck BMD (Figure 3, Table 2). However, two ML models, namely decision trees and KNN, exhibited significantly lower performance than the other ML or DL models. Comparisons with conventional clinical risk assessment tools such as OST, ORAI, and OSIRIS were also performed (Figure 4). OSIRIS exhibited the highest AUC among the three clinical tools (0.771). The AUC of OST was 0.759 and that of ORAI was 0.704 for NHANES, indicating that the developed DL models outperformed the conventional osteoporosis assessment tools.

The contribution rankings calculated using LIME for NHANES and KNHANES data are presented in Table 3 and Table 4, respectively. The results using Boruta and LASSO are presented in Multimedia Appendices 1-4. Based on the feature contribution rankings for NHANES and KNHANES, the high-ranked features were similar between the DL and ML models. These features include respondent sex, age, BMI, arm circumference, and prevalence of obesity. Other features also ranked in high positions, including economic status features such as occupation and family poverty income ratio (PIR), alkaline phosphatase, uric acid, depression, arthritis, and several nutrients (vitamin A, carotene, vitamin C, and vitamin D).

Figure 5 presents a visualization of the impact of the top 10 features according to the integrated contribution rankings of the femoral neck and total femur BMD. Each dot represents an individual respondent. Each feature’s values or categories were normalized and are expressed using different colors. For example, the smaller the BMI value, the more it contributes to the diagnosis of osteoporosis. The individual features and their relationships with the diagnosis of osteoporosis with nonnormalized values are also presented in Figure 5. For example, the graphs of BMI, arm circumference, vitamin D, and family PIR indicate that small values of these features contribute to osteoporosis diagnosis. In contrast, large values of alkaline phosphatase and parathyroid hormone increase the probability of diagnosis in the osteoporosis group.

Figure 6 presents the feature contributions to individual samples of respondents diagnosed with osteoporosis. The impact of each feature value on classifying these respondents into the osteoporosis group is represented in the graph. When comparing these three respondents, female sex increases the probability of being classified into the osteoporosis group. Additionally, the lower the BMI, the more likely an individual is to be diagnosed with osteoporosis. The explainable DL model for these individuals indicates that respondents with osteoporosis tend to have a low BMI and small body size. In terms of prophylactic use of risk factors, analysis results in younger respondents are most appropriate. Therefore, we present the analysis results for the age group of 50-60 years in Multimedia Appendices 5-8.

In this study, the XAI technique was applied to a DL model for osteoporosis risk screening using nationally representative samples from the NHANES and KNHANES data sets. Our DL model outperformed existing ML models, regression analysis, and clinical tools. Subsequently, we performed feature analysis by applying LIME to our DL model. The feature analysis results from LIME reconfirmed previously known features associated with osteoporosis risk, such as sex, age, and factors related to weight (eg, BMI, diagnosis of overweight and obesity). However, unexpected features that are easily ignored in clinical practice were also identified. Furthermore, our model provides individualized risk assessments for osteoporosis with an explanation of feature contributions.

To the best of our knowledge, this is the first study to use a DL model on medical big data to classify osteoporosis and LIME on the DL classifier and to select the most influential features. Previous studies have presented an application of DL or ML models to classify osteoporosis, estimate BMD, and predict osteoporotic fractures. However, previous studies utilizing DL have mainly focused on analyzing image data such as computed tomography [31-34] or X-ray images [35-39]. Most researchers have attempted to use ML models when analyzing tabular data. Additionally, some previous studies that examined important features related to osteoporosis with tabular data attempted to perform statistical analysis or ML based on statistical methods, including a Poisson regression model [40] and multivariate logistic regression [41] for feature selection. Related ML methods include correlation-based feature selection [42], backward elimination [43,44], LASSO applied to logistic regression models [45], use of Gini impurity on a gradient-boosting ML model [46], and outputs of interpretable ML models [47]. Only a few studies have adopted XAI techniques such as the Shapley additive explanation (SHAP), which is confined to ML classifiers. For example, Shim et al [44] used backward elimination to filter out features by considering their contribution to the outputs of a logistic regression model, and Chen et al [47] developed a hybrid model consisting of extreme gradient boosting (XGBoost) and an MLP, where the output values from XGBoost were used for feature selection. Another study that applied SHAP to select the most important features related to osteoporosis was conducted based on XGBoost, which is an ML model [48]. Similarly, Tanphiriyakun et al [49] trained seven ML models to predict the BMD response after the treatment of osteoporosis and analyzed the variable contribution using SHAP. Another study applied feature selection for a minority class method to ML models, KNN and random forest, to select useful information from the balance data [50].

Apart from these previous studies, we developed a DL classifier for screening osteoporosis using cross-sectional data and arranged features in order of importance for classifying respondents into osteoporosis groups using tabular LIME, which is an XAI technique. The strong performance of our model, as evaluated using AUC values, suggests that the DL classifier is also applicable to cross-sectional data analysis and osteoporosis diagnosis. Furthermore, concordance between the features selected through LIME and existing clinical risk factors indicated that the DL model utilized clinically significant features to classify respondents with osteoporosis. Despite the common concern that DL models are black-box models whose complexity makes it challenging to explain the decision-making process, our study demonstrates that interpretability can be achieved even for nonimage data if appropriate explanation techniques are adopted.

Existing studies support the major contributing features identified by LIME and their relationship to osteoporosis. Sex and age are widely known risk factors for osteoporosis, which coincides with several studies suggesting that female sex and low body weight are key risk factors for osteoporosis [51]. Reports on sex discrepancy can also be found for NHANES. For example, in the NHANES 2005-2006 data set, men had a lower prevalence of 30% in the osteopenia group and 2% in the osteoporosis group for femoral neck BMD compared to women with a prevalence of 49% in the osteopenia group and 10% in the osteoporosis group [52]. Other features worth noting for the NHANES and KNHANES data sets are economic status features, including occupation, family PIR, and house ownership. Some studies have observed that postmenopausal women in poverty or non-Hispanic white, Black, and Asian adults with low socioeconomic status have a lower BMD and higher risk of osteoporosis [53,54]. Measurements of alkaline phosphatase and uric acid were also included among the high-ranking features. It is known that patients with osteoporosis exhibit a high level of alkaline phosphatase [55] and patients with high uric acid levels tend to have lower BMD levels [56].

Generally, obesity has been considered to be a protective feature against osteoporosis by providing mechanical stimuli to the bones [57]. Our model indicates that several obesity-associated variables, including BMI, low-density lipoprotein cholesterol, triglycerides, arm circumference, prevalence of obesity, and dyslipidemia diagnosis, are significant features for osteoporosis and exhibit unified trends. In particular, for NHANES, arm circumference emerged as a high-ranking contributing feature. Research on the population living in the southern region of Stockholm indicated that among several anthropometric measurements, a small arm circumference was a stronger risk factor for osteoporosis than a low body weight [58]. Although recent mechanical, biochemical, and hormonal evidence has indicated an association between adipose tissue and deteriorating effects on the microarchitecture of the bone [59], our model identified linear trends for obesity and BMD.

Psychosocial behavior is often ignored as a risk factor for osteoporosis. However, our DL model indicates that depression-associated variables are significant risk factors. This may be associated with decreased physical activity, which is essential for maintaining bone strength, and negative effects of antidepressant medications on bone metabolism, including selective serotonin reuptake inhibitors [60]. Both arthritis and osteoporosis exhibit similar clinical features in terms of symptoms, treatment, and prognosis. However, they are different diseases with different etiologies. A potential role of systemic inflammatory diseases such as rheumatoid arthritis and ankylosing spondylitis in the development of osteoporosis has been proposed [61]. Additionally, several nutrients, including carotene and vitamins A, C, and D, ranked as significant features. Excluding vitamin D, the other nutrients have mainly been considered to be associated with deficient states. However, our model indicates the clinical significance of vitamin D. These diseases and nutrients should not be overlooked, and further research and consideration in clinical practice are required.

There are several limitations of this study that need to be addressed. First, the NHANES and KNHANES data sets utilized in this study were derived from cross-sectional surveys that collected data from a population at one specific point in time. Therefore, this study was confined to predicting the risk of osteoporosis at a specific point in time and we could not generate predictions for future occurrences. The DL model trained on cross-sectional data sets is limited to estimating the current osteoporosis status unless data sets with the label of future osteoporosis status, such as longitudinal cohort data, are utilized. Therefore, we plan to investigate training the DL model with longitudinal cohort data to predict the future risk of osteoporosis in follow-up work. Second, analysis of the BMD of the lumbar spine was not performed because the reference BMD for calculating the T-score for NHANES III did not contain lumbar spine BMD data. Moreover, both data sets lack information about the bone itself, such as bone turnover markers, with only BMD-related information provided. Given the clinical significance of osteoporosis, further use of data on the bone itself, including bone turnover markers, is expected to contribute to improving the model performance and reliability. Third, there was a high level of class imbalance in both the NHANES and KNHANES data sets, which occurred from a small number of respondents being diagnosed with osteoporosis. This imbalance may have caused the DL model to become biased toward the normal group with the highest number of respondents. We attempted applying techniques such as class weighting and the synthetic minority oversampling technique [62] to alleviate the class imbalance problem, but discontinued their use because they had a negative impact on model performance. Additionally, evaluation using the FRAX, which is one of the most widely used osteoporosis clinical assessment tools, was not performed in this study owing to significant differences between the formulae used in the FRAX for the United States and Korea, and the unavailability of an exact formula. Further, there exists a limitation of statistical analysis on the multicategorical features. Our DL model was trained using the original data format for most of the multicategorical features to minimize artificial factors unless unavoidable. However, some multicategorical features showed different results when subgroup analysis was used due to the nature of DL and LIME along with redundancy in features [63-65]. Finally, LIME is prone to instability because of randomness in the sampling step of its process [66]. Although we compared the selected features from LIME to commonly known clinical risk factors for osteoporosis and presented supporting references from previous research, it is still necessary to validate this process. Therefore, the analysis of the selected features from LIME should be performed carefully.

In conclusion, we developed a DL model for classifying osteoporosis using NHANES and KNHANES data, which outperformed conventional clinical assessment tools and ML models. Additionally, we discussed important features selected based on XAI technology. This implies that DL can be fully applied to medical big data for the risk analysis of certain diseases. Furthermore, our model enables individualized osteoporosis risk assessment with an explanation for each feature’s contribution to model results. Despite the limitations noted above, the performance of the DL classifier and feature analysis results are noteworthy and demonstrate the potential for applying DL and XAI techniques to medical research, including cross-sectional studies. Utilizing cohort data for external validation of our DL model and training with the cohort data to predict the future risk of osteoporosis are considered as future work.