Development and Validation of Machine Learning Models for Prediction of Fracture Risk in Patients with Elderly-Onset Rheumatoid Arthritis

Introduction

Rheumatoid arthritis (RA) is a chronic inflammatory joint disease, which can lead to cartilage and bone damage and disability.^1,2 In addition, local and systemic bone loss is one of the major extraarticular complications of RA and leads to an increased risk of brittle fracture, which further impairs functional ability, quality of life, and life expectancy.^3,4 Alarmingly, the incidence rate of osteoporosis in elderly-onset rheumatoid arthritis (EORA) patients has an upward trend in clinical practice.^5,6 The secondary fracture and osteoporosis-related pain in RA patients with osteoporosis can directly lead to the reduction of the quality of life of patients. Therefore, timely fracture risk management is urgently needed for RA patients with osteoporosis.

In the elderly, the presence of complications and the increase of drug-related adverse reactions pose specific treatment challenges.⁶ Although glucocorticoids can reduce clinical activity in the short term and may also reduce structural damage in the medium term, these benefits are offset by many adverse reactions (such as osteoporosis, fracture, etc.).⁷ Therefore, it is of positive significance for clinical treatment to clarify the fracture risk and related influencing factors of RA patients with osteoporosis. At least, to support the balance between over treatment and under treatment, it is necessary to study and propose many strategies to select patients who may omit fracture.

Currently, most risk stratification systems and prediction models are based on classical clinical parameters. For example, a genome-wide polygenic risk score (PRS) was more strongly associated with the risk of fracture, which could define patients with high-risk fracture.⁸ Susan et al established FRAX(^®) calculations to predict the risk of fracture in patients with RA.⁹ Samaneh et al developed the trabecular bone score (TBS) that predicts fracture risk better than bone mineral density (BMD). However, these risk stratification systems and prediction models do not seem to be suitable for accurately predicting fractures, especially in elderly RA patients with osteoporosis. It is worth mentioning that the factors related to bone health play an important role in preventing osteoporosis and its adverse consequences. For instance, a meta-analysis study shows that selenium(Se) deficiency may be one of the risk factors for osteoporosis.¹⁰ Given this situation, adding predictive indexes based on existing models to increase prediction performance and discrimination is crucial.

Machine learning (ML) is increasingly recognized as a useful tool in healthcare applications.¹¹ Nowadays, ML is widely used in the medical field, especially in cancer research, which shows great potential in all aspects, including the benchmark of cancer related issues, such as classification and prediction of cancer types, drug response and treatment strategies.^12,13 In addition, the use of ML can also help improve the reliability, performance and accuracy of specific disease diagnosis systems.^14,15 Compared with the traditional statistical model, ML-based integrated analysis is famous for ensuring the robustness of the model and improving the prediction accuracy through repeated iteration algorithms. In this study, we tried to apply ML-assisted decision-support models, using classical clinical parameters, to calculate the probability of fracture. This work may help to optimize the precise treatment and further improve the clinical efficacy of elderly RA patients with osteoporosis.

Materials and Methods

Study Population

We retrospectively analyzed the clinical records of 487 elderly patients with RA and osteoporosis who received blood biochemical measurements at the Central Hospital of Enshi Tujia and Miao Autonomous Prefecture of China from January 2015 to December 2021. The inclusion criteria were as follows: (i) RA diagnosis was clinically diagnosed by a committee-certified rheumatologist; (ii) Patient age ≥ 60 years old; (iii) Patient was accompanied with typical RA symptoms such as joint swelling, pain, deformity, and limited function. To minimize the confounding effects of circulating biomarkers, the exclusion criteria were as follows: (i)Patients with other autoimmune diseases (such as systemic lupus erythematosus, scleroderma, mixed connective tissue disease, etc.); (ii)Patients who continued to take drugs such as bisphosphate and estrogen within half a year; (iii) Patients with severe organic diseases (such as heart, liver, kidney, cerebrovascular diseases, etc.); (iv) The peripheral blood picture of the patient decreased significantly, such as the white blood cell count ≤ 3.0 × 10⁹/L, platelet count ≤ 8.0 × 10⁹/L, etc; (v) Patient was diagnosed with a malignant tumor. This study complies with the provisions of the Helsinki Declaration (revised in 2013) and was approved by the Institutional Review Committee of The Central Hospital of Enshi Tujia and Miao Autonomous Prefecture. The workflow for RA patient selection and model construction was summarized in Figure 1.

