Sage Journals: Discover world-class research

Abstract

Background:

Immune cells are involved in rheumatoid arthritis (RA), but the link between other blood cell indices and the disease activity of RA, along with the underlying mechanisms, is unclear.

Objective:

This study aimed to develop an interpretable machine learning model based on blood cell parameters to assess RA disease severity and assist in personalized treatment decisions.

Methods:

A retrospective case–control study was conducted with blood routine and biochemical detection data from 4401 patients at the First Affiliated Hospital of Guangxi Medical University, spanning from January 1, 2018, to January 1, 2024. The primary outcome was disease severity stratification. Recursive feature elimination was applied to identify key variables, and 10 machine learning algorithms were benchmarked on 55 clinical features with internal validation. Model interpretability was assessed with SHAP, while logistic regression and restricted cubic spline models were used to examine associations between blood cell indices and disease severity. In addition, Mendelian randomization analysis was performed to explore potential causal relationships.

Design:

This was a retrospective case–control study.

Results:

Blood cell indices were identified as the primary factors associated with RA severity. In model evaluation, the Random Forest achieved the best performance, with test set AUCs of 0.870 and 0.874. Mendelian randomization supported a causal relationship between blood cell indices and RA risk.

Conclusion:

These results reinforce the associations between blood cell indices and RA severity. The machine learning model demonstrates good predictive capabilities for RA severity and may assist clinicians in developing personalized treatment strategies.

Plain language summary

Blood tests predict rheumatoid arthritis severity to customize treatments

What is already known on this topic? Rheumatoid arthritis (RA) is a chronic disease where the immune system attacks joints, leading to pain, swelling, and potential permanent damage. Blood tests track inflammation via proteins like CRP, but they’re not precise for gauging RA severity. Immune cells drive RA, yet links between routine blood cell counts (e.g., red/white cells, platelets) and disease worsening—and their causes—are unclear.

What this study adds? We analyzed records from 4,401 RA patients at a Chinese hospital (2018–2024). We built a computer tool from blood tests to classify RA as mild or severe, aiding tailored treatments. Using 10 algorithms on 55 features, we selected top blood indicators (e.g., lymphocytes, platelets), explained predictions with SHAP, and probed causality via Mendelian randomization and nonlinear models.

What was found? Blood cells were top predictors of severity. Random Forest excelled (87% accuracy). Low lymphocytes and high platelet variation boosted severe RA odds (ORs 0.54–2.17). Causal ties confirmed: these changes drive RA risk. Nonlinear patterns showed extremes (too low/high counts) worsen outcomes.

What the results mean for patients and the public? This tool uses affordable blood tests for fast severity checks, enabling personalized meds to curb damage and pain. Patients get targeted care, fewer trial errors, better daily life. It promotes equitable RA management in resource-limited settings and inspires blood-focused prevention strategies.

Keywords

blood cell indices machine learning Mendelian randomization predictive model rheumatoid arthritis

Introduction

Rheumatoid arthritis (RA) is a chronic inflammatory autoimmune disease that frequently progresses to irreversible joint damage and significant physical disability.¹ RA primarily targets synovial joints, leading to progressive bone erosion and disability through chronic synovitis.² Blood biomarkers increase dramatically in response to tissue damage and inflammatory processes. Therefore, characterizing blood biomarker responses in RA critical to gain insight into the activity of the disease.^3,4 Most existing studies focus on the association between acute reactants and disease activity, but lack studies on markers associated with disease severity.

Blood biomarkers, such as C-reactive protein (CRP) and erythrocyte sedimentation rate (ESR), have long been used to monitor RA activity.^5,6 However, these markers lack specificity and may not fully capture the complex interplay between hematological and immune systems in RA progression. Emerging evidence suggests that routine blood cell indices including platelet count, neutrophil-to-lymphocyte ratio (NLR), and red cell distribution width (RDW) reflect systemic inflammation and oxidative stress, processes intimately linked to RA pathogenesis.^7–10 For instance, platelets contribute to synovial inflammation via cytokine release, while elevated NLR correlates with pro-inflammatory cytokine levels.^11,12 Despite these insights, the relationship between comprehensive blood cell indices and RA severity, the causal relationship between blood cell indices and RA, and the diagnostic value of blood cell indices in predicting RA severity remain underexplored.

In this study, the machine learning (ML) models built using blood cell indices can more accurately predict the severity of the disease. By analyzing a large retrospective case–control of incident RA patients, we establish associations between blood cell indices and RA severity. Utilizing interpretable ML models to predict disease severity stratification and identify key hematological thresholds with clinical relevance to support personalized management of patients with rheumatoid arthritis. Additionally, we employed Mendelian randomization (MR) analysis to assess causal associations.

Methods

Data sources

Patients diagnosed with RA were retrospectively examined for inclusion in the study. The reporting of this study conforms to the Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) statement.^13,14 All patients in this study demonstrated unequivocal radiographic evidence of joint erosions, accompanied by pain or swelling in the hands, wrists, or other small joints, and met the 2010 American College of Rheumatology/European League Against Rheumatism (ACR/EULAR) classification criteria for RA.¹⁵ A comprehensive dataset, encompassing demographic characteristics, hematological indices, biochemical indices, and derived hematological parameters, was extracted from electronic medical records. The dataset included demographic information (age and sex), hematological indices (red blood cell count (RBC), hemoglobin concentration (HGB), hematocrit percentage (HCT), mean corpuscular volume (MCV), mean corpuscular hemoglobin (MCH), mean corpuscular hemoglobin concentration (MCHC), red blood cell distribution width (RDW), white blood cell count (WBC), basophil count (BASO), basophil percentage (BASO%), eosinophil count (EO), eosinophil percentage (EO%), monocyte count (MONO), monocyte percentage (MONO%), neutrophil count (NEUT), neutrophil percentage (NEUT%), lymphocyte count (LYMPH), lymphocyte percentage (LYMPH%), platelet crit (PCT), mean platelet volume (MPV), platelet distribution width (PDW), and platelet count (PLT)), biochemical indices (total cholesterol (TC), direct bilirubin (DBIL), indirect bilirubin (IBIL), direct-to-total bilirubin ratio (DTR), total protein (TP), albumin (ALB), globulin (GLOB), albumin-to-globulin ratio (AGR), gamma-glutamyl transferase (GGT), total bile acids (TBA), aspartate aminotransferase (AST), alanine aminotransferase (ALT), aspartate-to-alanine aminotransferase ratio (AAR), alkaline phosphatase (ALP), prealbumin (PA), urea (UREA), creatinine (Cr), uric acid (UA), bicarbonate (HCO₃^-), endogenous creatinine clearance (Ccr), cystatin C (CysC), high-sensitivity troponin (hsTn), creatine kinase (CK), creatine kinase isoenzyme (CK-MB), lactate dehydrogenase (LDH), alpha-hydroxybutyrate dehydrogenase (HBDH), random blood glucose (RBG), and erythrocyte sedimentation rate (ESR)), and derived hematological parameters (neutrophil-to-lymphocyte ratio (NLR), monocyte-to-lymphocyte ratio (MLR), and platelet-to-lymphocyte ratio (PLR)). Each indicator’s corresponding abbreviation is provided in Table S1. Participants with other inflammatory or infectious diseases (N = 109), autoimmune disorders (N = 13), lacking follow-up data (N = 461), or presenting serious complications or tumors (N = 76) were excluded from the study. Ultimately, the primary analysis included 4401 participants (Figure 1). Prior to analysis, data quality control was performed. Continuous variables were examined for implausible outliers, which were handled using winsorization at the 1st and 99th percentiles. Missing values were rare (<5% for all variables), and all samples with any missing values were excluded from the analysis. As a sensitivity check, analyses were repeated using multiple imputation by chained equations, and the results were consistent with the primary analysis. Patients were eligible if they: (1) were aged 18 years or older; (2) had a confirmed diagnosis of RA according to the 2010 ACR/EULAR classification criteria; (3) had complete baseline laboratory and clinical data, including Disease Activity Score 28 based on ESR (DAS28-ESR) and the Clinical Disease Activity Index (CDAI). Exclusion criteria were: (1) concomitant autoimmune diseases other than RA; (2) active infection, hematologic disorders, or malignancy at baseline; (3) missing key clinical or laboratory information; (4) pregnancy or lactation.

Figure 1.

Flowchart of the study design.

Definition of disease severity

