A multicenter, retrospective, cross-sectional study was conducted, utilizing clinical data collected from three independent tertiary medical centers. Patients were consecutively enrolled between December 2012 and January 2025. All participants underwent systematic evaluation by rheumatologists, and definitive diagnoses were established based on integrated assessments of clinical manifestations, laboratory findings, and imaging examinations, in accordance with internationally recognized or widely accepted diagnostic criteria.

According to disease status, patients were classified into three groups. The first group comprised patients with AS alone. The diagnosis of AS was based on the ASAS (Assessment of SpondyloArthritis International Society) classification criteria and was independently confirmed by at least two experienced rheumatologists. Diagnostic evaluation incorporated typical clinical features, such as chronic inflammatory back pain and limitation of spinal or sacroiliac joint mobility, in line with the ASAS clinical criteria. Imaging evidence of sacroiliitis or structural changes on radiography, computed tomography (CT), or magnetic resonance imaging (MRI), consistent with the ASAS imaging criteria, were also considered. Additionally, relevant laboratory findings, including HLA-B27 testing, were evaluated as part of the diagnostic process. Patients in this group had no concomitant autoimmune diseases (11).

The second group included patients with autoimmune diseases other than AS, without clinical or imaging evidence of AS. These autoimmune conditions included psoriasis, inflammatory bowel disease, uveitis, dermatitis, rheumatoid arthritis, systemic lupus erythematosus, Sjögren’s syndrome, Crohn’s disease, amyopathic dermatomyositis, vasculitis, Behçet’s disease, systemic sclerosis, and other immune-mediated disorders. All diagnoses were made according to the corresponding internationally recognized or nationally accepted diagnostic criteria, including the CASPAR criteria for psoriatic arthritis, the 2010 ACR/EULAR criteria for rheumatoid arthritis, the 2019 EULAR/ACR criteria for systemic lupus erythematosus, the 2016 ACR/EULAR criteria for Sjögren’s syndrome, the 2012 Chapel Hill criteria for vasculitis, and the 1990 International Criteria for Behçet’s disease, among others. These criteria integrate clinical features, laboratory tests, imaging studies, and specific antibody markers to ensure accurate diagnosis and consistency.

The third group consisted of patients with AS comorbid with other autoimmune diseases. These patients fulfilled the diagnostic criteria for AS as well as for at least one additional autoimmune disease, with each diagnosis independently confirmed by specialist physicians.

Clinical and laboratory data were collected at the time of initial evaluation or hospital admission. Demographic variables included age and gender. Systemic inflammatory indices comprised the systemic immune-inflammation index (SII), systemic inflammation response index (SIRI), neutrophil-to-lymphocyte ratio (NLR), platelet-to-lymphocyte ratio (PLR), platelet-to-neutrophil ratio (PNR), and neutrophil-to-monocyte ratio (NMR).

Routine hematological parameters included white blood cell count (WBC), red blood cell count (RBC), hemoglobin (Hb), hematocrit (HCT), mean corpuscular volume (MCV), mean corpuscular hemoglobin (MCH), mean corpuscular hemoglobin concentration (MCHC), platelet count (PLT), mean platelet volume (MPV), platelet distribution width (PDW), and plateletcrit (PCT). Differential leukocyte counts included neutrophils (NE), lymphocytes (LY), monocytes (MO), eosinophils (EO), and basophils (BA).

Inflammatory activity was assessed using the erythrocyte sedimentation rate (ESR).

Liver function parameters included total bilirubin (TBIL), direct bilirubin (DBIL), indirect bilirubin (IBIL), the DBIL/TBIL ratio, total protein (TP), albumin (ALB), globulin (GLOB), the albumin-to-globulin ratio (ALB/GLOB), gamma-glutamyl transferase (GGT), total bile acid (TBA), aspartate aminotransferase (AST), alanine aminotransferase (ALT), the AST/ALT ratio, alkaline phosphatase (ALP), prealbumin (PREALB), and cholinesterase (CHE).

Renal function and related metabolic parameters included blood urea nitrogen (BUN), creatinine (Cr), uric acid (UA), bicarbonate (HCO₃), estimated glomerular filtration rate (eGFR), and cystatin C (CysC).

