The Medical Information Mart for Intensive Care (MIMIC)—IV database, specifically versions 2.2 and 3.0, was used for this retrospective analysis (9, 10). MIMIC-IV is a publicly accessible, large-scale intensive care database organized and maintained by the Laboratory for Computational Physiology at the Massachusetts Institute of Technology (MIT). This study utilized version 2.2, which includes medical data from approximately 300,000 patients treated at the Beth Israel Deaconess Medical Center (BIDMC) from 2008 to 2019, for model training and internal validation. Version 3.0, which includes data from 2020 to 2022 was utilized for model external testing. The use of MIMIC-IV data was ethically approved by the Institutional Review Boards of BIDMC and MIT. Since all personal data in the database are anonymized, informed consent was waived. The author (Heshan Cao) was granted access to the database (certification number: 63137030). Reporting of this study followed the Transparent Reporting of a Multivariable Prediction Model for Individual Prognosis or Diagnosis (TRIPOD+AI) guidelines (Supplementary Table 1) (11).

This study included patients aged 18 and above who were admitted to the ICU primarily for IS and had MV for more than 48 h. The primary outcome measure was the incidence of VAP. Patients with IS and VAP were identified in the MIMIC-IV database using the International Classification of Diseases, Ninth Revision (ICD-9), or Tenth Revision (ICD-10) codes. By restricting the diagnostic sequence to IS > VAP, we ensured that the diagnosis of IS was prioritized over that of VAP. Supplementary Table 2 lists the relevant ICD codes. Exclusion criteria included VAP diagnoses preceded IS and deaths within 7 days of ICU admission. For patients who had multiple ICU admissions, only the initial admission was considered. Figure 1 depicts a flowchart of inclusion and exclusion criteria.

The MIMIC-IV database was queried using Structured Query Language (SQL) to extract features such as demographics, comorbidities, vital signs, laboratory test indicators, and ventilator settings. Records within the first 24 h of MV were used to extract data on vital signs, laboratory indicators, and ventilator settings; variables with multiple records were averaged. Records were also acquired for antibiotic usage, suctioning procedures, and invasive catheter placements during the first 24 h of MV. We encoded comorbidities and VAP incidence as binary values.

A total of 59 features were obtained. Features with more than 20% missing data were first excluded to reduce missing data bias. Features with a missing data rate below 20% were addressed using multiple imputation methods, as indicated in Supplementary Figure 1. Following that, a correlation analysis was conducted on all features, and those with a correlation coefficient greater than 0.7 were excluded to prevent multicollinearity from impacting model performance (Supplementary Figure 2). Finally, the Boruta algorithm was applied to select the most relevant features. Boruta is an all-relevant feature selection method that uses a random forest (RF) classifier to compare the importance of original features with that of randomly permuted “shadow” features. By iteratively eliminating features that do not outperform their shadow counterparts, the algorithm robustly identifies truly informative features (12). In this study, Boruta was executed with a confidence level of 0.01, iterated 500 times to exclude rejected features, as presented in Figure 2.

The MIMIC-IV database data, which ranges from 2008 to 2019, was randomly divided into two sets: 80% for training and 20% for validation, using a stratified sampling strategy. To predict the risk of VAP in IS patients, 10 widely recognized ML models based on different principles were constructed: adaptive boosting (AdaBoost), category boosting (CatBoost), extra trees (ET), light gradient boosting machine (LightGBM), logistic regression (LR), multilayer perceptron (MLP), naive Bayes (NB), RF, support vector machine (SVM), and extreme gradient boosting (XGBoost). This comprehensive approach allowed us to identify the model that best balances performance for early VAP prediction in IS patients. To optimize the prediction models and avoid overfitting, the final hyperparameters for each model were determined using a combination of five-fold cross-validation and Bayesian search. For external validation, the trained models were tested on the MIMIC-IV database version 3.0, which covers the years 2020–2022.

The models’ performance was evaluated using metrics such as the area under the receiver operating characteristic (ROC) curve (AUC), accuracy, sensitivity, specificity, positive predictive value (PPV), and negative predictive value (NPV). The optimal cutoff value was determined by maximizing the Youden index (sensitivity + specificity − 1). Calibration and decision curves were used to evaluate the model’s calibration and clinical decision-making capability.

The SHAP method quantifies the contribution of each input feature to the final prediction by leveraging concepts from cooperative game theory, addressing the “black-box” nature of ML models (8). This approach incorporates global and local explanations. Global explanations reveal the impact of features on the overall model, whereas local explanations examine the contribution of features in specific samples. The decision-making process of the final model is visually depicted, with both global and local explanations provided via the SHAP method.