The RA patients’ treatment plan strictly followed the management of rheumatoid arthritis: summary of NICE guidance.¹⁶ Disease severity was primarily determined using validated and widely accepted disease activity indices, including the DAS28-ESR and CDAI. Patients with DAS28 ⩽3.2 and CDAI ⩽10 were classified as having mild disease, whereas those with DAS28 >5.1 and CDAI >22 were classified as having severe disease. To better contextualize the clinical characteristics of each group, we reviewed accompanying clinical information (e.g., care setting or treatment intensity). However, these variables were not used to define severity labels; instead, they served only as supportive descriptors to confirm overall clinical plausibility. This approach ensured that severity labeling was based solely on standardized disease activity indices, thereby minimizing circularity and reducing the risk of treatment-related confounding.

Development and performance benchmarking of ML models

All features, including demographic characteristics, hematological indices, biochemical indices, and derived hematological parameters, were included in the development of 10 machine learning models. Compare and evaluate the performance of 10 machine learning methods through benchmarks. Benchmarking is a common practice in machine learning to assess the performance of different algorithms. The “mlr3verse” R package allows for comparison of multiple learners (algorithms), tasks (data), and average performance across various resampling strategies (multiple data copies). Ensuring the fairness of tests is essential in benchmarking. Each algorithm must train on the same resampling training set and be evaluated on the same resampling test set. These tasks are automatically handled by the benchmark() function with a repeated 10-fold cross-validation approach designed to maximize performance metrics.^17–19 The optimal model was identified by benchmarking its performance on both the training and test sets using a comprehensive set of evaluation metrics. These included accuracy, sensitivity, specificity, precision, F1 score, and Brier score, all calculated based on the test set results with model-specific decision thresholds. To further assess model reliability, calibration curves were plotted to evaluate the agreement between predicted probabilities and observed outcomes. Additionally, decision curve analysis was performed to examine the clinical utility of each model across different threshold probabilities. The use of the independent test set ensured an objective and robust comparison, facilitating the selection of the best-performing model.

Feature selection and optimized ML model development

To ensure appropriate scaling for algorithms sensitive to feature magnitude, all continuous variables were standardized prior to model training. Specifically, z-score transformation (subtracting the mean and dividing by the standard deviation) was applied using parameters derived from the training set, and the same parameters were used to transform the validation and test sets. Standardized features were used for logistic regression, Support Vector Machine (SVM), neural network, K-Nearest Neighbors (KNN), and other distance- or gradient-based algorithms. Although tree-based models (Random Forest (RF), Extreme Gradient Boosting (XGB), LightGBM, CatBoost) are generally invariant to monotonic scaling, the same preprocessing pipeline was applied uniformly for consistency and reproducibility.

To identify key indices associated with RA severity, recursive feature elimination (RFE) was employed.²⁰ RFE is a feature selection algorithm that iteratively trains a model, ranks features based on their importance, and eliminates the least important features until an optimal subset is identified. In this study, RFE was implemented using a RF with the best performance classifier to evaluate feature importance, with 10-fold cross-validation to ensure robust performance assessment.²¹ The initial input included all variables. The RFE process was iterated until model performance convergence. The final selected features were used to reconstruct the predictive models. All analyses were performed using “pbapply,” “rlang,” “tidyverse,” “reshape2,” “openxlsx” R packages. To train and validate the models, the dataset was randomly split into 70% for training and 30% for testing. Ten machine learning algorithms were implemented, including Logistic Regression, SVM, Gradient Boosting Machine (GBM), Neural Network, RF, XGBoost, KNN, AdaBoost, LightGBM, and CatBoost. Final hyperparameters for each model were optimized using a combination of grid search and manual fine-tuning.^22,23 To improve reproducibility, we provide a detailed description of the hyperparameter tuning procedure. First, for each model, we defined a predefined hyperparameter search space (e.g., number of trees = 100–1000 and maximum tree depth = 3–15 for RF; learning rate = 0.01–0.3 and max_depth = 3–10 for XGBoost; C = 0.1–10 and kernel = radial basis/linear for SVM). Grid search was performed using a 5-fold cross-validation strategy on the training set to identify the optimal combination of hyperparameters within these ranges. Model performance was evaluated using average area under the curve (AUC) and Brier score across folds. After grid search, a second step of manual fine-tuning was performed only for the top-performing models (RF, XGBoost, LightGBM). This process involved adjusting a small set of sensitive parameters (e.g., learning rate, tree depth, or minimum child weight) within a narrower range around the grid-search-derived values. Adjustments were guided by cross-validated AUC and calibration performance. Because treatment factors such as conventional synthetic disease-modifying antirheumatic drugs (DMARDs), biologic DMARDs, targeted synthetic DMARDs, and surgery may influence both RA severity and hematological indices, we carefully considered treatment confounding. In the primary analysis, disease severity was defined independently of treatment factors, based solely on DAS28-ESR and CDAI. Additionally, we performed a sensitivity model excluding ESR and re-evaluated discrimination and calibration.

Explainable machine learning model

In this study, Shapley Additive Explanations (SHAP) was applied to gain insights into feature importance, interactions, and relationships.^24–26 SHAP supports targeted feature engineering by identifying the most impactful features for ML models.²⁷ It assigns a level of importance to each individual patient (i.e., data point) based on the corresponding Shapley values. The SHAP value of a clinical variable reflects its average contribution across all possible combinations of clinical variables. These values can be positive or negative, representing an increased or decreased likelihood of achieving a clinical diagnosis, respectively.^24,28 Model interpretation reveals not only the importance of features but also the direction of their impact on the predicted outcomes. This insight can help physicians understand potential associations between features and results.^29,30 SHAP analysis can reveal the differences in features between different patients, helping doctors identify the heterogeneity of patient groups and thus perform stratified treatment.²⁹ Two commonly used visualization tools for model interpretation are SHAP waterfall charts and SHAP dependency graphs. SHAP waterfall charts are used to explain the prediction results of a machine learning model for a specific sample. It shows how the model’s prediction changes from the expected value (E[f(X)]) to the final predicted value (f(X)).³¹ SHAP dependency graphs are used to show the impact of a single feature on the model’s prediction. This helps us understand how the feature affects the model’s output by plotting the relationship between the value of the feature and its corresponding SHAP value.³² Interpretable ML methods were implemented via Python packages based on “NumPy,” “pandas,” “scikit-learn,” and “SHAP.” The code used in this study is available from the corresponding author upon reasonable request in accordance with the institutional confidentiality agreement.

Decision threshold determination

To enable clinical use of the model, a probability decision threshold was derived to convert predicted probabilities into binary severity classifications. The optimal threshold was determined on the training set using the Youden index (maximizing sensitivity + specificity − 1). For robustness, thresholds maximizing the F1 score and thresholds achieving a minimum sensitivity of 0.80 were also evaluated. The selected threshold was then applied to the held-out test set, and performance metrics including sensitivity, specificity, positive predictive value (PPV), negative predictive value (NPV), accuracy, and F1 score were calculated. The final threshold and its corresponding diagnostic metrics are provided in Table S2.

Covariates

At baseline, relevant covariates were chosen based on clinical expertise and previous research, including demographic factors such as age and sex, along with various blood cell indices and biochemical markers (see Table S3 for the complete list). To examine multicollinearity among these covariates, the variance inflation factor (VIF) was computed, and all VIF values were <10, signifying the absence of multicollinearity (Table S3 in the supporting information). Baseline laboratory tests were defined as those obtained within ±14 days of the index visit, which corresponded to the time of severity assessment. This time window was chosen to ensure temporal consistency between exposure and outcome measurements.

Statistical analysis

Baseline characteristics were described according to RA status, with continuous variables presented as medians (interquartile ranges (IQRs)) and categorical variables as frequencies (percentages). A flow diagram depicting the study design is shown in Figure 1. The top 15 significant clinical variables were ranked based on the SHAP feature importance plot. The Age, ESR, PLR, NLR, and ALP were excluded. Ten key blood cell indicators were identified: MPV, RBC, MCV, PLT, LYMPH, LYMPH%, HCT, PDW, NEUT%, and PCT. These 10 indicators were considered primary variables for their association with the severity of incident rheumatoid arthritis. Before performing binary logistic regression, the blood cell indices were standardized into Z scores (Z = [value − mean]/SD) to ensure that the odds ratio (OR) represented the effect of each standard deviation increase. The non-linearity was tested using cubic splines with three knots set as 10th, 50th, and 90th percentiles to facilitate comparison against the models with linear terms. Model I, adjusted for demographic variables, served as the baseline model. Model II made additional adjustments for other blood cell indices. In the sensitivity analysis, Model II was applied to examine the robustness of the results. All statistical analyses were carried out using IBM SPSS Statistics, version 27 (IBM Corp., Armonk, NY, USA).

