The MIMIC-IV is a large, free, and open database; the latest version is MIMIC-IV (version 1.0) (Johnson et al., 2021), which contains comprehensive information on approximately 250,000 patients hospitalized from 2008 to 2019, and provides strong data support for clinical studies (Yang et al., 2020; Wu et al., 2021). The database was approved by the Massachusetts Institute of Technology (Cambridge, Mass.) and the Beth Israel Deaconess Medical Center (Boston, Mass.), and consent was obtained for collection of the original data (Zhang et al., 2021). The database also anonymizes patient information, and so informed consent did not need to be obtained. The researchers needed to complete corresponding courses and obtain certificates to access and extract data from this database.

ICU patients diagnosed with VAP in the MIMIC-IV database were included. For patients admitted to the ICU more than once, only information on the first admission was obtained. Patients without a detailed record of ventilator use or who stayed in the ICU for less than 24 h were excluded. Patients were divided based on whether or not they had undergone bronchoscopy during IMV.

Data with less than 20% missing values in the database were extracted using Structured Query Language. The included demographic information were age, sex, body mass index (BMI), and race. Other information included the first care unit, disease severity [ Acute Physiology Score III (APSIII)], interventional therapy [vasopressor use, continuous renal replacement therapy (CRRT) use, and duration of IMV], major comorbidities [sepsis, myocardial infarction (MI), congestive heart failure (CHF), hypertension, cerebrovascular disease (CD), chronic pulmonary disease (CPD), liver disease (LD), renal disease (RD), diabetes, and malignancy]; microbiology; results of the first laboratory tests after ICU admission [white blood cells (WBC), neutrophils, lymphocytes, basophils, eosinophils, monocytes, red blood cells (RBC), hematocrit, hemoglobin, mean corpuscular volume (MCV), mean corpuscular hemoglobin (MCH), mean corpuscular hemoglobin concentration (MCHC), red blood cell distribution width (RDW), platelets, anion gap (AG), bicarbonate, calcium total, calcium free, magnesium phosphate, chloride, sodium, potassium, base excess (BE), calculated total CO₂, pH, PaCO₂, PaO₂, lactate, creatinine, blood urea nitrogen (BUN), glucose, international normalized ratio (INR), prothrombin time (PT), partial thromboplastin time (PTT), aspartate aminotransferase (AST), alanine aminotransferase (ALT), alkaline phosphatase (AP), lactate dehydrogenase (LDH), total bilirubin, and albumin], and vital signs within 24 h of ICU admission [mean heart rate (mHR), mean value of mean arterial pressure (mMAP), mean respiratory rate (mRR), mean temperature (mT), mean SpO2(mSpO2), urine output].

The outcomes of this study were ICU and in-hospital mortality.

In this study, Pearson’s correlation coefficient was adopted to determine the correlation between characteristic variables and to remove strongly correlated variables (Zhang et al., 2022). The following variables had coefficients higher than 0.6, so were removed since they were considered strongly correlated: RBC, hematocrit, MCV, chloride, BE, calculated total CO₂, pH, BUN, PT, AST, and LD (Supplementary Figure S1). For other variables with less than 20% of missing values, the “mice” package of R software was used to perform multiple imputation.

The data for continuous variables were represented as mean and standard-deviation or median and interquartile (IQR) values, while those for categorical data were represented as frequencies. The Mann-Whitney U test was used for continuous variables, and the χ² test or Fisher’s exact test was used for categorical variables.

The Kaplan-Meier (KM) method draws a cumulative incidence curve, indicating the occurrence of ICU and in-hospital deaths in different groups of patients, and the differences in risk between the groups were compared using log-rank tests. Furthermore, after adjusting different covariates, two Cox proportional-hazards models were constructed to analyze the influence of the relationship between bronchoscopy and outcomes. There were no adjustments for covariates in model I. In model II, all variables, including age, sex, race, BMI, first care unit, APSIII, vasopressors, CRRT, duration of IMV, sepsis, MI, CHF, hypertension, CD, CPD, LD, RD, diabetes, malignancy, WBC, neutrophils, lymphocytes, basophils, eosinophils, monocytes, hemoglobin, MCH, MCHC, RDW, platelet, AG, bicarbonate, calcium total, calcium free, magnesium phosphate, sodium, potassium, PaCO2, PaO2, lactate, creatinine, glucose, INR, PTT, ALT, AP, bilirubin total, albumin, mHR, mMAP, mRR, mT, mSpO2, urine output were adjusted for confounding factors. In the multivariate COX regression, we also evaluated the multicollinearity between the variables using variance inflation factors (VIF). Supplementary Table S1 revealed that the VIF of each variable was less than 4, indicating that there was no multicollinearity between them.

