This study used data from the electronic medical record system of Shanghai Pulmonary Hospital, a tertiary Grade A hospital. We enrolled patients with pulmonary sarcoidosis (stage I and II) and nonsarcoidosis diagnoses (including MLNT, lung cancer, pneumoconiosis, and lymphoma) hospitalized between May 2014 and May 2021, all exhibiting atypical imaging features upon their first admission.

The dataset was randomly split into a training set and a validation set in a 7:3 ratio. The parameters for the data preprocessing steps on the validation set (including missing value imputation and standardization) were all derived from the training set to prevent data leakage and ensure an unbiased evaluation of the validation results. Feature selection was performed on the training set. The standardized formula was:

where min(X) denoted the minimum value of the data and max(X) denoted the maximum value of the data. Model variable selection was performed solely on the training set, and the selected features were fixed before any validation set evaluation, following the process: first, univariate analysis was performed to identify variables with a significance level of P<.10. These variables were then subjected to importance analysis using the Boruta algorithm, which rigorously identifies features that are genuinely associated with the target variable by iteratively comparing the importance of original features with that of randomly permuted shadow features. This allows us to retain variables that are statistically relevant to sarcoidosis while filtering out noise, providing a reliable initial feature set [28]. The variables obtained from the Boruta algorithm were subjected to feature compression via least absolute shrinkage and selection operator (LASSO) regression (implemented via the glmnet package). LASSO applies L1 regularization to further shrink the coefficient vector. Building on the Boruta output, LASSO eliminates redundant or multicollinear variables, resulting in a more parsimonious model.

The selected feature set from the training set was used to develop models using 6 ML algorithms: logistic regression (LR), random forest (RF), support vector machine (SVM), extreme gradient boosting (XGBoost), naïve Bayes (NB), and k-nearest neighbors (KNNs). To counteract the class imbalance issue in the male data (sarcoidosis cases: 42/224, 18.7%), we used a model weighting approach when modeling the male subgroup. The weight calculation formula is:

where n is the total sample size, and n_i is the sample size for group i.

LR was implemented using the stats package, RF models were implemented using the randomForest package for the male cohort and the ranger package for the female cohort, SVM and NB using the e1071 package, XGBoost using the xgboost package, and KNN using the kknn package. In the male group modeling, weight parameters were incorporated into the code. Specifically, the core of the KNN algorithm lies in its distance-weighted classification: for each test sample, the KNNs are identified based on Euclidean distance, and different weights are assigned according to their distances—closer neighbors receive higher weights. The final predicted class was determined by weighted voting. For the NB algorithm, weighting was achieved by adjusting the prior probabilities. Under class imbalance, setting prior probabilities inversely proportional to the class proportions enhances the algorithm’s ability to recognize the minority class. The validation set was reserved for model evaluation. To achieve optimal model performance while minimizing the risk of overfitting, hyperparameter tuning was performed using nested cross-validation and grid search.

Each model was comprehensively evaluated for its discrimination, calibration, and clinical net benefit (using decision curve analysis [DCA]). The discriminatory ability of the 6 models was assessed using the AUC, sensitivity, specificity, positive predictive value, negative predictive value, F₁-score, and accuracy, while the DCA evaluated the practical utility and clinical application value of each model.

We performed 10-fold cross-validation on the training set using the best model to demonstrate the model’s stability during training. Then, we additionally assessed the robustness and generalizability of the model in the validation cohort by evaluating its discriminative performance between sarcoidosis and each individual disease separately.

SHAP is a method that ranks feature importance and interprets model outputs. It provides both global and local explanations for a model. Globally, it reveals overall feature-disease associations, while locally, it illustrates the reasoning behind individual predictions [17,23]. To facilitate clinical adoption and enhance understanding for clinicians, the final models were deployed as an interactive Shiny web application. By inputting a patient’s corresponding metric values, the application provides the probability score of the patient being diagnosed with pulmonary sarcoidosis and generates an individual force plot for case-specific interpretation.