MR analysis of blood cell indices on RA risk

We first retrieved the single nucleotide polymorphisms (SNPs) from the Integrative Epidemiology Unit (IEU) genome-wide association study (GWAS) database, and selected those SNPs that showed a p-value smaller than the GWAS significant threshold (p-value <5e-6). SNPs were then clumped based on the “TwoSampleMR” package in statistical software R, with the default linkage disequilibrium threshold r² <0.001, distance >10,000 kb, and 1000 genomes European reference panel. Palindromic and ambiguous SNPs were further removed after harmonizing the data between exposure and outcome datasets (i.e., the default option in TwoSampleMR Rpackage). SNPs with F-statistic <10 were excluded to avoid weak instrument bias.³³ We employed multiple MR methods as summarized in Table S4.^34–37 In brief, random-effect inverse variance weighted (IVW) provides unbiased estimates when all three MR assumptions are satisfied. The presence of pleiotropic SNPs that influence the outcome through pathways that bypass the exposure violates the exclusion restriction assumption of MR and undermine the validity of the estimates. Weighted median is robust to such outliers provided that the majority of selected genetic variants are valid instruments.³⁶ Fixed-effect inverse variance weighted (FE-IVW) provides unbiased and precise estimates when all three Mendelian randomization assumptions are satisfied, and there is no heterogeneity among the genetic variants. The presence of pleiotropic SNPs that influence the outcome through pathways bypassing the exposure violates the exclusion restriction assumption of MR, leading to biased estimates. FE-IVW is highly sensitive to heterogeneity, and its validity is undermined if the assumption of a common effect size across all variants is violated.³⁸ MR-Egger relaxes the assumption of exclusion restriction in the situation where the strengths of the instruments are uncorrelated with their direct effects on the outcome (the InSIDE assumption).³⁹ The Simple Mode method in Mendelian randomization clusters genetic variants to estimate causal effects, minimizing bias by assuming most variants are valid.⁴⁰ Weighted mode clusters genetic variants to estimate causal effects, allowing pleiotropic effects, and ensures consistency using kernel-smoothed density.⁴⁰ In this study, MR-Egger, weighted median (WM), IVW, simple mode, and weighted mode were used to conduct MR analysis to better evaluate the causal relationship between RA and blood cell count indices. When the level of pleiotropy of instrumental variables (IVs) does not exist, IVW has the highest statistical power; therefore, it is used as the main analysis method.³⁸ Different types of genetic pleiotropy are considered by other MR methods and are based on potentially different assumptions to check the robustness of the results. For example, WM can still provide unbiased causal estimation, when more than 50% of the information comes from invalid IVs.³⁴ The MR-Egger method can also provide an unbiased estimation of causality when all IVs do not meet the core MR hypothesis.³⁹ The effect of exposure on outcome was expressed using an OR and 95% confidence interval (95% CI). All the MR analyses conducted in this study were performed using “TwoSampleMR” packages in the R version 4.4.3.

Results

Basic characteristics of participants

The baseline characteristics of the study cohort, stratified by the severity of RA (mild vs severe), are presented in Table 1. This study included 4401 patients with RA, of whom 28.6% were male (n = 1260) and 71.4% were female (n = 3141). The median age of the overall cohort was 58 years, with the severe RA group being significantly older than the mild RA group. Significant differences were observed between the mild and severe RA groups in multiple hematological and biochemical indices, including RBC, hemoglobin, PLT, and ESR (all p < 0.05), as detailed in Table 1.

Table 1.

Baseline characteristics of patients with rheumatoid arthritis stratified by disease severity.

Characteristics	Overall (n = 4401)	Mild RA (n = 2573)	Severe RA (n = 1828)	p-Value
Sex (%)
Male	1260 (28.6)	717 (27.9)	543 (29.7)	0.184
Female	3141 (71.4)	1856 (72.1)	1285 (70.3)
Age, median (IQR)	58 (50, 68)	56 (46, 67)	61 (54, 68)	<0.001
Blood count indices
WBC, median (IQR)	6.770 (5.285, 8.740)	6.72 (5.33, 8.54)	6.83 (5.23, 9.10)	<0.001
RBC, median (IQR)	3.97 (3.48, 4.45)	4.11 (3.69, 4.59)	3.73 (3.17, 4.18)	<0.001
HGB, median (IQR)	110.5 (96.2, 123.1)	113.0 (101.0, 126.0)	106.0 (91.0, 119.8)	<0.001
HCT, median (IQR)	0.340 (0.297, 0.376)	0.352 (0.316, 0.386)	0.320 (0.264, 0.361)	0.034
MCV, median (IQR)	86.50 (80.30, 91.46)	86.50 (80.71, 91.40)	86.48 (79.55, 91.68)	0.005
MCH, median (IQR)	28.19 (25.09, 30.20)	28.20 (25.10, 30.16)	28.16 (24.93, 30.20)	0.003
MCHC, median (IQR)	324.0 (314.0, 332.0)	324.0 (313.0, 332.4)	324.0 (314.9, 331.3)	<0.001
PLT, median (IQR)	282.7 (224.3, 352.9)	263.4 (215.3, 324.0)	320.2 (251.0, 408.4)	0.023
MPV, median (IQR)	8.37 (7.39, 9.40)	8.60 (7.60, 9.60)	7.92 (7.12, 9.10)	0.176
PDW, median (IQR)	16.0 (11.0, 17.0)	16.0 (10.0, 16.0)	16.0 (16.0, 17.0)	0.059
NEUT, median (IQR)	4.25 (3.00, 6.05)	4.12 (2.95, 5.72)	4.52 (3.08, 6.51)	0.002
NEUT%, median (IQR)	0.640 (0.549, 0.727)	0.624 (0.532, 0.709)	0.662 (0.578, 0.758)	<0.001
LYMPH, median (IQR)	1.53 (1.16, 2.02)	1.67 (1.28, 2.15)	1.37 (0.99, 1.83)	0.395
LYMPH%, median (IQR)	0.227 (0.159, 0.308)	0.254 (0.191, 0.335)	0.185 (0.133, 0.260)	0.557
MONO, median (IQR)	0.53 (0.39, 0.71)	0.52 (0.40, 0.69)	0.54 (0.38, 0.74)	<0.001
MONOpct, median (IQR)	0.078 (0.061, 0.100)	0.078 (0.062, 0.100)	0.078 (0.061, 0.100)	<0.001
EO, median (IQR)	0.120 (0.065, 0.220)	0.130 (0.070, 0.210)	0.120 (0.057, 0.230)	<0.001
EO%, median (IQR)	0.019 (0.009, 0.034)	0.019 (0.009, 0.033)	0.019 (0.008, 0.034)	<0.001
BASO, median (IQR)	0.03 (0.02 0.04)	0.03(0.02, 0.04)	0.03 (0.02, 0.04)	<0.001
BASO%, median (IQR)	0.004(0.003, 0.006)	0.005 (0.003, 0.006)	0.004 (0.003, 0.006)	<0.001
RDW, median (IQR)	15.0 (14.0, 18.0)	15.0 (14.0, 18.0)	16.0 (14.0, 18.0)	<0.001
PCT, median (IQR)	0.240 (0.190, 0.301)	0.230 (0.180, 0.280)	0.260 (0.200, 0.335)	<0.001
Biochemical indices
TC, median (IQR)	4.46 (3.75, 5.24)	4.40 (3.74, 5.14)	4.51 (3.76, 5.39)	<0.001
DBIL, median (IQR)	2.80 (2.00, 4.00)	2.90 (2.05, 4.10)	2.70 (1.90, 3.90)	<0.001
IBIL, median (IQR)	5.50 (3.60, 8.10)	5.80 (3.90, 8.50)	5.10 (3.30, 7.60)	<0.001
DTR, median (IQR)	0.34 (0.27, 0.42)	0.33 (0.27, 0.41)	0.34 (0.27, 0.43)	0.034
TP, median (IQR)	66.5 (61.0, 71.9)	66.4 (61.2, 71.8)	66.7 (60.6, 72.1)	0.005
ALB, median (IQR)	35.1 (31.4, 39.1)	35.4 (32.1, 39.6)	34.5 (30.3, 38.4)	0.003
GLOB, median (IQR)	30.8 (26.9, 35.1)	30.3 (26.5, 34.4)	31.4 (27.5, 36.1)	<0.001
AGR, median (IQR)	1.10 (0.90, 1.40)	1.20 (1.00, 1.40)	1.10 (0.90, 1.30)	0.023
GGT, median (IQR)	22.0 (14.0, 37.0)	20.0 (13.0, 33.0)	24.0 (16.0, 44.0)	0.176
TBA, median (IQR)	5.30 (3.10, 9.10)	5.20 (3.10, 8.90)	5.50 (3.10, 9.30)	0.059
AST, median (IQR)	23.0 (18.0, 30.0)	23.0 (18.0, 30.0)	23.0 (18.0, 32.0)	0.002
ALT, median (IQR)	15.0 (11.0, 24.0)	15.0 (11.0, 24.0)	16.0 (11.0, 24.0)	<0.001
AAR, median (IQR)	1.50 (1.10, 2.00)	1.40 (1.00, 2.00)	1.50 (1.10, 2.10)	0.395
ALP, median (IQR)	80.0 (62.0, 104.0)	77.0 (60.0, 102.0)	83.0 (65.0, 107.0)	0.557
PA, median (IQR)	205.3 (154.7, 256.7)	212.1 (166.2, 262.6)	195.7 (139.8, 247.0)	<0.001
UREA, median (IQR)	4.75 (3.68, 6.17)	4.71 (3.68, 6.03)	4.82 (3.69, 6.33)	<0.001
Cr, median (IQR)	61.0 (50.0, 75.0)	59.0 (49.0, 73.0)	63.0 (51.0, 79.0)	<0.001
UA, median (IQR)	283 (221,362)	280.0 (220.0, 354.0)	287.0 (221.0, 369.0)	<0.001
HCO₃^-, median (IQR)	23.9 (22.0 26.0)	24.0 (22.1, 26.0)	23.8 (21.8, 26.0)	<0.001
Ccr, median (IQR)	71.5 (53.9, 90.2)	75.0 (58.5, 92.3)	65.9 (47.5, 84.5)	<0.001
CysC, median (IQR)	1.05 (0.86, 1.36)	1.01 (0.84, 1.26)	1.14 (0.91, 1.55)	<0.001
hsTn, median (IQR)	11.30 (4.93, 25.80)	10.60 (4.86, 24.05)	11.80 (5.53, 26.94)	<0.001
CK, median (IQR)	57.0 (35.0, 92.0)	56.0 (33.0, 89.0)	58.0 (36.0, 96.0)	<0.001
CK-MB, median (IQR)	13.0 (9.5, 16.0)	12.3 (10.0, 16.0)	13.0 (9.0, 16.6)	<0.001
LDH, median (IQR)	193.0 (160.0, 236.0)	187.0 (156.0, 227.0)	204.0 (165.0, 255.0)	<0.001
HBDH, median (IQR)	141.0 (117.0, 176.0)	137.0 (115.0, 169.0)	147.0 (120.0, 189.0)	<0.001
RBG, median (IQR)	5.52 (4.71, 6.90)	5.63(4.75, 7.08)	5.47 (4.69, 6.71)	<0.001
ESR, median (IQR)	34.0 (16.0, 62.0)	26.0 (13.0, 49.0)	51.0 (23.0, 82.0)	<0.001
NLR, median (IQR)	2.72 (1.77, 4.35)	2.44 (1.59, 3.70)	3.17 (2.05, 5.37)	<0.001
MLR, median (IQR)	0.34 (0.24, 0.49)	0.31 (0.22, 0.44)	0.38 (0.26, 0.56)	<0.001
PLR, median (IQR)	180.1 (128.7, 260.6)	156.7 (116.8, 213.4)	232.2 (160.3, 341.1)	<0.001