Figure 1 The flow chart of patient selection and data process. Abbreviations: RFC, random forest classifier; SVM, support vector machine; DT, decision tree; ANN, artificial neural network; XGboost, eXtreme gradient boosting; PMOF, probability of major osteoporotic fractures; PHF, probability of hip fracture.

Data Collection and Quality Assessment

We extracted eighteen features including epidemiological, demographic, clinical, laboratory, radiological, and outcome data from an electronic case recording system, including gender composition, age, body mass index (BMI), waist circumference, RA course, fracture history, hip fracture of parents, smoking, long-term glucocorticoid treatment history (previous oral glucocorticoids for more than 3 months), daily alcohol intake ≥ 3 units, calcium intake and 28 joint disease activity (DAS28), calcium intake was judged according to the patient’s diet (for example, the daily intake of elemental calcium for adults is 800mg). For variables with missing values, the median was typically used. If ≥10% of values were missing for a given variable, it was excluded from variable screening for the final model. In addition, all data entries were entered and proofread by trained professionals. To reduce the risk of bias in the included data, all serological indicators were subject to the data collected for the first time after admission. In addition, the definition of RA disease activity: DAS28 score ≤ 2.6 is remission, 2.6–3.2 is low activity, 3.2–5.1 is moderate activity, and >5.1 is high activity.

Measurement of Bone Mineral Density(BMD)

A dual-energy X-ray bone densitometer (purchased from GE company of the United States, model: lu43616cn) was used for anteroposterior photography to measure the BMD of the patient’s lumbar spine, hip, and femoral neck. The T value was calculated according to the BMD reference value of healthy adults of the same sex and nationality.¹⁶ T value ≥ −1.0 indicates normal bone mass; −2.5<t value <-1.0 indicates bone loss; T value ≤ −2.5 indicates osteoporosis.

Fracture Risk Assessment

We input relevant risk indicators and BMD values from the FraX interface to calculate the 10-year probability of major osteoporotic fractures (PMOF) and hip fractures (PHF).¹⁷ Fracture risk: PMOF>20% indicates high risk, 10%~20% indicates medium risk, and <10% indicates low risk; PHF ≥ 3% indicates medium and high risk, and <3% indicates low risk.

Development and Validation of ML-Based Models

The data were randomly divided into a training set (70%) and a verification set (30%) to verify the prediction model. We trained the model with random variables and the results of 70% of patients. During model training, we installed ten-fold cross-validation to fine-tune model parameters without including baseline ML-based predictive models into the training set queue, that is, random forest classifier (RFC) model, artificial neural network (ANN), support vector machine(SVM), eXtreme gradient boosting (XGBoost), decision tree (DT).¹⁸ The increased weight of minority categories in the model can increase the identification of classification errors during training and improve the ability of the model to identify minority categories. The model variables (ie characteristic variables) were screened according to the principle of “OOB error”,¹⁹ as follows:

If the Gini index is smaller, the probability that the selected samples in the set will be mixed will be small, that is, the higher the purity of the set is, on the contrary, the more impure the set is. However, if all the samples in the set are of the same class, the Gini index approaches zero.²⁰

Prediction Efficiency Evaluation of ML-Based Models

The optimal subset variables of the modeling were obtained based on the intersection of variable sets. The predictive performance of the model was evaluated by the receiver operating characteristic (ROC) curve, and evaluation indexes including area under the ROC curve (AUC), accuracy, sensitivity, specificity, positive predictive value (PPV), and negative predictive value (NPV). The discrimination ability of each model was quantified by the decision curve analysis (DCA) and clinical impact curve (CIC).^21,22 That is, we call the ggDCA program, analyze all the included variables, and obtain the curve of the best prediction model through gradual screening. The farther away from the “threshold curve”, the best performance of the model will be considered.

Statistical Analysis