Prior to model development, clinical data were preprocessed to ensure data quality. Patient records with missing values were excluded, and only patients with complete clinical and laboratory data were included in the analysis. The data were then organized and analyzed separately according to their originating centers.

Continuous variables were uniformly processed before analysis to minimize the influence of different measurement scales on model training, while categorical variables were encoded into numerical formats suitable for machine learning algorithms. The same preprocessing procedures were applied across all centers to ensure comparability and robustness of the analytical results (12, 13).

After data screening, the final study population comprised three disease groups across the three centers. In the AS-alone group, 326 patients were included from Center 1, 108 from Center 2, and 110 from Center 3. In the group with autoimmune diseases other than AS, 603 patients were enrolled from Center 1, 201 from Center 2, and 201 from Center 3. In the AS with autoimmune comorbidity group, 343 patients were included from Center 1, 114 from Center 2, and 115 from Center 3.

After data preprocessing, multiple machine learning–based binary classification models were developed using clinical and laboratory variables to discriminate between different disease states. The algorithms applied included Naïve Bayes, discriminant analysis, logistic regression, LASSO, partial least squares, k-nearest neighbors, support vector machines, artificial neural networks, random forest, gradient boosting and other boosting methods, as well as ensemble learning models such as LightGBM and CatBoost.

Model training was primarily performed on the training dataset with repeated five-fold cross-validation to assess model stability. Model performance was evaluated in both the training and testing datasets using the area under the receiver operating characteristic curve (AUC), sensitivity, specificity, accuracy, F1 score, and Youden index. Receiver operating characteristic (ROC) curves and decision curve analysis (DCA) were used to compare discriminative ability and potential clinical utility across models (14, 15). The model with the highest AUC in the testing dataset was selected for subsequent analyses.

To enhance model interpretability, SHAP were applied to the optimal machine learning model selected based on test-set performance. SHAP was used to quantify the contribution of individual clinical and laboratory features to model predictions at both global and individual levels, thereby elucidating the model’s decision-making process (8).

For tree-based models (e.g., LightGBM and CatBoost), SHAP values were computed using model-specific built-in methods, while a model-agnostic kernel-based SHAP approach was used for the remaining models. Feature importance was assessed using mean absolute SHAP values, and SHAP summary and dependence plots were generated to visualize the direction and magnitude of feature effects (16). In addition, waterfall or force plots were used to provide intuitive explanations for representative individual predictions.

In this study, we performed feature engineering and feature selection through the following steps: (1) Data Preprocessing: Missing values were handled, and continuous variables were standardized to mitigate the effects of different scales. Categorical variables (e.g., gender) were converted into numerical format. (2) Feature Selection Method: A total of 48 variables were initially considered. Machine learning algorithms were applied to assess the importance of features in clinical and laboratory data. The most important features for distinguishing AS from autoimmune comorbidities were selected. (3) SHAP Analysis: To enhance model interpretability, SHAP (Shapley Additive Explanations) was used to quantify the contribution of each feature to the model’s predictions. This analysis helped confirm the clinical and biological relevance of the selected features. (4) Feature Selection Criteria: The selected features demonstrated strong statistical correlation and were aligned with clinical relevance. Feature selection ensured that the model was both effective in prediction and resistant to overfitting.

After merging data from the three centers, a total of 2,121 patients were included in the analysis, comprising 544 patients with AS alone, 1,005 patients with autoimmune diseases without AS, and 572 patients with AS combined with other autoimmune diseases. Comparisons of demographic characteristics, systemic inflammatory indices, hematological parameters, liver function, and renal function–related variables among the three groups are summarized in Tables 1, 2, 3, 4. The complete research flowchart is shown in Figure 1.

Significant differences in demographic characteristics were observed among the three groups (Table 1). Patients in the AS-alone group were younger and predominantly male, whereas a higher proportion of female patients and older age were observed in the autoimmune disease group and the AS comorbidity group (all P < 0.01).

Systemic inflammatory indices differed significantly across groups (Table 1). Compared with the autoimmune disease–alone group, patients with AS alone exhibited higher levels of ESR, SII, SIRI, NLR, and PLR. Patients with AS and autoimmune comorbidities demonstrated distinct inflammatory profiles, with several indices differing from both the AS-alone and autoimmune disease–alone groups.