Data presented as median (IQR) for continuous variables and number (%) for categorical variables.

AAR, aspartate-to-alanine aminotransferase ratio; AGR, albumin-to-globulin ratio; ALB, albumin; ALP, alkaline phosphatase; ALT, alanine aminotransferase; AST, aspartate aminotransferase; BASO, basophil count; Ccr, endogenous creatinine clearance; CK, creatine kinase; CK-MB, creatine kinase isoenzyme; Cr, creatinine; CysC, cystatin C; DBIL, direct bilirubin; DTR, direct-to-total bilirubin ratio; EO, eosinophil count; ESR, erythrocyte sedimentation rate; GGT, gamma-glutamyl transferase; GLOB, globulin; HBDH, alpha-hydroxybutyrate dehydrogenase; $H C O_{3}^{-}$ , bicarbonate; HCT, hematocrit percentage; HGB, hemoglobin concentration; hsTn, high-sensitivity troponin; IBIL, indirect bilirubin; IQR, interquartile range; LDH, lactate dehydrogenase; LYMPH, lymphocyte count; LYMPH%, lymphocyte percentage; MCH, mean corpuscular hemoglobin; MCHC, mean corpuscular hemoglobin concentration; MCV, mean corpuscular volume; MLR, monocyte-to-lymphocyte ratio; MONO, monocyte count; MPV, mean platelet volume; NEUT, neutrophil count; NEUT%, neutrophil percentage; NLR, neutrophil-to-lymphocyte ratio; PA, prealbumin; PCT, platelet crit; PDW, platelet distribution width; PLR, platelet-to-lymphocyte ratio; PLT, platelet count; RA, rheumatoid arthritis; RBC, red blood cell count; RBG, random blood glucose; RDW, red blood cell distribution width; TBA, total bile acids; TC, total cholesterol; TP, total protein; UA, uric acid; UREA, urea; WBC, white blood cell count; Bold values indicate statistical significance (p < 0.05).

Performance evaluation and feature selection of machine learning models for rheumatoid arthritis severity prediction

The RF model, trained on 55 clinical variables without any feature selection or extraction, exhibited outstanding predictive performance on the test set, with an AUC of 0.870. Figure 3 presents a comparison of model performance, calibration curves, and decision curves. Additional evaluation metrics are provided in Table S5 of the supporting information. Given the superior performance of the Random Forest model compared to other models, its RFE effectively identified an optimal feature set comprising 15 variables (Figure 2). Beyond 15 features, the cross-validation error (0.20) and accuracy (0.75–0.78) stabilize, indicating diminishing returns and supporting the use of a concise feature subset. Variables such as Age, LYMPH%, PDW, PLR, ESR, RBC, HCT, PLT, MPV, NEUT%, PCT, LYMPH, MCV, ALP, and NLR were identified as key factors driving the model’s final output. Figure 3 shows a comparison of the performance, calibration curves, and decision curves for each final model using these improved variable sets. The RF models achieved the highest AUC of 0.874 on the test set. Other evaluation indicators are listed in Table S6 in the supporting information. We conducted a sensitivity analysis by excluding ESR and reassessed the model’s discrimination and calibration, with results remaining consistent (Figure S1 and Table S7). On the basis of these results, we chose the RF model as the optimal model for further analysis.

Figure 2.

Cross-validation performance and feature importance from recursive feature elimination with Random Forest. The left panel displays the 10-fold CV error versus the number of features, with the dashed red line indicating the minimum error threshold. The middle panel shows the 10-fold CV accuracy versus the number of features, with the dashed red line marking the maximum accuracy threshold. The right panel presents the top 15 feature importances from RFE, ranked on the y-axis by significance, with higher-ranked features indicating greater importance.

Figure 3.

The performance and comparison of 10 different predictive models. (a) ROC curves for the training set using all features. (b) Calibration curves for the training set using all features. (c) Decision curve analysis for the training set using all features. (d) ROC curves for the test set using all features. (e) Calibration curves for the test set using all features. (f) Decision curve analysis for the test set using all features. (g) ROC curves for the training set using the top 15 most important features. (h) Calibration curves for the training set using the top 15 most important features. (i) Decision curve analysis for the training set using the top 15 most important features. (j) ROC curves for the test set using the top 15 most important features. (k) Calibration curves for the test set using the top 15 most important features. (l) Decision curve analysis for the test set using the top 15 most important features.

Explainable machine learning model