Two independent data entry clerks performed double entry and proofreading to ensure accuracy and reliability. Except for lipid variables and the CD4⁺/CD8⁺ ratio, all other variables with missing values were imputed using the median method [29]. Specifically, the median calculated from the training set was used to impute missing values in the training set, and the same median value was then applied to the validation set to prevent data leakage. Categorical variables are presented as percentages and compared using the chi-square test. Normally distributed continuous data are expressed as mean (SD) and analyzed with the 2-sided t test, while nonnormally distributed data are summarized as median (IQR) and compared using the Mann-Whitney U test. The optimal cutoff value was established by maximizing the Youden index (sensitivity+specificity–1). Statistical significance was defined as a 2-sided P value <.05. All analyses, including ML model construction and SHAP value analysis, were performed using the R software (version 4.5.1; R Foundation for Statistical Computing). The R code used for the correlation analyses is provided in Multimedia Appendix 2 to facilitate reproducibility.

This study was conducted in accordance with the amended Declaration of Helsinki and the ethical guidelines of the institutional review board of Tongji University (project 221 approval K18-142). Clinical data were anonymized.

The entire dataset included 224 male cases and 196 female cases. After the data were divided into training and validation sets, the male and female groups consisted of 162 cases and 141 cases in the training set, and 62 cases and 55 cases in the validation set, respectively. Among male patients with pulmonary sarcoidosis, the median age was 44.0 (IQR 39.0-49.8) years; among female patients, the median age was 53.0 (IQR 48.0-60.0) years. In nonsarcoidosis cases, the median age was 59.0 (IQR 46.0-66.0) years for the male group and 55.5 (IQR 45.0-64.2) years for the female group. A higher proportion of bilateral hilar lymphadenopathy was observed in patients with pulmonary sarcoidosis compared to patients without sarcoidosis in both male and female groups. Table 1 describes the demographic information for male and female groups in the entire dataset, as well as the variables included in the models. Tables S2 and S3 in Multimedia Appendix 1 summarize the missing value patterns, while Tables S4 and S5 in Multimedia Appendix 1 present the statistical results for all variables.