To facilitate its application in clinical contexts, the final prediction model was integrated and released into a web application built with the Streamlit Python library. When values are entered into the required features settings, the application generates a risk score for VAP in individual IS patients, as well as a force plot showing the effect of each feature on the risk assessment.

Data preprocessing, model construction, performance evaluation, and result visualization were conducted using Python (version 3.9.18) and R (version 4.3.2). The variance inflation factor (VIF) was employed to assess potential multicollinearity among the selected features, as detailed in Supplementary Table 3. Continuous variables with a normal distribution were presented as means and standard deviations, while non-normally distributed data were reported as medians (m) and inter-quartile range (IQR). Categorical variables were presented as numbers (n) and percentages (%). Differences in continuous variables were compared using Student’s t-tests or Wilcoxon rank-sum tests, while categorical variables were analyzed using Chi-square tests or Fisher’s exact tests. Statistically significant differences were defined as two-tailed p values <0.05.

Using ICD-9/10 codes and the defined inclusion–exclusion criteria, we extracted data on 401 and 118 IS patients from versions 2.2 and 3.0 of the MIMIC-IV database, respectively. In the derivation cohort, 149 IS patients developed VAP, while 64 patients in the test cohort experienced VAP. Table 1 outlines the baseline characteristics of both cohorts. The baseline characteristics of the VAP and non-VAP groups in the derivation cohort are detailed in Supplementary Table 4.

Six of the initial 59 potential predictive features were removed due to a missing rate of more than 20%, while 7 were excluded due to a correlation coefficient greater than 0.7, as shown in Supplementary Table 5 and Supplementary Figure 2. Additionally, 39 features were removed during the Boruta algorithm phase. After considering clinical practicality, 7 features were eventually included to construct the ML models, including Systolic Blood Pressure (SBP), Diastolic Blood Pressure (DBP), International Normalized Ratio (INR), Length of Stay before Mechanical Ventilation (LOS Before MV), Dysphagia, Antibiotic counts, and Suctioning counts, as illustrated in Figure 2.

The performance of the 10 ML models is presented in Table 2 and Supplementary Table 6. In the derivation cohort, the NB model performed best with an AUC of 0.790, followed by the RF and LightGBM models, both of which had an AUC of 0.776. In the test cohort, the RF model exhibited the best generalization capability with an AUC of 0.644. Figure 3A compares the AUC values of the 10 ML models in the internal validation and external test sets, while Figures 3B–D present the ROC curve, calibration curve, and decision curve of the final model.

The SHAP method was employed to interpret the final model. Figure 4A presents a bar plot of features ranked by their mean absolute SHAP values. Figure 4B shows a beeswarm plot, illustrating the relationship between each feature’s value and the predicted risk of VAP. These plots indicate that antibiotic counts, LOS before MV, and INR are the top three contributors. Figure 4C shows SHAP dependence plots that further reveal the distribution of each feature and its global relationship with VAP risk.

In addition, the SHAP method was used to conduct local interpretation for the final model. Figure 5 shows detailed local interpretations using SHAP waterfall plots and force plots. As shown in Figure 5A, for patient who did not eventually develop VAP, the SHAP analysis indicates that higher INR, the administration of antibiotics twice, a shorter LOS before MV, normal SBP and DBP, the absence of dysphagia, and fewer suctioning operations negatively contribute to the model’s prediction of VAP, resulting in a low-risk classification. In contrast, as depicted in Figure 5B, for patient who eventually developed VAP, the SHAP analysis reveals that multiple suctioning operations, higher SBP, and the presence of dysphagia positively support the prediction of VAP, while the administration of antibiotics four times and a shorter LOS before MV have a negative impact. The aggregated contributions of these features incline the model toward predicting a high risk of VAP for the patient.

Based on the final RF model, an interactive web-based tool was developed to facilitate clinical application.¹ Clinicians can input patient-specific clinical parameters to obtain an individualized prediction of VAP risk, along with a SHAP force plot clearly depicting each feature’s contribution. As illustrated in Figure 6, red features on the left side like LOS Before MV, suctioning counts, INR, and dysphagia push the prediction toward “VAP,” while the blue features on the right side like antibiotic counts, DBP, and SBP drive the prediction toward “non-VAP.” For the selected scenario shown in Figure 6, the model predicts a 68.51% probability of VAP occurrence, indicating a high risk for VAP.