The SHAP feature importance plot and SHAP summary plot (Figure 4) elucidate the interpretability of the RF model, highlighting the 15 most influential features. These key determinants of the model’s final output predominantly encompass blood cell indices and their potential derivatives, alongside demographic characteristics. Notably, 10 critical hematological parameters—LYMPH%, PDW, RBC, HCT, PLT, MPV, NEUT%, PCT, LYMPH, and MCV—emerge as central contributors, supplemented by the effects of Age, ESR, PLR, ALP, and NLR. Among these, LYMPH% and PDW exhibit the highest mean SHAP values, underscoring their significant impact on the model’s predictions. Therefore, we further analyzed the key blood cell indicators.

Figure 4.

Global interpretation of the RF model and feature importance. The SHAP summary plot shows key features affecting the RF model’s predictions for rheumatoid arthritis severity. Features are ranked on the y-axis by significance, with higher-ranked features indicating greater importance. Each dot represents a patient, with its color indicating the feature value—blue for lower and red for higher values. The dot’s position on the x-axis (SHAP value) reflects the feature’s impact on the model’s output, indicating the likelihood of a clinical prediction, with higher SHAP values corresponding to a higher probability of the more severe rheumatoid arthritis.

From the waterfall plot of randomly selected patients (Figure 5), it can be seen that the order of influence of each parameter varies between different patients. The predicted value of sample 11 (Figure 5(a)) is 0.03, which is lower than the preliminary predicted value of 0.425, mainly due to the combined effect of negative contributions such as “PLT” (−0.06), “PLR” (−0.06), and “PDW” (−0.04). The predicted value of sample 1248 (Figure 5(b)) is 0.09, and the preliminary predicted value is 0.425, which is mainly affected by “PLT” (−0.05), “LYMPH%” (−0.05), “ESR” (−0.05) and the predicted value of sample 3341 (Figure 5(c)) is 0.64, and the preliminary expected value is 0.425, which is mainly driven by the positive and negative contributions of “PDW” (+0.25) and “LYMPH%” (−0.05). The predicted value of sample 4145 (Figure 5(d)) is 0.93, and the preliminary expected value is 0.425, which is mainly driven by positive contributions of “PLR” (+0.16) and “LYMP%” (+0.09). The SHAP waterfall chart analysis shows that indicators such as LYMPH%, PDW, ESR, and RBC have a significant impact on the prediction results. For further study, we use SHAP dependency diagram to analyze.

Figure 5.

SHAP waterfall plots demonstrating feature contributions to model predictions for individual samples. Panels (a) to (d) correspond to samples with 11, 1248, 3341, and 4145 patients, respectively. Each panel displays features including LYMPH%, ESR, PDW, PLT, and others, with associated SHAP values. Red bars indicate positive SHAP values, increasing the likelihood of the predicted clinical outcome, while blue bars represent negative SHAP values, decreasing the likelihood. The x-axis shows the SHAP value (E[f(X)]), reflecting the impact of features on the prediction, with bar length indicating the magnitude of influence. The f(X) value at the top of each panel denotes the model’s predicted output for the sample, and the E[f(X)] value at the bottom represents the baseline expected output.

The influence of variables on the model’s output is shown in SHAP dependence plots (Figure 6). LYMPH% (Figure 6(a)), lower values (0.4) show negative SHAP values, suggesting a reduced likelihood of severe RA, with a transition around 0.2–0.4 where the influence shifts. PDW (Figure 6(b)) exhibits a positive contribution to severe RA likelihood, with SHAP values increasing sharply beyond 15, peaking around 20–25, indicating that higher PDW values are strongly associated with severe RA. ESR (Figure 6(c)) shows a non-linear trend: SHAP values are positive and increase with ESR up to approximately 60–80, suggesting a higher likelihood of severe RA, but plateau or slightly decline beyond 100, indicating a complex relationship with severity. RBC (Figure 6(d)) has a negative influence on severe RA likelihood, with SHAP values decreasing as RBC values rise above 4, particularly becoming strongly negative beyond 5, reflecting a protective effect of higher RBC levels. These patterns reflect key thresholds in the model’s interpretation of RA severity. Regarding the clinical significance and practical application of the selected variables, we emphasize that the final model relies on routinely available, low-cost hematological and demographic indicators (e.g., NEUT, LYMPH, Age, HGB), which are commonly assessed in clinical settings and may support severity stratification for RA in resource-limited environments. Furthermore, to improve the clinical utility of the model, we calculated additional diagnostic evaluation metrics, including Threshold, Accuracy, Sensitivity, Specificity, Precision, F1, and Brier, in both internal validation sets. In the all-variable model evaluation, the RF model achieved a F1 of 76.8% and a Brier of 14.8%. And in the evaluation restricted to the top 15 variables, it achieved a F1 of 77.5% and a Brier of 14.0% (Tables S4 and S5). These results highlight the strong predictive performance of the model and underscore its potential for future clinical integration, pending confirmation in prospective validation studies.

Figure 6.

The SHAP dependence plot shows the effect of individual features on the prediction model’s output. Panels (a) through (d) depict features like LYMPH%, PDW, ESR, and RBC. Positive SHAP values indicate a higher probability of a clinical diagnosis. The x-axis represents the feature value, and the y-axis shows its impact on the predicted outcome. Each point corresponds to a prediction for a specific patient.

Partial dependence plots offer valuable insights into the nonlinear relationships effects within predictive models. As illustrated in Figure 6(a)–(d), the features LYMPH%, PDW, ESR, and RBC demonstrate nonlinear influences on RA severity. These four main indices were chosen for analysis due to their significant impact on the model’s prediction results. Figure 7(a)–(c) elucidates the interactions among LYMPH, NEUT%, PLT, and RBC: in Figure 7(a), the predicted probability of RA is elevated (0.48–0.50) when LYMPH is low (1.0–1.5 × 10⁹/L) and NEUT% is low (0.45–0.50), but it decreases to 0.38–0.42 when LYMPH rises to 2.0–2.5 × 10⁹/L and NEUT% exceeds 0.60; Figure7(b) shows an increased probability (0.60–0.63) with low LYMPH (1.0–1.5 × 10⁹/L) and high PLT (400–500 × 10⁹/L), dropping to 0.36–0.42 when LYMPH is high (2.0–2.5 × 10⁹/L) and PLT is low (200–250 × 10⁹/L); Figure 7(c) reveals a probability range of 0.57–0.63 when PLT is high (400–500 × 10⁹/L) and NEUT% is low (0.45–0.50), decreasing to 0.36–0.39 when PLT is low (200–250 × 10⁹/L) and NEUT% is high (0.70–0.80). Figure 7(d)–(f) further illustrates additional interaction effects: Figure 7(d) indicates a higher probability (0.64–0.68) with low LYMPH (1.0–1.5 × 10⁹/L) and high RBC (4.5–5.0 × 10¹²/L), falling to 0.36–0.40 when LYMPH is high (2.0–2.5 × 10⁹/L) and RBC is low (3.0–3.5 × 10¹²/L); Figure 7(e) demonstrates the highest probability (0.75–0.80) with high PLT (400–500 × 10⁹/L) and high RBC (4.5–5.0 × 10¹²/L), dropping to 0.35–0.40 when both are low (200–250 × 10⁹/L and 3.0–3.5 × 10¹²/L); Figure 7(f) shows a probability of 0.64–0.68 with low LYMPH (1.0–1.5 × 10⁹/L) and high RBC (4.5–5.0 × 10¹²/L), decreasing to 0.36–0.40 when LYMPH is high (2.0–2.5 × 10⁹/L) and RBC is low (3.0–3.5 × 10¹²/L). Collectively, the synergistic interaction between PLT and RBC, particularly at PLT levels of 400–500 × 10⁹/L and RBC levels of 4.5–5.0 × 10¹²/L, emerges as the significant predictor of RA severity, with lower LYMPH and NEUT% levels further enhancing this effect. Consequently, monitoring and managing PLT, RBC, LYMPH, and NEUT% levels may hold critical importance for mitigating RA severity. In addition to probability-based discrimination metrics, we identified an optimal decision threshold based on the Youden index to facilitate clinical interpretation. Application of this threshold yielded high sensitivity and specificity on the internal test set, demonstrating the model’s potential utility for risk stratification in clinical settings. Detailed threshold values and associated metrics are provided in Table S2.

Figure 7.

Synergistic effects of hematological indices on rheumatoid arthritis severity. PDPs elucidate the synergistic interactions between hematological indices and the severity of rheumatoid arthritis, encompassing (a) NEUT% and LYMPH, (b) PLT and LYMPH, (c) NEUT% and PLT, (d) NEUT% and RBC, (e) RBC and PLT, and (f) RBC and LYMPH. The color bar on the right provides an intuitive visualization of the predicted severity values for rheumatoid arthritis derived from model training. The x-axis and y-axis represent clinical data for the respective hematological indices, while the contour plots vividly illustrate the combined influence of the two indices on the predicted severity of rheumatoid arthritis.