For descriptive analysis, median (IQR) and frequencies (%) were assessed for continuous and categorical variables, respectively. Bonferroni corrected probability values are used to compare the qualitative data.²³ Wilcoxon rank-sum test or chi-square test was used to compare the differences between different groups. For parameters whose variables do not conform to normal distribution, Mann Whitney U-test was used for comparison between groups. All analysis was performed using the R Project for Statistical Computing (version 4.0.4, http://www.r-project.org/). All P values were two-tailed, and P <0.05 was considered statistically significant.

Results

Baseline Characteristics of the Study Population

The baseline characteristics of 487 hospitalized RA patients with osteoporosis were summarized in Table 1. For internal validation, patients were randomly divided into a training set (N=340, 70%) and a validation set (N=147, 30%) using the caretpackage. There were significant differences in age, BMI, RA course, long-term glucocorticoid treatment history, DAS28 score, ESR, RF, 25 hydroxyvitamin D3, PMOF, and PHF between the osteoporosis group and non-osteoporosis group (P<0.05). However, there were no significant differences in gender, calcium intake, waist circumference, fracture history, CRP, and anti-CCP antibody (P>0.05).

Table 1 Characteristics of Baseline Demographic and Clinical Test Indicators of Patients

Comparison of Fracture Risk Stratification Between Two Groups

According to PMOF combined with PHF, 487 elderly RA patients with osteoporosis were classified as medium and high risk; Among 116 non-osteoporosis patients, 38(32.8%) had normal bone mass without fracture risk, and 68 had decreased bone mass, including 53(45.7%) low-risk patients and 25(21.6%) medium and high-risk patients. As shown in Supplementary Table 1, there was a significant difference in fracture risk between the osteoporosis group and the non-osteoporosis group (P<0.05). To further evaluate the fracture prediction efficiency of PMOF and PHF, the area under the ROC curve(AUC) showed that the fracture prediction efficiency of PMOF and PHF in the training set was 0.725(95% CI: 0.674–0.776) and 0.695(95% CI: 0.644–0.746), respectively (Table 2).

Table 2 The ROC Curve Analyses for Predicting Fracture Risk in Each ML-Based Model

Selection of Candidate Variables

Candidate covariates of each algorithm were filtered and twenty-two were included in the correlation analysis between outcome and independent variables. The correlation matrix revealed that fracture significantly correlated with image factors and some clinical variables (Figure 2A). Additionally, each meaningful candidate variable, including PMOF, PHF, serum selenium, and clinical indices, that is, glucocorticoid, and erythrocyte sedimentation rate(ESR) were important factors in the ML-based model (Figure 2B).

Figure 2 Variable screening and weight allocation. (A) Correlation matrix analysis of candidate features. (B) The weight distribution of the candidate variables of each ML-based model. Abbreviations: RFC, random forest classifier; SVM, support vector machine; DT, decision tree; ANN, artificial neural network; XGboost, eXtreme gradient boosting.

Construction of ML-Based Fracture Predictive Model

Random forest classifier (RFC) and decision tree (DT) are commonly used ML-based algorithms in supervised learning. In this study, The RFC model was constructed using the formula I(X=xi)=−log2P(xi), where I(X) is the information for candidate variables and P(xi) is the probability of xi (Figure 3A). Twenty-two variables were ordered according to the mean decrease in the Gini index (Supplementary Table 2); the top seven ranked variables were used to construct the optimal RFC prediction model, which included PMOF, PHF, BMD, RF, and Se. Similarly, the Se also served as irreplaceable weight at DT branches (Figure 3B). Using the iterative algorithm of supervised learning, both RFC and DT models were used for fracture prediction. Meanwhile, the ANN model showed more robust prediction efficiency than other models but was inferior to the RFC (Figure 4).

Figure 3 Predictive model visualization based on ML-based algorithm. (A) RFC model. (B) DT model. Notes: The candidate factors associated with fracture risk were ordered via RFC algorithm (A and B) prediction node and weight were allocated via DT algorithm.

Figure 4 Predictive model visualization based on ANN algorithm. (A) ANN model. (B) Variable importance using connection weight. Notes: The candidate factors associated with fracture risk were ordered via ANN algorithm (A and B) prediction node and weight were allocated via ANN algorithm.

Comparison Among ML-Based Models

Based on the iterative analysis of baseline characteristics, we used five ML-based supervised learning models and two common predictive tools (PMOF, PHF) for fracture risk assessment and to optimize predictive performance. As expected, the RFC model was better able to predict fracture risk. The AUCs of the RFC model reached a plateau when seven variables were introduced, indicating that the RFC model had the highest predictive accuracy, followed by DT, artificial neural network (ANN), support vector machine (SVM), and eXtreme gradient boosting (XGBoost) models (Figure 5). The predictive performance of ML-based models was summarized in Table 2 and Supplementary Table 3. Interestingly, the prediction performance of all supervised learning prediction models was significantly better than that of PMOF and PHF prediction tools.

Figure 5 Prediction performance of candidate models based on ML-based algorithm. (A) DCA for five ML-based models in the training set. (B) DCA for five ML-based models in the testing set. Abbreviations: RFC, random forest classifier; SVM, support vector machine; DT, decision tree; ANN, artificial neural network; XGboost, eXtreme gradient boosting.

Internal Validation of the Optimal Predictive Model

To further verify the optimal prediction model, CIC was used to verify the RFC model with the training set and test set respectively. This was supported by fracture risk factors identified in the training focus, as shown in Supplementary Figure 1. In addition, these were also consistent with the results of validation cohorts, indicating that RFC had the best performance across the metrics of discrimination, calibration, and overall performance, especially the candidate clinical indices that were highly relevant to fracture.

Discussion

In this study, we performed a longitudinal cohort study to examine the risk of incident fracture among RA patients with osteoporosis. To our knowledge, few previous studies have assessed the risk of fracture in EORA. Therefore, this study detected the fracture incidence of elderly RA patients with osteoporosis for the first time, we successfully constructed an early warning classification model for guiding diagnosis and treatment.

In this study, we found RA patients with osteoporosis had an increased prevalence of risk factors for fracture. The advantages of this study included a large cohort of patients (487 cases) with an average follow-up of 68.5 years, the use of fracture prediction, in which the definition of exposure (glucocorticoid, RA, BMD code) has been validated, the examination before fracture, and the ability to adjust for other measured risk factors of osteoporosis, including combination medication, BMI and ESR.^24–27In addition, since our estimates of rheumatoid arthritis are similar to previous studies, the inclusion of a group of patients with rheumatoid arthritis for internal comparison can provide validity to the results.^28,29

Previous studies have investigated the effects of Se on osteoporosis, however, the results are inconclusive.¹⁰ Notably, in this study, we found that Se expression was associated with a higher risk of fracture in RA patients with osteoporosis because the expression was significantly lower in patients with fractures than in patients without fractures. Previous studies have shown that selenoprotein expressed in human fetal osteoblasts seems to protect bones from oxidative stress, which may lead to osteoporosis by inhibiting the osteogenic differentiation of bone marrow stromal cells.^30,31 We speculate that the trace element selenium is an important component of selenoprotein and may play an important role in the relationship between Se and bone mineral density. Consistent with the expected results, Se accounts for a large weight in the multiple ML-based prediction models, which indicated that Se has potential value in predicting fracture risk.^32,33

Supervised machine learning algorithm has always been the mainstream method in the field of data mining.³⁴ Recently, the use of health data for disease prediction has shown the potential application of these methods.^35,36 This study demonstrated that the use of machine learning models can accurately predict fracture in RA patients with osteoporosis. As an important branch of supervised learning, the RFC model has been successfully applied to high-dimensional and multi-source data reduction of many diseases.³⁷ The popular supervised learning classifiers, including support vector machines, random forest, convolutional neural networks, and decision trees, have been gradually applied in clinical practice. Consistent with the results of previous research reports,^38,39 we found that the RFC model has more advantages in feature selection and classification, and shows better prediction performance than PMOF and PHF.

Our study also has limitations. First, This is a single-center study that needs to be validated with data from other sources. The variables that are input to machine learning algorithms are usually variables that can be obtained or evaluated in most cases, therefore, it is necessary to conduct repeated validation using data from other more clinical medical centers. Second, the added value of roles of serum selenium in RA patients with osteoporosis remains unknown, and further in vivo and in vitro experiments are needed to reveal their functions. Third, all of the samples from this study is retrospective, and future validation should be performed in a prospective multicentre cohort.

Conclusion

In conclusion, fracture, especially osteoporotic fracture, is a major health problem, leading to adverse consequences. Osteoporosis is largely underestimated in the elderly population. Additionally, we have successfully constructed the RFC with the highest prediction accuracy. Further research is needed to develop these models so that they can be used in daily practice and care for elderly RA patients with osteoporosis.