Hematological parameters also showed significant intergroup differences (Table 2). WBC, RBC, Hb, HCT, PLT, and related indices differed significantly among the three groups (most P < 0.05). In addition, differential leukocyte counts, including NE, LY, and MO, exhibited distinct distributions across disease groups, with the AS comorbidity group showing patterns that differed from both the AS-alone and autoimmune disease–alone groups.

Liver function parameters differed significantly among the three groups (Table 3). TBIL, DBIL, TP, ALB, GLOB, ALB/GLOB, AST, ALT, and TBA showed statistically significant differences across groups (most P < 0.05). Compared with the AS-alone group, patients with AS and autoimmune comorbidities exhibited more pronounced alterations in several liver-related indices.

Renal function and metabolic parameters also varied significantly among the three groups (Table 4). BUN, Cr, UA, HCO₃, eGFR, and CysC differed significantly across groups (all P < 0.01). The AS comorbidity group displayed renal and metabolic profiles that were not entirely consistent with those of the AS-alone or autoimmune disease–alone groups.

To evaluate inter-center consistency, demographic and laboratory variables were compared across the three centers within each disease group, including AS alone, autoimmune diseases alone, and AS with autoimmune comorbidities (Supplementary Tables 1-3). No significant inter-center differences were observed within the same disease category, indicating minimal center-related heterogeneity. Accordingly, data from the three centers were pooled for subsequent analyses.

Based on the merged multicenter dataset, multiple machine-learning–based binary classification models were constructed to discriminate patients with AS alone from those with AS combined with autoimmune diseases. Model performance, potential clinical utility, and interpretability were systematically evaluated (Figure 2).

Comparisons of overall model performance demonstrated substantial differences among algorithms in terms of accuracy, sensitivity, specificity, F1 score, and AUC (Figures 2A1, A2). Overall, ensemble learning models—including LightGBM, random forest, gradient boosting, and CatBoost—outperformed traditional linear models and several single classifiers in both the training and validation datasets. Among these, LightGBM consistently exhibited superior and more stable performance across multiple evaluation metrics.

ROC curve analyses are shown in Figures 2B1, B2. LightGBM achieved ROC curves that were consistently above those of other models in both the training and validation sets, with AUC values ranking in the top tier, indicating strong discriminative ability. AUC forest plots further confirmed that LightGBM demonstrated favorable AUC values with relatively narrow confidence intervals, suggesting robust and stable model performance (Figures 2C1, C2).

To evaluate potential clinical utility, DCA was performed for all models (Figures 2D1, D2). Across a wide range of threshold probabilities, LightGBM provided greater net benefit than other models, supporting its potential usefulness in distinguishing AS alone from AS with autoimmune comorbidities.

Confusion matrices are presented in Figures 2E1, E2. LightGBM showed excellent fitting performance in the training set and maintained good classification performance in the test set, indicating acceptable generalizability. Clinical impact curves (Figures 2F1, F2) further demonstrated good agreement between the number of patients predicted to be at high risk and the observed outcomes across different threshold probabilities.

Calibration curves for the training and validation datasets are shown in Figures 2G1, G2. The curves were generally close to the ideal reference line, indicating good calibration of predicted probabilities generated by the LightGBM model.

After selecting LightGBM as the optimal model, SHAP were applied to interpret model predictions (Figures 2H–J). Global feature importance analysis revealed that age, Gender, eGFR, DBIL, CysC, BUN, UA, and ALB were the primary contributors to discriminating AS alone from AS with autoimmune comorbidities (Figure 2H).

Integration of SHAP summary plots with the statistical comparisons in Tables 1-4 demonstrated good consistency between feature expression patterns and model-predicted directions. Specifically, older age and female sex were associated with higher SHAP values and tended to drive predictions toward AS with autoimmune comorbidities, with female patients with AS more frequently identified by the model as having additional autoimmune diseases. Regarding renal function–related variables, decreased eGFR and increased levels of CysC, BUN, and UA showed positive contributions to predicting AS with autoimmune comorbidities, whereas patients with AS alone generally exhibited relatively preserved renal function. In addition, elevated DBIL and reduced ALB levels also contributed positively to the prediction of AS with autoimmune comorbidities.