Associations between blood cell indices and the severity of incident rheumatoid arthritis

The relationships between 10 significant blood cell indices and the severity of RA were examined using binary logistic regression model I. In the primary analysis, all indicators were significantly associated with the severity of incident RA (Table 2). Increased PDW (OR 2.17, 95% CI: 2.01–2.35, p < 0.001), PLT (OR 1.96, 95% CI: 1.81–2.13, p < 0.001), PCT (OR 1.61, 95% CI: 1.49–1.73, p < 0.001), and NEUT% (OR 1.47, 95% CI: 1.32–1.63, p < 0.001) were associated with a higher risk of severity of RA, while increased HCT (OR 0.11, 95% CI: 0.05–0.24, p < 0.001), RBC (OR 0.59, 95% CI: 0.53–0.65, p < 0.001), MPV (OR 0.70, 95% CI: 0.65–0.75, p < 0.001), LYMPH% (OR 0.54, 95% CI: 0.50–0.59, p < 0.001), LYMPH (OR 0.62, 95% CI: 0.57–0.67, p < 0.001), and MCV (OR 0.91, 95% CI: 0.84–0.97, p < 0.001) were protective against severity of RA. Sensitivity analyses confirmed that the relationship between all indicators and severity of RA remained significant in Model II (Table 2).

Table 2.

Associations between seven core blood count indices selected based on the SHAP feature importance plot and the severity of rheumatoid arthritis.

Blood count indices	Model I		Model II
Blood count indices	OR (95% CI)	p-Value	OR (95% CI)	p-Value
PDW	2.17 (2.01–2.35)	< 0.001	2.13 (1.96–2.32)	<0.001
HCT	0.11 (0.05–0.24)	<0.001	0.12 (0.05–0.25)	<0.001
RBC	0.59 (0.53–0.65)	<0.001	0.60 (0.54–0.67)	<0.001
MPV	0.70 (0.65–0.75)	<0.001	0.74 (0.69–0.80)	<0.001
PLT	1.96 (1.81–2.13)	<0.001	1.86 (1.70–2.04)	<0.001
LYMPH%	0.54 (0.50–0.59)	<0.001	0.58 (0.53–0.63)	<0.001
PCT	1.61 (1.49–1.73)	<0.001	1.52 (1.40–1.65)	<0.001
LYMPH	0.62 (0.57–0.67)	<0.001	0.64 (0.59–0.71)	<0.001
NEUT%	1.47 (1.32–1.63)	<0.001	1.31 (1.17–1.46)	<0.001
MCV	0.91 (0.84–0.97)	0.008	0.85 (0.78–0.92)	<0.001

Model I: Adjusted for age, sex, hemoglobin concentration, mean corpuscular hemoglobin concentration, neutrophil count, monocyte count, monocyte percentage, eosinophil count, eosinophil percentage, basophil count, basophil percentage, red blood cell (erythrocyte) distribution width. Model II: Model I with additional adjustment for total cholesterol, direct bilirubin, indirect bilirubin, direct to total bilirubin ratio, total protein, albumin to globulin ratio, gamma-glutamyl transferase, total bile acids, aspartate aminotransferase, alanine aminotransferase, aspartate aminotransferase to alanine aminotransferase ratio, alkaline phosphatase, prealbumin, urea, creatinine, uric acid, bicarbonate, endogenous creatinine clearance, cystatin C, high-sensitivity troponin, creatine kinase, creatine kinase isoenzyme, alpha-hydroxybutyrate dehydrogenase, random blood glucose, and erythrocyte sedimentation rate.

CI, confidence interval; HCT, hematocrit percentage; LYMPH, lymphocyte count; LYMPH%, lymphocyte percentage; MCV, mean corpuscular volume; MPV, mean platelet volume; NEUT%, neutrophil percentage; OR, odds ratio; PCT, platelet crit; PDW, platelet distribution width; PLT, platelet count; RBC, red blood cell count; SHAP, Shapley Additive Explanations; Bold values indicate statistical significance (p < 0.05).

Using restricted cubic spline, we found that PDW, LYMPH%, HCT, PCT, RBC, LYMPH, MPV, and MCV blood cell indicators had significant nonlinear relationships with the severity of RA (p < 0.05), consisting of J-shape, inversion of a J-shape, and U-shape curve associations, but NEUT% and PLT had linear relationships with the severity of RA (Figure 8). The results from binary logistic regression and restricted cubic splines (RCS) align with the SHAP interpretability outcomes of the machine learning model, particularly in terms of feature importance, directionality, and nonlinear patterns, with key indicators such as LYMPH%, LYMPH, PDW, and PLT showing consistent trends. While logistic regression and RCS cannot directly validate SHAP results, their similar patterns provide corroborative evidence and complement the machine learning approach, fostering a more comprehensive understanding of data relationships. This integration enhances model interpretability and offers practical utility for stratification of RA severity, especially in resource-constrained clinical settings.

Figure 8.

Nonlinear associations of hematological indices with rheumatoid arthritis severity fitted by logistic regression. Panels present RCS plots for MPV, RBC, MCV, PLT, LYMPH, LYMPH%, HCT, PDW, NEUT%, and PCT, with p-values indicating the significance of nonlinearity (p < 0.05). Red dashed lines represent the fitted curves, black dashed lines denote reference levels, and shaded areas indicate 95% confidence intervals. The x-axes display values of the hematological indices, while the y-axes represent the probability of severity.

Causal associations between blood cell indices with RA risk

MR analysis revealed causal associations between 10 blood cell indices and RA risk, as detailed in Table S8. Protective effects were observed for MCV (OR = 0.90, p < 0.0001 for both IVW-FE and IVW-RE), MPV (OR = 0.94, p = 0.0020 IVW-FE; p = 0.0195 IVW-RE), NEUT% (OR = 0.81, p < 0.0001 IVW-FE; p = 0.0136 IVW-RE), as well as slight protection from LYMPH, LYMP%, and PDW (OR = 0.90–0.94, p < 0.05 IVW-FE), though high heterogeneity and non-significant IVW-RE results suggest weak evidence requiring additional SNPs to address variability; PCT also showed potential risk reduction (weighted median OR = 0.84, p = 0.0384) but with limited support. Conversely, risk-increasing effects were linked to HCT (OR = 1.20, p = 0.0005 IVW-FE; p = 0.0361 IVW-RE), RBC (OR = 1.12, p < 0.0001 IVW-FE; p = 0.0098 IVW-RE), and PLT (OR = 1.08, p = 0.0010 IVW-FE; p = 0.0398 IVW-RE). These findings partially align with machine learning, binary logistic regression, and RCS analyses of RA severity, demonstrating consistent protective roles for MCV, PLT, LYMPH, and LYMPH%, yet directional discrepancies in HCT, RBC, and PDW underscore distinct pathological mechanisms between RA risk and severity. SHAP and RCS analyses further highlighted nonlinear effects across index ranges, complementing the linear MR estimates. All instruments exhibited F-statistics >10 for validity, with Cochran’s Q tests indicating no overall heterogeneity except localized elevations in LYMPH, LYMP%, and PDW; MR-Egger intercept and MR-PRESSO analyses confirmed the absence of horizontal pleiotropy, supporting the robustness of these causal inferences.

Discussion

RA is a chronic, progressive, immune-mediated inflammatory disease that primarily affects synovial joints, causing bone erosion and significant disability.^41,42 Accurate assessment of disease severity remains challenging, as it typically requires composite indices, laboratory biomarkers, and imaging evaluations, which are costly and not always available in resource-limited settings.^43,44 In this context, peripheral blood indices, which are inexpensive, widely accessible, and routinely measured, represent promising tools for evaluating RA severity.

Existing studies have established associations between hematological markers, such as NLR and RDW, and RA activity.^45,46 However, these investigations often focus on isolated markers or acute-phase reactants like CRP and ESR, lacking comprehensive integration of blood cell indices for severity prediction. Moreover, causal relationships and nonlinear associations remain underexplored, limiting their clinical utility. Building on these insights, our study employed ML, SHAP, binary logistic regression, RCS, and MR to elucidate the role of blood cell indices in RA progression.

