Introduction
Kawasaki disease (KD) is an acute and immune-mediated vasculitis, principally influencing children under 5 years old. The coronary artery lesions (CALs) are the most common and detrimental complication, resulting in KD being emerged as the major cause of acquired heart disease in developed countries [1]. The effective treatment with intravenous immunoglobulin (IVIG) and aspirin remarkably decreases the incidence of CALs from 25% to 4% [2], while there are approximately 10%–20% of rate of IVIG-resistance after initial therapy. Remarkably, IVIG-resistance is widespread regarded as the risk factor for CALs. Therefore, it is crucial to construct a reliable model to predict IVIG-resistance and then conduct the suitable regimes as early as possible in order to decrease the occurrence of CALs.
In the past decade, numerous investigators worldwide have been exploring reliable predictive models or single indicators of-IVIG-resistance. Examples of such single predictors include peripheral blood cell parameters, C-reactive protein (CRP), procalcitonin, erythrocyte sedimentation rate (ESR), CRP/albumin (ALB) ratio [3], prognostic value of pretreatment prognostic nutritional index [4], serological indicators [5,6] and others. Additionally, various predictive models based on Japanese population, such as Egami [7], Kobayashi [8], and Sano [9] have been extensively studied in other countries, demonstrating low sensitivity or specificity in other countries [10-13]. Obviously,the overall applicability of these models was limited in the general population.
Machine learning (ML) is an innovative field in medicine that holds promise for the future of biomedical research. ML has demonstrated its effectiveness in various diseases, including cancer [14], acute kidney injury [15] and malnutrition [16]. Some studies have shown that ML can accurately predict IVIG-resistance [17,18]. Extreme gradient boosting (XGBoost), logistic, support vector machine (SVM), light gradient boosting machine (LightGBM) and multiple layers perception (MLP) cloud enhance the optimization of the mathematical relationship between learning covariates and outcomes [19]. Furthermore, ML-based models demonstrate superior performance in managing nonlinear data relationships, positioning them as an effective surrogate for predicting outcomes [20]. However, the black box characteristic of ML algorithm may limit the trust of patients and clinicians, thereby impeding their widespread application in clinical settings. SHapley Additive exPlanations (SHAP) is a method utilized for interpreting individual predictions by determining the optimal Shapley value from game theory and TreeSHAP has been extensively utilized to facilitate the interpretation of ML algorithms, thereby overcoming the limitations of black box models. It can provide insights into the complex relationship between features and predictions, offering intuitive explanations, with significant advancements observed in the outcomes analysis of diseases [21]. This can be particularly helpful for clinicians in understanding how the model accurately predicts outcomes and improving early interventions.
Therefore, the objectives of this study were threefold: first, to identify the associated risk factors of IVIG-resistance in the Mian Yang region; second, to develop and validate a high-performing ML model; and third, to interpret the ML model using SHAP.
Methods
1. Patients
The clinical data of eligible patients admitted to Mian yang Central Hospital between December 2014 and November 2022 were screened. The diagnostic criteria for KD followed the American Heart Association 2004 guidelines. IVIG-resistance was defined as the presence of persistent or recurrent fever lasting for at least 36 hours after the initial IVIG infusion. All patients received standard treatment, which included IVIG infusion (2g/kg) and oral aspirin (30–50 mg/kg), within 10 days of disease onset. The exclusion criteria were as follows: (1) patients who did not receive IVIG therapy during hospitalization, (2) patients who received glucocorticoids or other immunosuppressive drugs during or before the initial IVIG therapy, (3) patients who were rehospitalized due to recurrence of KD, and (4) patients with a significant amount of missing clinical data (exceeding 20% of the total). During the acute fever period before IVIG treatment, various clinical data were collected, including white blood cell count, neutrophil count, lymphocyte (LYM) count, eosinophil (EOS) count, monocyte count, basophil count, platelet (PLT) count, hemoglobin (Hb) level, ESR, CRP level, gamma-glutamyl transpeptidase level, lactate dehydrogenase level, direct bilirubin (DBIL) level, total bilirubin (TBIL) level, globulin (GLO) level, ALB level, aspartate transaminase (AST)level, and alanine aminotransferase level. Throughout the entire study process, including clinical information collection, patients' identity information was filtered to ensure patients' privacy. This study was approved by the Ethics Committee on Human Subjects at Mianyang Central Hospital (No. KY2024062), in accordance with the 1975 Helsinki Declaration revised in 1996. Written informed consent was waived due to the retrospective nature and anonymity of the data.
2. Multiple ML algorithms
Each algorithm exhibits unique characteristics and capabilities within specific contexts. SVM is a versatile tool for linear/nonlinear classification tasks and is particularly useful for ML challenges with limited samples, albeit its performance is contingent upon parameter tuning and function selection [22]. XGBoost, an ensemble learning algorithm, is adept at efficiently handling missing data to generate precise prediction models [23]. The MLP is capable of learning nonlinear models and processing high-dimensional data, however, it necessitates rigorous standardization and preprocessing of the data, and is associated with a relatively lengthy training period [24]. LightGBM demonstrates exceptional performance in processing extensive structured datasets at a rapid training pace, although it may be influenced by the number of features and sample size [25,26].
3. Model interpretation
SHAP methodology was utilized to elucidate the optimal predictive model. Feature prioritization was achieved through the calculation of SHAP values, with features ranked based on the mean absolute value of their respective SHAP values. The integration of ML with SHAP has the potential to offer a clear rationale for the prediction model.
4. Statistical analysis
R ver. 4.2.3 (R Foundation for Statistical Computing, Vienna, Austria) and Python ver. 1.11.3 (Python Software Foundation) were chosen for statistical analysis. For ML models, Python (XGBoost, LightGBM, and scikit-learn) packages was used to construct models. Python ver. 0.43.0 was employed for SHAP analysis. Missing data (less than 20% of the total) was filled with random forest algorithms. Continuous data are expressed as median with interquartile range. Mann-Whitney U test and chi-square test were performed for continuous and categorical variables, respectively. The multivariable logistic regression analysis with backward stepwise selection was performed on the significant factors in univariate logistic analysis. A probability value of P<0.05 was required for entry into the model and P>0.05 for elimination. Relative risks were evaluated with odds ratios and 95% confidence intervals (CIs). In order to evaluate the dose association of independent risk factors with IVIG-resistance, a restricted cubic spline (RCS) was employed.
The whole data was randomly broken down into training and testing set based on the ratio of 8:2 in order to select the best prediction model. This random division was recured until the patient data were evenly distributed between the 2 groups. Five types of ML algorithm, including logistic, XGBoost, LightGBM, SVM, and MLP were developed to construct predictive model. Under the receiver operating characteristics curve, precision recall (PR) curves, area under PR (AP) curve and decision curve analysis (DCA) were employed to assess the predictive performance of each model in order to select the best-performing model. In addition, the accuracy, F1-score value, Kappa coefficient, sensitivity and specificity of each model also were used to further assess the models. Best hyperparameter and 10 cross-validation were used to avoid the risk of overfitting and heighten model accuracy. The selected optimal model was trained, verified, and tested using a 10-fold cross-validation technique. Additionally, SHAP analysis was conducted to determine the significance of individual features in the selected ML model. A 2-sided P value of less than 0.05 was considered statistically significant
Results
1. Baseline characteristics
Table 1 presents the demographic and laboratory parameters of the study participants. A total of 602 children with KD were recruited, out of which 68 (11.3%) experienced IVIG-resistance. The IVIG-resistant group exhibited higher levels of CRP, DBIL, TBIL, GLO, and AST, while showing lower counts of LYM, EOS, PLT, and lower levels of Hb, Na+, and ALB (all P<0.05) when compared to the IVIG-responsive group.
2. Association between clinical indicators and IVIG-resistance
In univariate logistic analysis, lower LYM, EOS, Na+, Hb and ALB, as well as higher CRP, DBIL and GLO were associated with IVIG-resistance (Table 2). However, multivariate logistic regression analysis showed that Na+, Hb, CRP, and GLO were independently correlated to IVIG-resistance (Table 2).
3. Dose-response associations
The association between risk factors and IVIG-resistance based on RCS was analyzed, as shown in Fig. 1. Hb and Na+ exhibited a linear and negative association with IVIG-resistance, with reference points at 105 g/Land 137 mmol/L, respectively (Fig. 1A and B). GLO showed a nonlinear and positive relationship with IVIG-resistance, displaying a V-shaped curve (all values of Poverall and Pnonlinear <0.05) (Fig. 1C). The risk of IVIG-resistance increased rapidly when GLO levels exceeded approximately 22.5 g/L. CRP also showed a linear and positive association with IVIG-resistance (all values of Poverall <0.05, and all values of Pnonlinear >0.05) (Fig. 1D), with a risk threshold concentration of 75 mg/L.
4. Model constructing and evaluation
The baseline data of the training set and validation set was presented in Supplementary Table 1. There was no significant difference between the 2 groups. The XGBoost, logistic, SVM, LightGBM, and MLP models based on 4 variables (Hb, CRP, Na+, and GLO) were built within the training sets. The results showed that the area under curve (AUC) of XGBoost was the highest in both the training set and validation set (Fig. 2A and B) (Supplementary Tables 2 and 3). Furthermore, XGBoost had the highest sensitivity, positive predictive value, negative predictive value, F1 score and Kappa in the training set and the validation set. Moreover, LightGBM has the highest accuracy and specificity, followed by XGBoost in the training set and the validation set, respectively (Supplementary Tables 2 and 3). Besides, DCA, calibration curves, and PR were also used to evaluate the performance efficiency of the 5 models. The DCA showed that logistic had better net benefits, followed by XGBoost (Fig. 3C). Calibration curves displayed the best accuracy of XGBoost models (Fig. 3D). XGBoost model indicated the optimal performance, as supported by the highest AP value in both the training and validation sets (Fig. 3E and F). Collectively, these analyses demonstrated that XGBoost was considered to be the optimal model.
5. The optimal model building and evaluation
XGBoost and 10-fold cross-validation were performed on the training set. The average AUC for the training set, verification set, and testing set were 0.935 (95% CI, 0.904–0.966), 0.774 (95% CI, 0.565–0.964), and 0.821 (95% CI, 0.698–0.944), respectively, as shown in Fig. 2A–C. The AUC index of the validation set did not exceed that of the test set, or the exceedance ratio was less than 10%, indicating successful model fitting. The learning curve demonstrated a strong fit and high stability of the training set and the verification set (Fig. 2D). These results indicate that the XGBoost model is suitable for the classification modeling task with the current data.
6. XGBoost model interpretation with SHAP method
The SHAP method was utilized to explain how the selected variables predicted the likelihood of IVIG-resistance in the model. In Fig. 4A, the 4 most important features were displayed in the current model. Each feature had a significant line, and different colored dots represented the attribution of outcomes for all patients. The red and blue points indicated high and low risk values, respectively. The risk of IVIG-resistance was found to increase with decreased levels of Hb and Na+, and increased levels of GLO and CRP. Fig. 4B demonstrated that the 4 predictors were ranked based on their average absolute SHAP value. Furthermore, 2 typical cases were used to highlight the interpretability of the XGBoost model. The patient with IVIG-resistance had a high SHAP value (0.69) (Fig. 4C), while another patient without IVIG-resistance showed a low SHAP value (0.02) (Fig. 4D).
Discussion
In this retrospective study, we investigated the dependent risk factors of IVIG-resistance. We developed and validated 5 ML algorithms to predict IVIG-resistance. Na+, Hb, CRP and GLO were identified as independent risk factors for IVIG-resistance. Additionally, we observed a nonlinear relationship between GLO and IVIG-resistance. Among these models, the XGBoost model demonstrated superior performance compared to logistic, SVM, MLP, and LightGBM.
The existing prediction models for IVIG-resistance have not performed satisfactorily in other regions. In the past, logistic regression mode was commonly used to establish prediction models for IVIG-resistance, but they were not suitable for dealing with high-latitude, large-volume data. In contrast, ML algorithms could exhibit excellent performance than traditional logistic model [27,28]. In the present study, XGBoost model exhibited the best performance than other models (logistic, SVM, MLP, and LightGBM) in the internal validation set, as assessing the PR curves, calibration curves, AUC and DCA. Previous studies indicated that LightGBM model was superior to other models (Logistic regression, SVM, XGBoost) to predict IVIG-resistance [17,29], implying the excellent predictive performance of ML algorithms. Although the optimal ML model was different between our study and previous studies [17,29], it may be duo to different selection of feature variables.In this study, in order to control confounding factors, we identified risk factors based multivariate logistic regression, and then construct IVIG-resistant ML models with these risk factors. Additionally, Wang et al. [30] found that XGBoost also displayed satisfactory performance to predict 3-year all-cause mortality in patients with heart failure compared with SVM, MLP, logistic regression, k-nearest neighbors, and naive Bayesian models. The optimal model we developed could assist pediatricians in identifying high-risk patients. If the model discriminated patients with high probability of IVIG-resistance, appropriate regiments can be implemented during the initial treatment to decrease the risk of CALs. However, further investigations should be taken to confirm the performance of this model.
The utilization of ML is often impeded by the inexplicable and unobservable attributes, leading to doubts regarding the applicability of such algorithms in forecasting disease occurrences within clinical environments. As such, we estimated SHAP value of XGBoost model to interpret and visualize the prediction results, whereby establishing a score system to predict IVIG-resistance. Likewise, the application of it was found in survival prognosis of cancer [21], which served as the theoretical framework and practical foundation for related research endeavors. In present study, Hb was the most important index in IVIG-resistance, followed by GLO, CRP and Na+. All of them were independently associated with IVIG-resistance. Actually, A reduction in Hb levels may result in insufficient oxygen delivery, particularly in compromised blood vessels, thereby exacerbating vascular damage and inflammatory reactions, consequently heightening the susceptibility to cardiovascular events. Furthermore, Hb serves as a crucial diagnostic marker for anemia, the presence of which was associated with KD [31] and decreased Hb was positively associated with IVIG-resistance [32]. Our results seems to be in line with previous study indicating that low Na+ levels [33] and high CRP [7,33]. were associated with IVIG-resistance. In addition, there was a dose-dependent relationship between GLO and IVIG-resistance, whereby an increased risk of IVIG-resistance was observed when GLO exceeded 22.5 g/L. In fact, numerous studies have showed that elevated GLO levels are associated with infection [34] and mortality [35,36], thereby supporting our findings, to some extent. The comprehensive data outlined in the results and elucidation of risk factors offer doctors a deeper understanding, enabling them to make well-informed decisions rather than relying solely on algorithmic outputs. Moreover, individualized explanations can aid in the comprehension of the rationale behind the model's recommendations for critical decisions. However, further research is need to evaluate the practicality of the SHAP method.
However, there are several limitations that should not be ignored in this study. Firstly, the sample size in our study was relatively small. Secondly, selection bias was inevitable due to the retrospective nature of the study. Thirdly, The ML model we constructed was limited to one institution, which may restrict its practicality and needs validation in external datasets. Finally, although our XGBoost model can be interpreted by SHAP, it would be beneficial to explore other popular interpretable models for pediatricians to consider.
In conclusion, our results indicated that Na+, Hb, CRP, and GLO were the independent risk factors for IVIG-resistance. We constructed an IVIG-resistance prediction model using a ML algorithm, with the XGBoost model demonstrating the best performance. The interpretability method was employed to explain the XGBoost model, which has significant clinical value and can assist pediatricians in effectively managing KD patients by maximizing the identification of high-risk patients.