Table 2 lists the detailed parameters for each model after tuning. Figure S1 in Multimedia Appendix 1 illustrates the variable selection process for the 2 models. In the female model, features corresponding to the minimum lambda value were selected, including lymphocyte count (Lym#), sACE, TC, triglyceride to high-density lipoprotein cholesterol ratio (TG/HDLC), glucose, CD3⁺, and CD8⁺. The features incorporated into the male model were age, Lym#, red cell distribution width standard deviation (RDW-SD), D-dimer, sACE, CD4⁺/CD8⁺, erythrocyte sedimentation rate (ESR), TG/HDLC, absolute basophil count (BASO#), and TP.

When evaluating the overall performance of the 6 ML models based on multiple performance metrics, including AUC, accuracy, specificity, F₁-score, and Brier score, SVM performed best among the male models (training set: AUC 0.989, 95% CI 0.975-1.000; validation set: AUC 0.896, 95% CI 0.746-1.000), while LR performed best among the female models (training set: AUC 0.848, 95% CI 0.781-0.915; validation set: AUC 0.762, 95% CI 0.632-0.892). Figure 2 presents the receiver operating characteristic curves, calibration curves, and DCA curves for all 6 ML algorithms in the validation set. The corresponding results for the training set are provided in Figure S2 in Multimedia Appendix 1. The discrimination metrics for each model are detailed in Table 3 and Table S6 in Multimedia Appendix 1.

Figure S3 in Multimedia Appendix 1 presents the detailed results of the 10-fold cross-validation. The selected models, SVM for the male group and LR for the female group, achieved a 10-fold cross-validation AUC of 0.922 (SD 0.079) and 0.825 (SD 0.110) within the training sets, respectively. Table S7 in Multimedia Appendix 1 shows the cross-validation results as well as the discriminative performance for sarcoidosis versus each specific comparator disease.

The SHAP method explains the final model’s output by calculating the contribution of each variable to the prediction, aiding clinicians in understanding and trusting the diagnostic model. The global interpretation, depicted in the SHAP summary plot (Figure 3B), illustrates the collective impact and importance of all features. The mean SHAP value was used to assess each feature’s overall contribution to the model. In the male model, Lym# and TP were the most important features. In the female model, the most significant predictor was sACE, followed by the TG/HDLC (Figure S4 in Multimedia Appendix 1).

SHAP dependence plots help visualize how a single feature influences the model’s output across its range of values, showing the relationship between the feature value and its SHAP value (Figure 3A). A SHAP value greater than 0 corresponds to the positive class in the model’s prediction, which, in this study, indicates the pulmonary sarcoidosis class. Figure 3D is a beeswarm plot. Each point represents one sample, and the x axis indicates the SHAP value of the corresponding feature for that sample. A positive SHAP value means the feature increases the model’s output (eg, the predicted probability of sarcoidosis), while a negative SHAP value means it decreases the output. The color of each point represents the original value of that feature in the sample. By examining the distribution of points and the color gradient, one can infer the relationship between the feature value and its SHAP contribution. For example, if high feature values cluster in the positive SHAP region, higher values of that feature tend to push the prediction toward the positive class; conversely, if high values cluster in the negative SHAP region, higher values are associated with a negative prediction. Force plots (or waterfall plots) detail the prediction for an individual sample, demonstrating how each feature value contributes to pushing the model’s base value toward the final prediction output for that specific case (Figure 3C and E).

The final predictive model was implemented in a web-based application to enhance its utility in clinical practice, as illustrated in Figure 4 and Figure S5 in Multimedia Appendix 1. Upon entering the required feature values, the application automatically generates a diagnostic output. Additionally, an individualized force plot is displayed for each patient, indicating which features contributed to the sarcoidosis prediction: features on the right in red push the diagnosis toward “nonsarcoidosis,” while those on the left in yellow support a “pulmonary sarcoidosis” classification. The web applications for the male and female models are publicly accessible for research and demonstration purposes [30,31].

This study developed diagnostic models to differentiate early-stage pulmonary sarcoidosis with normal sACE activities from cases of other atypical pulmonary diseases. To our knowledge, this is the first ML-based diagnostic model for pulmonary sarcoidosis. Our study indicated that, although sarcoidosis itself is a highly heterogeneous disease, sex-stratified ML models may provide useful discrimination between pulmonary sarcoidosis and other similar diseases, with relatively better performance observed in the male group.

The model for male patients incorporated the following indices: age, Lym#, RDW-SD, sACE, CD4⁺/CD8⁺, ESR, TG/HDLC, BASO#, TP, and D-dimer. The model for female patients used Lym#, sACE, TC, TG/HDLC, glucose, CD3⁺%, and CD8⁺%. We developed a web-based tool using the selected ML algorithms—the SVM model for male group (training set: AUC 0.989, 95% CI 0.975-1.000; validation set: AUC 0.896, 95% CI 0.746-1.000) and the LR model for female group (training set: AUC 0.848, 95% CI 0.781-0.915; validation set: AUC 0.762, 95% CI 0.632-0.892)—to facilitate the clinical application of this early pulmonary sarcoidosis differential diagnosis model. Additionally, we visualized the model using the SHAP method to enhance clinicians’ understanding.

Our model showed promising performance across internal cross-validation and analyses of individual diseases, although further validation in independent cohorts is needed to confirm its generalizability. This study fully considered sex differences in sarcoidosis and adopted a sex-stratified modeling strategy. The National Institutes of Health policy on Sex as a Biological Variable, issued in 2015 [32], calls for the inclusion of sex as a biological factor in biomedical research. Multiple studies have shown that there are significant immunological, clinical, and genetic differences between male and female patients with sarcoidosis [33-35]. A study on the pulmonary functional phenotype of sarcoidosis found that there are differences in pulmonary functional phenotypes between the sexes [36]. In addition, the association between certain human leukocyte antigen markers and sarcoidosis also shows sex differences [34], suggesting that immune pathways play different roles in the male and female groups. In this study, the selected model was SVM for the male group and LR for the female group, a difference that can be attributed to distinct intrinsic data characteristics between the 2 groups. The indices included in the male model may exhibit nonlinear relationships or threshold effects with pulmonary sarcoidosis. SVM can effectively capture complex nonlinear relationships among variables without requiring predefined functional assumptions. They are particularly suitable for small-to-moderate sample sizes and high-dimensional datasets and have demonstrated strong robustness and favorable generalization performance in various classification tasks. The female group had a smaller sample size and fewer features; LR is suitable for small-sample data and is less prone to overfitting, and its low-variance property makes it more robust. Moreover, the sets of included indices differ slightly between sexes (eg, D-dimer, RDW-SD, and BASO# in the male group; TC, glucose, and CD3⁺% in the female group), suggesting that sex-specific modeling helps remove confounding bias and improve predictive performance.

To address the class imbalance in the male data, we applied sample weighting during model development, which effectively increased the influence of the minority class (sarcoidosis). As a result, the model showed preliminary discriminative ability on the validation set (male SVM: AUC=0.896; sensitivity=0.700; specificity=0.942). The male model showed relatively better performance across several evaluation metrics than the female model, although the differences should be interpreted cautiously, given the limited sample size and wide CIs. Studies have suggested that estrogen exposure may contribute to the delayed onset of sarcoidosis in the female group by influencing the Th1/Th2 balance [37,38]. In this study, estrogen levels in the female patients fluctuate considerably with the menstrual cycle and age. This may result in greater heterogeneity or atypical patterns in immune-related biomarkers among female patients, leading to a more ambiguous and noisier relationship between features and outcomes, thereby reducing the predictive power of the model. Additionally, the male group in this study had a slightly larger sample size (n=224 cases) than the female group (n=196 cases), and the male model included more features (n=10) than the female model (n=7). These factors may also partly explain the more stable performance of the male model.

It should be noted that although the RF algorithm uses bootstrap sampling and ensemble learning (bagging) to mitigate overfitting, the AUC on the training set in this study still reached 1.000. This suggests that the current hyperparameter tuning range (eg, tree depth, minimum node size, and regularization parameters) may still be insufficient to adequately penalize model complexity, or alternatively, the combination of a small sample size and a high-dimensional feature space may have allowed the model to easily memorize noise in the training data. Specifically, while bootstrap sampling reduces variance, when there are many features or strong discriminative predictors in the data, individual trees can still achieve near-perfect splits on the training samples; after aggregating multiple trees, the predictions on the training set become highly consistent, resulting in an AUC of 1. The relatively limited sample size in this study further exacerbates this issue. Future research should adopt stricter regularization strategies (eg, limiting tree depth, reducing the number of features, and increasing pruning strength) alongside larger sample sizes to enhance model generalizability.

All indices incorporated in our model are routinely available in clinical practice, including complete blood count, inflammatory markers, T-cell subsets, biochemical tests, and lipid profiles, offering good operability and generalizability. Moreover, these indices collectively reflect multiple pathophysiological dimensions—immune activation, chronic inflammation, metabolic disturbance, and nutritional status—closely aligning with the complex mechanisms of pulmonary nodular granulomatous disease. Traditional approaches relying on a single index or manually exploring combinations of indices have often been inefficient and difficult to translate into clinical use. To address this, we used the Boruta algorithm combined with LASSO regression for feature selection, efficiently identifying the above representative and jointly effective indices.

Both models included sACE and T-cell subset indices, which are the most common biomarkers for sarcoidosis. In this research, despite sACE levels of all the cases being within the normal ranges, sACE ranked fourth in importance in the male model but first in the female model, possibly due to the higher prevalence of pneumoconiosis—which also elevates sACE—among the male group [39]. However, sACE activity has been found to be associated with population characteristics. In certain populations, sACE levels fall within the high end of the normal range, which is associated with the angiotensin-converting enzyme gene insertion/deletion polymorphism (ACE I/D) [40]. Research by Csongrádi et al [41] found that combining sACE activity testing with ACE I/D genotyping improved the specificity of sarcoidosis diagnosis to 100%, albeit with a relatively low sensitivity of 42.5%. Both male and female groups showed a reduction in Lym# in peripheral blood, which may be associated with more impaired cellular immune function, for example, dysregulation of T cell subsets, in sarcoidosis than in other differentiated diseases [42]. Elevated CD4⁺/CD8⁺ ratios in bronchoalveolar lavage fluid have been demonstrated in multiple studies to hold diagnostic value for sarcoidosis [43]. Nevertheless, the invasive risk for bronchoalveolar lavage fluid sampling limits its clinical applicability. Our study demonstrated that the peripheral blood CD4⁺/CD8⁺ ratio in male, and percentages of CD8⁺ and CD3⁺ in female, when combined with other readily available indicators, can yield a well-performing diagnostic model for sarcoidosis.

Patients with sarcoidosis exhibit dyslipidemia. An integrative observational and Mendelian randomization study on blood lipids and sarcoidosis confirmed that high-density lipoprotein cholesterol and TC are closely associated with the disease [26], which may be attributed to abnormal lipid metabolism in macrophages during sarcoidosis pathogenesis [44]. Consistently, in this study, 2 indices reflecting lipid metabolic status [45], TC and TG/HDLC ratio, were selected. Currently, the TG/HDLC ratio has been extensively studied in metabolic and cardiovascular diseases [46,47]; its relationship pattern with sarcoidosis warrants further investigation. Our model, leveraging easily obtainable samples and providing rapid results, offers a promising solution for the early differential diagnosis of pulmonary sarcoidosis. In addition, the web-based applications developed in this study are research prototypes intended solely for demonstration and peer review purposes. Translation into clinical practice would require implementation of the following security measures: (1) deploy the application within the hospital’s internal network and enforce user authentication or authorization via Lightweight Directory Access Protocol or Shiny Proxy; (2) use HTTPS to encrypt all data transmission; (3) do not store any patient identifiable information; input data are used only for the current prediction and are not persisted; (4) maintain audit logs (timestamp, user, and deidentified input values) for traceability; (5) establish model version control and an approval workflow for updates; and (6) include a clear disclaimer that the tool is for decision support only, and the final diagnosis must be made by a clinician.

This study has several limitations. First, the single-center nature of our study represents a key limitation, as it introduces potential selection bias and may restrict the model’s generalizability to other demographic and clinical environments. Although internal cross-validation was robust, this constraint necessitates future external validation. Multicenter collaborations are needed to test, refine, and ultimately ensure the broad clinical applicability of the model. Second, although we split an independent validation set and performed robustness analyses, the inherent limitation of sample size remains, which reduces the reliability of the DCA and calibration curves. Future studies with larger sample sizes and external validation are needed to further confirm the robustness and generalizability of the models. Third, class imbalance exists in the male data. Although we applied sample weighting during modeling and acceptable predictive capacity on the validation set, increasing the male sample size and using other methods to address class imbalance are still necessary in the future to enhance model performance and generalizability. Fourth, the cohort splitting strategy was performed using random sampling rather than stratified sampling. Although the distributions of disease categories were comparable between the training and validation cohorts in the present study, future multicenter studies with larger cohorts should consider stratified sampling based on disease categories and external validation to further assess the robustness, clinical applicability, and generalizability of the proposed model. Fifth, although both models demonstrated acceptable discriminative ability, the CIs for sensitivity remained relatively wide. This is partly due to the small sample size and partly because some features of sarcoidosis (eg, elevated sACE and decreased lymphocytes) may also be present in other diseases. Future studies should adopt larger sample sizes and explore novel blood biomarkers. Sixth, the data incorporated into our model consisted of peripheral blood biomarkers, which are subject to inherent variability. However, all blood samples were collected as fasting morning draws, and we excluded patients with a history of medication use (eg, corticosteroids) prior to admission to minimize confounding factors. Seventh, we included patients who had used cephalosporins within 2 weeks, because these drugs are commonly used as empirical antibiotics in clinical practice. To preserve sample size and reflect real-world clinical scenarios, we retained such cases; moreover, short-term use of cephalosporins has a minimal impact on blood parameters [48] (eg, the half-life of ceftriaxone is approximately 1.5-2 hours) [49]. Nevertheless, we acknowledge that incomplete control of antibiotic exposure remains a potential confounder. We suggest that future studies with stricter inclusion criteria or prospective designs may further address this issue. Eighth, although we performed robustness and generalizability analyses for sarcoidosis versus each single comparator disease, the sample sizes for some diseases were too small to allow meaningful statistical analysis. Therefore, the model’s generalizability to rare diseases such as lymphoma remains unclear. Future multicenter studies with larger sample sizes are needed to comprehensively validate the model’s generalizability across all comparator diseases. Ninth, while the webpage is flexible and user-friendly, it must undergo prospective validation and obtain ethics approval before being implemented in clinical practice.

This study developed an interpretable model for the differential diagnosis of early pulmonary sarcoidosis versus other radiologically atypical pulmonary diseases with normal sACE levels using conventional peripheral blood markers. Given the wide CIs reflecting limited sample size, these findings are preliminary and require validation in larger independent cohorts. If successfully validated, the model holds promise to support clinical decision-making, facilitate the management of early-stage pulmonary sarcoidosis, and provide a scientific basis for clinicians when considering invasive diagnostic examinations.