Our retrospective case–control analysis of 4401 RA patients identified significant associations between 10 key blood cell indices (e.g., LYMPH%, PDW, MCV) and RA severity. Binary logistic regression revealed that elevated PDW, PLT, and NEUT% were associated with higher severity, while HCT, RBC, and MCV were protective. RCS analyses confirmed nonlinear relationships, such as J-shaped or U-shaped curves for PDW and LYMPH%. The optimized RF model, using 15 selected features, achieved an AUC of 0.874 on the test set, with robust metrics including F1 (77.5%) and Brier (14.0%). SHAP analyses highlighted LYMPH% and PDW as top contributors, revealing patient-specific and synergistic effects.

Notably, MR analysis indicated that genetically higher HCT (OR = 1.20, p < 0.001) and RBC (OR = 1.12, p < 0.001) increased RA risk, contrasting with their protective role in severity models. This directional discrepancy likely stems from differing contexts: MR captures lifelong susceptibility, whereas severity models reflect baseline associations in established RA, where anemia of chronic inflammation reduces HCT and RBC in severe cases.^47,48 These findings underscore distinct mechanisms between RA onset and progression.

The observed associations may involve immunological and inflammatory pathways. Elevated NEUT% could exacerbate severity through proinflammatory cytokines (e.g., IL-6, TNF-α), while higher LYMPH% may promote immune balance. Platelet indices like PDW and PLT contribute to synovial inflammation via cytokine release,⁴⁹ and red cell parameters influence tissue oxygenation.⁵⁰ These align with prior evidence linking hematological indices to RA activity.

Clinically, abnormalities in indices such as elevated PDW or reduced LYMPH% could serve as cost-effective markers for RA severity stratification. Integrated with clinical symptoms and biomarkers (e.g., ESR, CRP), they may guide tailored treatments, optimizing resource allocation in constrained settings. Future prospective studies should validate these markers across diverse populations, incorporating disease controls for specificity.

A limitation of this study is the potential construct overlap between predictors (e.g., ESR, CRP, and blood cell ratios) and the outcome measures (DAS28-ESR and CDAI), which may introduce same-source bias. While our analyses suggest that the main results are robust, future studies should exclude inflammation-based indices from predictor sets or employ alternative severity definitions to further minimize bias. Another limitation of this study is the lack of RA-specific serological biomarkers, such as rheumatoid factor and anti-CCP, as well as standardized composite indices (e.g., SDAI). These markers are widely used in clinical practice and are important for diagnosis and prognosis in RA. Their absence may limit the generalizability of our findings and the direct applicability of our models to clinical decision-making. Future studies incorporating these RA-specific biomarkers are warranted to further validate and enhance the predictive value of blood cell indices. A potential limitation of this study is treatment confounding. Therapeutic interventions such as conventional synthetic DMARDs, biologic DMARDs, targeted synthetic DMARDs, and surgery can affect hematological indices and are also related to disease severity. Although we minimized this issue by defining severity solely on DAS28-ESR and CDAI, residual confounding by treatment cannot be completely excluded. Future prospective studies should explicitly model treatment variables to further disentangle their effects. In addition, the single-center and retrospective nature of this study inherently limits the external generalizability of our findings. All patients were recruited from a single tertiary hospital, and the demographic as well as clinical characteristics may not fully represent the broader RA population in other geographic regions or healthcare systems. Moreover, although the machine learning models demonstrated excellent internal performance, their predictive ability may not directly translate to other institutions with different patient profiles, laboratory platforms, or clinical practice patterns. Future work will therefore include external validation in multicenter cohorts, prospective data collection, and model recalibration to improve robustness and applicability across diverse clinical settings.

Conclusion

This study developed and internally validated a prediction model for RA severity. Interpretable machine learning techniques were used to reveal important associations between key factors and RA severity, and the final model showed strong predictive accuracy. In addition, this population-based retrospective study explored the association between blood cell indices and RA severity, providing new insights into the molecular mechanisms behind RA disease severity stratification.

Supplemental Material

sj-doc-2-tab-10.1177_1759720X251407065 – Supplemental material for Associations of blood cell indices with the severity of rheumatoid arthritis: a retrospective case–control and machine learning modeling study

Supplemental material, sj-doc-2-tab-10.1177_1759720X251407065 for Associations of blood cell indices with the severity of rheumatoid arthritis: a retrospective case–control and machine learning modeling study by Rongqing He, Boli Qin, Xiaopeng Qin, Jiarui Chen, Jiang Xue, Tianyou Chen, Songze Wu, Zhuo Chen, Xinli Zhan and Chong Liu in Therapeutic Advances in Musculoskeletal Disease

Supplemental Material

sj-docx-3-tab-10.1177_1759720X251407065 – Supplemental material for Associations of blood cell indices with the severity of rheumatoid arthritis: a retrospective case–control and machine learning modeling study

Supplemental material, sj-docx-3-tab-10.1177_1759720X251407065 for Associations of blood cell indices with the severity of rheumatoid arthritis: a retrospective case–control and machine learning modeling study by Rongqing He, Boli Qin, Xiaopeng Qin, Jiarui Chen, Jiang Xue, Tianyou Chen, Songze Wu, Zhuo Chen, Xinli Zhan and Chong Liu in Therapeutic Advances in Musculoskeletal Disease

Supplemental Material

sj-docx-4-tab-10.1177_1759720X251407065 – Supplemental material for Associations of blood cell indices with the severity of rheumatoid arthritis: a retrospective case–control and machine learning modeling study

Supplemental material, sj-docx-4-tab-10.1177_1759720X251407065 for Associations of blood cell indices with the severity of rheumatoid arthritis: a retrospective case–control and machine learning modeling study by Rongqing He, Boli Qin, Xiaopeng Qin, Jiarui Chen, Jiang Xue, Tianyou Chen, Songze Wu, Zhuo Chen, Xinli Zhan and Chong Liu in Therapeutic Advances in Musculoskeletal Disease

Research Data

sj-xlsx-1-tab-10.1177_1759720X251407065 – for Associations of blood cell indices with the severity of rheumatoid arthritis: a retrospective case–control and machine learning modeling study

sj-xlsx-1-tab-10.1177_1759720X251407065 for Associations of blood cell indices with the severity of rheumatoid arthritis: a retrospective case–control and machine learning modeling study by Rongqing He, Boli Qin, Xiaopeng Qin, Jiarui Chen, Jiang Xue, Tianyou Chen, Songze Wu, Zhuo Chen, Xinli Zhan and Chong Liu in Therapeutic Advances in Musculoskeletal Disease

This article is distributed under the terms of the Creative Commons Attribution 4.0 License (http://www.creativecommons.org/licenses/by/4.0/) which permits any use, reproduction and distribution of the work without further permission provided the original work is attributed as specified on the SAGE and Open Access pages (https://us.sagepub.com/en-us/nam/open-access-at-sage).

Footnotes

Acknowledgements

We are very grateful to Prof. Chong Liu and Prof. Xinli Zhan (Spine and Osteopathy Ward, The First Affiliated Hospital of Guangxi Medical University) for support in all stages of this study.

Declarations

ORCID iD

Chong Liu

Supplemental material

Supplemental Material for this article is available online.

References

Smolen

Aletaha

Barton

, et al. Rheumatoid arthritis. Nat Rev Dis Primers 2018; 4: 18001.

Gonzalez-Rey

Gonzalez

Varela

, et al. Human adipose-derived mesenchymal stem cells reduce inflammatory and T cell responses and induce regulatory T cells in vitro in rheumatoid arthritis. Ann Rheum Dis 2010; 69(1): 241–248.

Smolen

Aletaha

McInnes

. Rheumatoid arthritis. Lancet 2016; 388(10055): 2023–2038.

Stojanovic

Stamenkovic

Cvetkovic

, et al. Matrix metalloproteinase-9 level in synovial fluid-association with joint destruction in early rheumatoid arthritis. Medicina (Kaunas) 2023; 59(1): 167.

Orr

Najm

Young

, et al. The utility and limitations of CRP, ESR and DAS28-CRP in appraising disease activity in rheumatoid arthritis. Front Med (Lausanne) 2018; 5: 185.

Smolen

Landewé

RBM

Bijlsma

JWJ

, et al. EULAR recommendations for the management of rheumatoid arthritis with synthetic and biological disease-modifying antirheumatic drugs: 2019 update. Ann Rheum Dis 2020; 79(6): 685–699.

Cheng

Cai

Xiang

, et al. Platelet-to-lymphocyte ratios as a haematological marker of synovitis in rheumatoid arthritis with normal acute phase reactant level. Ann Med 2024; 56(1): 2346546.

Zhou

, et al. The neutrophil-lymphocyte ratio predicts all-cause and cardiovascular mortality among U.S. adults with rheumatoid arthritis: results from NHANES 1999–2020. Front Immunol 2023; 14: 1309835.

Alghamdi

. Red blood cell distribution width: a potential inexpensive marker for disease activity in patients with rheumatic diseases; scoping review. Open Access Rheumatol 2023; 15: 173–180.

10.

Chen

Lei

Liu

, et al. Association between red blood cell distribution width and genetic risk in patients with rheumatoid arthritis: a prospective cohort study and Mendelian randomization analysis. BMC Rheumatol 2025; 9(1): 1.

11.

Sargin

Senturk

Yavasoglu

, et al. Relationship between neutrophil-lymphocyte, platelet-lymphocyte ratio and disease activity in rheumatoid arthritis treated with rituximab. Int J Rheum Dis 2018; 21(12): 2122–2127.

12.

Liu

Sun

, et al. The association between neutrophil-to-lymphocyte ratio and disease activity in rheumatoid arthritis. Inflammopharmacology 2023; 31(5): 2237–2244.

13.

von Elm

Altman

Egger

, et al. The Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) statement: guidelines for reporting observational studies. Ann Intern Med 2007; 147(8): 573–577.