To ensure that the results were stable and reliable, we further adjusted for covariates using propensity score matching (PSM) and propensity score-based inverse probability of treatment (IPTW) after analyzing the original population. Multivariate logistic regression model was used to estimate patient propensity scores by using one-to-one nearest neighbor matching with a caliper width of 0.05. The IPTW model was created using estimated propensity scores as weights. Differences in baseline levels between the two groups were evaluated using p values. Then, COX regression was performed on the matched population and the weighted population, respectively.

Causal mediation analysis (CMA) can distinguish the total effect of treatment into direct effects and indirect effects, If the independent variable X has a certain influence on the dependent variable Y through a certain variable M, then M is called the mediating variable of X and Y (Imai et al., 2010). In the present study, we hypothesized that changes in a particular indicator were mediating variables, that is, we assumed that bronchoscopy might lead to changes in that indicator and that such changes were associated with the prognosis of VAP patients. The average causal mediating effect (ACME), average direct effect (ADE) and total effect obtained through CMA can help us verify the above conjecture. Finally, WBC, platelets and RDW, three indicators that were further investigated for the presence of mediating effects because of their low absence rate after the cessation of IMV. Their change level was expressed as the difference between the corresponding first examination results and recorded as delta-WBC, delta-RDW, delta-platelet.

We also analyzed the impact of bronchoscopy on the prognoses of patients from different subgroups. The subgroups included age (dichotomized at 65 years), sex, type of first care unit, APSIII (dichotomized at the median of 69), and all comorbidities. The interactions between subgroups were further analyzed.

A two-tailed probability value of p < 0.05 was considered statistically significant. All statistical analyses in this study were performed using R software (version 4.1.0).

This study enrolled 1,560 patients with VAP: 1,355 in the no-bronchoscopy group and 205 in the bronchoscopy group (Figure 1). The median (IQR) ages of the patients in these two groups were 64.00 years (53.00–75.00 years) and 62.00 years (50.00–74.00 years), respectively; there were more male patients in the two groups (62.7 and 65.9%, respectively); the median (IQR) APSIII were 69.00 (51.00–90.00) and 70.00 (56.00–91.00), respectively; the median (IQR) durations of IMV were 136.00 h (60.50–253.50 h) and 145.00 h (59.00–286.00 h), respectively; and sepsis was the most common complication, present in 96.5 and 97.6% of patients in the two groups, respectively. The remaining baseline characteristics of the patients are listed in detail in Table 1.

The KM survival curve in Figure 2 indicated a significant difference in survival probability between the two groups. The survival probability in the ICU and in hospital was significantly higher in the bronchoscopy group than in the no-bronchoscopy group (Figure 2). The results of the log-rank test indicated that the mortality risks in the ICU and in hospital differed between the two groups.

Table 2 lists the results for the Cox proportional-hazards model. In two models with no adjustment for confounders, or adjustment for all confounders, the hazard ratio (HR) was significantly lower than one for the bronchoscopy group compared with the no-bronchoscopy group. In other words, compared with patients in the no-bronchoscopy group, those in the bronchoscopy group had a lower risk of ICU and in-hospital mortality. In model II, after adjusting for all covariates as confounding factors, the HR (95% CI) values for ICU and in-hospital mortality in the bronchoscopy group were 0.33 (0.20–0.55) and 0.40 (0.26–0.60), respectively, indicating that the risks of ICU and in-hospital mortality were 0.67 and 0.60 lower than in the no-bronchoscopy group. (Table 2).

After PSM and IPTW, there was no significant difference in baseline levels between the two groups (Supplementary Table S2). Kaplan-Meier survival curves of matched and weighted populations are consistent with the original population (Figure 2). Univariate and multivariate COX regressions were then performed on the matched and weighted populations, respectively, yielding results consistent with the original population (Table 2). After multivariate COX regression, HR (95% CI) values for ICU mortality in the bronchoscopy group were 0.33 (0.17,0.64) and 0.26 (0.14,0.50), and in-hospital mortality the HR (95% CI) values of 0.33 (0.20,0.58) and 0.35 (0.21,0.57), respectively.

In the CMA analysis, after analyzing the three indicators, only the changes of RDW were significant. Figures 3A, B showed that in terms of in-hospital mortality of VAP patients, delta-RDW mediated 8% (95% CI:1–26%; p = 0.02) of the beneficial effect of bronchoscopy (ACME:p = 0.02). At the same time, Figures 3C, D implied that in terms of ICU mortality, delta-RDW mediated 7% (95% CI: 1–21%; p = 0.02) of the beneficial effects of bronchoscopy (ACME:p = 0.02).

There was no significant difference between the subgroups of in-hospital mortality, meaning that no interaction was present. However, bronchoscopy was associated with a reduced risk of ICU mortality in MI and CPD populations, with HRs (95% CI) of 0.21 (0.05–0.86) and 0.13 (0.05–0.34), respectively. Meanwhile, patients without LD were more likely to benefit from bronchoscopy than those with LD (HR = 0.29, 95% CI = 0.16–0.54) (Table 3).