SHAP dependence plots further indicated nonlinear relationships between several key variables and model output, suggesting that the model captured complex feature patterns that are not readily reflected by conventional statistical comparisons (Figure 2I). Moreover, individual-level SHAP explanations for a representative patient are illustrated using waterfall and force plots (Figures 2J1, J2) to demonstrate the cumulative contribution of individual features to the final prediction. Individual SHAP explanations for additional patients are provided in Supplementary Figure 2.

Based on the merged multicenter dataset, multiple machine learning–based binary classification models were further developed to discriminate patients with autoimmune diseases alone from those with AS combined with autoimmune diseases. Model discriminative performance, potential clinical utility, and interpretability were systematically evaluated (Figure 3).

Comparisons of overall model performance demonstrated marked differences among algorithms in accuracy, sensitivity, specificity, F1 score, and AUC (Figures 3A1, A2). Overall, ensemble learning models outperformed traditional linear models in both the training and validation datasets, among which LightGBM exhibited relatively superior and stable predictive performance across multiple evaluation metrics.

ROC curve analyses are presented in Figures 3B1, B2. LightGBM achieved ROC curves that were consistently above those of other models in both the training and validation sets, with relatively high AUC values, indicating strong discriminative ability for distinguishing autoimmune diseases alone from AS with autoimmune comorbidities. AUC forest plots further showed that LightGBM demonstrated favorable AUC values with relatively narrow confidence intervals, suggesting robust model performance (Figures 3C1, C2).

To assess potential clinical applicability, decision curve analysis (DCA) was performed for all models (Figures 3D1, D2). Across a broad range of threshold probabilities, LightGBM yielded greater net benefit than other models, supporting its potential utility in identifying patients with AS combined with autoimmune diseases.

Confusion matrices are shown in Figures 3E1, E2. LightGBM displayed strong fitting performance in the training set and maintained satisfactory classification performance in the test set, indicating acceptable generalizability. Clinical impact curves (Figures 3F1, F2) further demonstrated good agreement between the number of patients predicted to be at high risk and the observed outcomes across different risk thresholds.

Calibration curves for the training and validation datasets are presented in Figures 3G1, G2. The curves were generally close to the ideal reference line, indicating good calibration of predicted probabilities generated by the LightGBM model.

After selecting LightGBM as the optimal model, SHapley Additive exPlanations (SHAP) were applied to interpret model predictions (Figures 3H–J). Global feature importance analysis revealed that eGFR, GLOB, ALP, Cr, HCO₃^-, CysC, PREALB, CHE, ESR, and UA were among the most influential variables for discriminating autoimmune diseases alone from AS with autoimmune comorbidities (Figure 3H).

Integration of SHAP summary plots with the statistical comparisons in Tables 1, 2, 3, 4 demonstrated good consistency between feature distributions and model-predicted directions. Specifically, compared with patients with autoimmune diseases alone, those with AS and autoimmune comorbidities exhibited more pronounced abnormalities in renal function–related parameters, including decreased eGFR and increased levels of Cr, CysC, and UA, which showed positive contributions to predicting AS with autoimmune comorbidities in the SHAP analysis (Table 4). In addition, protein metabolism and liver function–related variables, such as GLOB, PREALB, CHE, and ALP, also showed substantial contributions in SHAP analysis, consistent with the intergroup differences observed in Table 3.

Regarding inflammatory markers, ESR contributed to model discrimination and was generally higher in patients with AS and autoimmune comorbidities, consistent with the statistical differences shown in Table 1. SHAP dependence plots further indicated nonlinear relationships between several key variables and model output, suggesting that the model captured complex feature patterns not fully reflected by conventional statistical comparisons (Figure 3I). Moreover, individual-level SHAP explanations for a representative patient are illustrated using waterfall and force plots (Figures 3J1, 2) to demonstrate the cumulative contribution of individual features to the final prediction. Individual SHAP explanations for additional patients are provided in Supplementary Figure.