14.

von Elm

Altman

Egger

, et al. Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) statement: guidelines for reporting observational studies. BMJ 2007; 335(7624): 806–808.

15.

Aletaha

Neogi

Silman

, et al. 2010 rheumatoid arthritis classification criteria: an American College of Rheumatology/European League Against Rheumatism collaborative initiative. Ann Rheum Dis 2010; 69(9): 1580–1588.

16.

Deighton

O’Mahony

Tosh

, et al. Management of rheumatoid arthritis: summary of NICE guidance. BMJ 2009; 338: b702.

17.

Binder

Pfisterer

Lang

, et al. Mlr3pipelines-flexible machine learning pipelines in r. J Mach Learn Res 2021; 22(184): 1–7.

18.

Bischl

Sonabend

Kotthoff

, et al. Applied machine learning using mlr3 in R. Boca Raton, FL: CRC Press, 2024.

19.

Binder

Pfisterer

Schneider

, et al. mlr3pipelines: Preprocessing Operators and Pipelines for ‘mlr3’. R package version 0.3 1, 2020.

20.

Matsuo

Kamada

Imamura

, et al. Machine learning-based prediction of relapse in rheumatoid arthritis patients using data on ultrasound examination and blood test. Sci Rep 2022; 12(1): 7224.

21.

Kohavi

(ed.). A study of cross-validation and bootstrap for accuracy estimation and model selection. Montreal, QC: IJCAI, 1995.

22.

Yang

Liu

L-L

Huang

, et al. Machine learning-based landslide susceptibility assessment with optimized ratio of landslide to non-landslide samples. Gondwana Res 2023; 123: 198–216.

23.

Amin

Imtiaz

Haroon

, et al. Assessment of machine learning algorithms for land cover classification in a complex mountainous landscape. J Geovis Spat Anal 2024; 8(2): 34.

24.

Yesilyaprak

Kumar

Agrawal

, et al. Predicting long-term clinical outcomes of patients with recurrent pericarditis. J Am Coll Cardiol 2024; 84(13): 1193–1204.

25.

Lundberg

. A unified approach to interpreting model predictions. arXiv preprint. arXiv:170507874. 2017.

26.

Scott

Su-In

. A unified approach to interpreting model predictions. Adv Neural Inf Process Syst 2017; 30: 4765–4774.

27.

Kraev

Koseoglu

Traverso

, et al. Shap-select: lightweight feature selection using SHAP values and regression. arXiv preprint. arXiv:241006815. 2024.

28.

Mihirette

Tan

(eds.). SHAP algorithm for healthcare data classification. International conference on hybrid artificial intelligence systems 2022; 363–374.

29.

Vellido

. The importance of interpretability and visualization in machine learning for applications in medicine and health care. Neural Comput Appl 2020; 32(24): 18069–18083.

30.

Apley

Zhu

. Visualizing the effects of predictor variables in black box supervised learning models. J R Stat Soc Series B Stat Methodol 2020; 82(4): 1059–1086.

31.

Nohara

Matsumoto

Soejima

, et al. Explanation of machine learning models using shapley additive explanation and application for real data in hospital. Comput Methods Programs Biomed 2022; 214: 106584.

32.

Ponce-Bobadilla

Schmitt

Maier

, et al. Practical guide to SHAP analysis: explaining supervised machine learning model predictions in drug development. Clin Transl Sci 2024; 17(11): e70056.

33.

Wootton

Lawn

Millard

LAC

, et al. Evaluation of the causal effects between subjective wellbeing and cardiometabolic health: Mendelian randomisation study. BMJ 2018; 362: k3788.

34.

Bowden

Davey Smith

Haycock

, et al. Consistent estimation in Mendelian randomization with some invalid instruments using a weighted median estimator. Genet Epidemiol 2016; 40(4): 304–314.

35.

Burgess

Thompson

. Interpreting findings from Mendelian randomization using the MR-Egger method. Eur J Epidemiol 2017; 32(5): 377–389.

36.

Slob

EAW

Burgess

. A comparison of robust Mendelian randomization methods using summary data. Genet Epidemiol 2020; 44(4): 313–329.

37.

Gkatzionis

Burgess

Newcombe

. Statistical methods for cis-Mendelian randomization with two-sample summary-level data. Genet Epidemiol 2023; 47(1): 3–25.

38.

Burgess

Butterworth

Thompson

. Mendelian randomization analysis with multiple genetic variants using summarized data. Genet Epidemiol 2013; 37(7): 658–665.

39.

Bowden

Davey Smith

Burgess

. Mendelian randomization with invalid instruments: effect estimation and bias detection through Egger regression. Int J Epidemiol 2015; 44(2): 512–525.

40.

Hartwig

Davey Smith

Bowden

. Robust inference in summary data Mendelian randomization via the zero modal pleiotropy assumption. Int J Epidemiol 2017; 46(6): 1985–1998.

41.

Gao

Zhang

Liu

. Rheumatoid arthritis: pathogenesis and therapeutic advances. MedComm 2024; 5(3): e509.

42.

Guo

Wang

, et al. Rheumatoid arthritis: pathological mechanisms and modern pharmacologic therapies. Bone Res 2018; 6(1): 15.

43.

Yerima

Akintayo

Adebajo

. Implications of a diagnosis of rheumatoid arthritis in resource-poor countries. Best Pract Res Clin Rheumatol 2022; 36(1): 101725.

44.

O’Neil

Alpízar-Rodríguez

Deane

. Rheumatoid arthritis: the continuum of disease and strategies for prediction, early intervention, and prevention. J Rheumatol 2024; 51(4): 337–349.

45.

Kushwaha

Kaushik

Kakkar

, et al. Red cell distribution width and neutrophil-lymphocyte ratio as inflammatory markers in patients with rheumatoid arthritis. Reumatologia 2023; 61(1): 13–20.

46.

Uslu

Küçük

Şahin

, et al. Two new inflammatory markers associated with Disease Activity Score-28 in patients with rheumatoid arthritis: neutrophil-lymphocyte ratio and platelet-lymphocyte ratio. Int J Rheum Dis 2015; 18(7): 731–735.

47.

Ali

Jabbar

Mohammed

. A comparative study of interleukin 6, inflammatory markers, ferritin, and hematological profile in rheumatoid arthritis patients with anemia of chronic disease and iron deficiency anemia. Anemia 2019; 2019: 3457347.

48.

Youssef

Hassan

Morad

, et al. Erythroferrone expression in anemic rheumatoid arthritis patients: is it disordered iron trafficking or disease activity? J Inflamm Res 2021; 14: 4445–4455.

49.

Xian

Yang

. Mendelian randomization and mediation analysis reveal the role of immune cell subsets in the causal pathways between blood cell perturbation responses and rheumatoid arthritis. Clin Rheumatol 2025; 44(4): 1537–1548.

50.

Cabrales

Tsai

Intaglietta

. Modulation of perfusion and oxygenation by red blood cell oxygen affinity during acute anemia. Am J Respir Cell Mol Biol 2008; 38(3): 354–361.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.09 MB

3.40 MB

0.08 MB

1.85 MB