Using flexible methods to determine risk factors for ventilator-associated pneumonia in the Netherlands

Seven hospitals participated in the Dutch national surveillance for ventilator-associated pneumonia (VAP) and its risk factors. We analysed time-independent and time-dependent risk factors for VAP using the standard Cox regression and the flexible Weighted Cumulative Effects method (WCE) that evaluates both current and past exposures. The prospective surveillance of intensive care patients aged ≥16 years and ventilated ≥48 hours resulted in the inclusion of 940 primary ventilation periods, comprising 7872 ventilation days. The average VAP incidence density was 10.3/1000 ventilation days. Independent risk factors were age (16–40 years at increased risk: HR 2.42 95% confidence interval 1.07–5.50), COPD (HR 0.19 [0.04–0.78]), current sedation score (higher scores at increased risk), current selective oropharyngeal decontamination (HR 0.19 [0.04–0.91]), jet nebulizer (WCE, decreased risk), intravenous antibiotics for selective decontamination of the digestive tract (ivSDD, WCE, decreased risk), and intravenous antibiotics not for SDD (WCE, decreased risk). The protective effect of ivSDD was afforded for 24 days with a delay of 3 days. For some time-dependent variables, the WCE model was preferable over standard Cox proportional hazard regression. The WCE method can furthermore increase insight into the active time frame and possible delay herein of a time-dependent risk factor.

Some of these patient and treatment characteristics are time-dependent. Longer exposure or treatment will modify the (cumulative) risk, but it is less clear and often not evaluated how the association of these time-dependent risks develop during and following exposure [9,10,12,16,17,[19][20][21]. In this manuscript we present the VAP surveillance results and evaluate the risk factors. For the time-dependent risk factors we use both standard Cox regression and the flexible Weighted Cumulative Effects (WCE) approach that evaluates both current and past exposures [22]. The WCE approach allows estimation of the timeframe during which a risk factor is (still) relevant.

Study setting and set up
PREZIES (Prevention of HAI through Surveillance) is the Dutch national surveillance network for healthcare associated infections, hosted by the National Institute for Public Health and the Environment (RIVM), in which hospitals participate voluntarily. Since 1996, PREZIES offers Dutch hospitals the possibility to participate in the surveillance of hospital-acquired infections with attendant benchmarks. The VAP surveillance module was offered from January 2004 -December 2011.
The VAP surveillance protocol was developed by a working group of relevant professionals (intensivist, infectious disease specialist, pulmonologist, anaesthesiologist, epidemiologist and infection control professional) and was based on international literature, the preceding ICUsurveillance conducted by PREZIES [23] and a pilot study. The VAP definition used (Fig 1) was a simplified version of that used by the former Centers for Disease Control and Prevention/National Nosocomial Infections Surveillance system, currently the National Healthcare Safety Network definitions. In Dutch clinical practice, a VAP diagnosis is most often based on positive cultures of tracheal aspirates, in combination with clinical symptoms. With a negative culture, the diagnosis of a clinical pneumonia was still possible.
All patients ventilated invasively for two days (48 hours) or more, aged 16 years or older, and present in the ICU were included in the study. Successive ventilation periods of at least two days were recorded until a VAP occurred or until the end of follow-up, 28 days after the start of each ventilation period. Each ICU admission was assigned a unique identifier that could not be used to link admissions at the patient level. Therefore, patients admitted more than once to the ICU were included as separate patients (termed "admissions"). For this manuscript, we considered only data from the first ventilation period of each admission. Data were recorded prospectively, with infection control professionals checking the patients and patient records in the ICU on average twice a week. Suspected VAPs were usually discussed with a dedicated radiologist or intensivist.
• Per admission were recorded: sex, age, admission and discharge dates of hospital and ICU, Apache II score, specialism, type of ICU and reason for end of follow-up.
• Per infection were recorded: infection date, infection criteria and a maximum of three micro-organisms. When a clinical or possible pneumonia was later followed by respectively a possible and/or ('confirmed') pneumonia (for the same infectious episode), the latter were recorded as final diagnosis. In this paper, we do not distinguish between these diagnostic categories. According to Dutch legislation, written consent from each individual patient was not required because the data from the PREZIES network is anonymized and was collected as a legal task of the National Institute for Public Health and the Environment.

Statistical analysis
Data are expressed as median (interquartile range (IQR)) and absolute and relative frequencies, as appropriate. A non-parametric test (Kruskal-Wallis) was used to compare median ventilation durations.
We first calculated cause-specific hazards using univariate Cox regression models, adjusted for hospital as a fixed effect. Time-fixed, non-linear, or continuous covariates as well as duration of participation were modelled categorically in separate models using dummy variables. Separately for each time-dependent covariate, we fitted four regression models to determine if, and in which form, the covariate should be included in a multivariate analysis. The four models included three Cox models (current effects and delayed effects of one or two days) and one flexible Cox model using the WCE approach [22]. The WCE approach estimates not only the current effect of a covariate, but also the cumulative delayed effects of past exposures, and provides the timeframe for which a covariate is significantly associated with the studied event. See Box 1 for more details [25]. We assessed which of these four models fitted best using the Akaike information criteria (AIC) [26]. A difference of less than four suggested that the models fit the data equally well, a difference between four and ten suggested a slight difference, and a difference greater than ten suggested a major difference in model fit [27]. We chose the WCE model as the best-fitting model when it yielded a slight or major improvement in model fit to all of the three non-cumulative models.
All risk factors with p-value<0.2 in the univariate analysis were included in the initial multivariate Cox PH model, using the best-fitting univariate models. The final multivariate model was selected manually by backward selection using the likelihood ratio test. At each iteration, we removed from the model the variable that was associated with the highest p-value>0.05, except when the 95% CI for that variable did not include the null.
The WCE model requires complete data for the entire follow-up of each admission. Since some time-dependent data were missing for only 32 ventilation days (0.4%) in 24 admissions (2.6%), we used the last observation carried forward approach to fill in these values for missing days in our regression analyses. The sensitivity of our results to this approach for completing missing data was assessed by removing admissions with missing time-varying data from the dataset and rerunning the multivariate model. Admissions with data missing on time-fixed patient or ventilation period characteristics were excluded from the specific analysis.
In the Netherlands, three ICU levels are discerned, ranging from I (relatively low complex care) to III (high complex care). The participating hospitals had ICU level I (1), II (3) and III (3) (see Table B in S1 File for hospital-specific results). Median ventilation duration differed significantly between hospitals (p<0.0001), ranging from four to eight days. Four of the five hospitals that participated for two years or more demonstrated a reduction in VAP (Fig 2), of which two significantly so.

Ventilation periods
The surveillance included 940 first ventilation periods of 940 ICU-admissions, all in mixed medical/surgical ICUs, comprising 7872 ventilation days, including all calendar days. The median ventilation duration for all admissions was 6 days (interquartile range (IQR) 4-10 days) (Table B in S1 File); for patients without VAP, 6 days (4-11); and for those that developed VAP, 5 days (3-7) until VAP. During follow-up, the number of patients that were still on the ventilator each day declined exponentially (Figure A in S1 File). Of the 940 admissions, 81 developed a VAP during follow-up (8.6%) and 23 were still on the ventilator on day 28. The average VAP incidence density was 10.3/1000 ventilation days (range between hospitals 0.0 to 20.1).

Box 1. The WCE method in more detail
The WCE model uses the Cox PH framework and time-varying covariates to generate, for each covariate, a function that describes the delayed and immediate effect of (past) exposures/levels on the outcome [22].The WCE model requires as input (1) the timewindow in which past exposures are considered to have an impact on the risk of the outcome, (2) a pre-specified number of internal knots, which determine the flexibility of the cubic B-spline, and (3) whether the impact of past exposures reaches zero at the earliest point of the time-window (i.e., constraining the effect of the covariate to the null at that point). When insufficient prior information is available to make an informed choice on these inputs, as in our case, the data may be used to determine which inputs provide the best model fit. The approach used to select the optimal WCE model for each factor involved fitting multiple WCE models using all possible time-windows (up to 2, 3, 4, . . ., 28 days back); 1, 2 or 3 internal knots; and with the effect of the exposure at the most historical time point included in the time-window either unconstrained or constrained to the null.
From these 162 (27 x 3 x 2) models, the best-fitting WCE model for each factor was selected as the WCE model with the lowest Akaike information criteria (AIC). Since we selected the best of multiple models, p-values for WCE univariate models are likely to be artificially low. For each factor, we therefore simulated 1000 datasets in which there was no association, keeping the exposure patterns and outcome times consistent with the original dataset. For each of the 1000 simulated datasets, we ran the same set of WCE models (with alternative time-windows, numbers of internal knots, and weight-function constraints), and selected the best-fitting model using the AIC, as above. The distribution of the 1000 p-values was plotted and the proportion of the 1000 p-values that were smaller than the p-value of the optimal WCE model was recorded as the p-value corrected for multiple testing [25]. and men developed VAP more often (10.1%; 95% confidence interval (CI) (8.2-12.0%) versus 6.5% (4.9-8.1%)). The median age was 69 (IQR 59-77) years for all admissions and the median Apache II score 21 (IQR [16][17][18][19][20][21][22][23][24][25][26][27]. See Tables 1 and 2 and Table B in S1 File. The median age and Apache II scores were comparable for admissions with and without a VAP. The routines and regimens offered to patients differed between hospitals. Apart from one hospital, hospitals preferred either MDI or nebulizers, four of them exclusively. IvSDD was used in three hospitals, intestinal prophylaxis in three, and SOD in five (see Figure B in S1 File). In hospital G, SOD was given for one year, followed by a complete SDD regimen (oropharyngeal and intestinal prophylaxis and four days ivSDD), as part of a study. In hospital E, SOD was initially given occasionally and subsequently extended to all admissions.

Univariate results
Of the time-independent covariates, age, length of hospital stay before start of the ventilation, COPD, intubation department, and duration of participation with the surveillance (in years) were significantly associated with VAP (Table 1). All time-dependent variables were significantly associated with the risk of acquiring VAP except for type of inhalation therapy ( Table 2). Fig 3 shows the WCE results for the univariate analyses. The overall hazard ratio of specific exposure patterns compared with uniformly non-exposure can be calculated by multiplying the hazard ratios for all days where, for example, systemic antibiotics ( Fig 3C) were administered. Suppose a patient was ventilated for 6 days and treated with ivAB on the first four ventilation days (day -5 to -2), but not on the last two ventilation days (day -1 and 0). The hazard ratio for this patient on the present day (day 0) is a multiplication of the hazards on days -5, -4, -3, and -2. Note that the hazard ratio before day -4 in Fig 3C is assumed to be 1. The WCE model for ivAB suggests a harmful effect of antibiotics taken on the present day (day 0), which is probably a result of reverse causation since system antibiotics were likely taken on the day of a VAP to treat the pneumonia. Further, the narrower confidence intervals shown at the left (e.g. Fig 3A) may be seen as counterintuitive since the number of patients declines with   increasing follow-up duration ( Figure A in S1 File), but are an artefact of the analysis since the best-fitting model for ivAB was one where the hazard ratio was constrained to the null at this point.

Multivariate results
After backward selection to step-wise exclude covariates from the model, the final multivariate model included age, COPD, current sedation score, current SOD, inhalation therapy (WCE), ivSDD (WCE), and ivAB (WCE) ( Table 3). The model suggests that, compared to patients between 40 and 60 years of age, patients � 40 had an over two-fold higher risk of developing VAP. Patients were also at an increased risk when their sedation scores were high. Patients who were on SOD, ivSDD or ivAB and patients with COPD had a lower VAP risk. The model further demonstrates a delayed protective effect of ivSDD by three days. This protective effect is afforded for 24 days (up to 27 days back) (Fig 4A). IvAB was protective for two days, with the protection being afforded after a delay of one day. The effect from inhalation therapy was minimal, with jet nebulizers showing a delayed protective effect for days 6 to 8 back. Our results were not sensitive to using the last observation carry forward approach for completing missing data. Similar results were obtained when we removed admissions with missing time-varying data from the dataset and reran the multivariate model.

Discussion
In this manuscript, we presented data generated by the VAP surveillance module of the Dutch PREZIES network. The VAP rates reported by the hospitals participating in the surveillance (0-20.1/1000 ventilation days) are in the expected range [1][2][3]. The observed variation between hospitals could be partly due to the inter-observer reliability of a VAP diagnosis, which is known to be low [31,32]. However, this applies less to intrahospital comparisons where only one intensivist or radiologist was usually dedicated to surveillance. Four of the five hospitals that participated �2 years demonstrated a reduction in VAP, two of which significantly so. Apart from a possible effect of the surveillance itself [33], various changes during the follow-up could have caused this reduction. Several hospitals (C-G) implemented interventions (SOD �� The number of ventilation days after reducing the numbers of missings for time-dependent covariates with the 'last observation carried forward (LOCF) approach ��� Selected as the model with the lowest AIC among the current exposure-risk, 1-day delay exposure-risk, and 2-day delay exposure-risk models in Table C in S1 File.
���� Preference given to a non-cumulative model. WCE model selected if the WCE model had an AIC that was 4 units lower than that of the best-fitting non-cumulative model. ����� The p-value was estimated using 1000 bootstrapped data sets to account for multiple testing when selecting the best-fitting WCE model. and/or ivSDD, Evac cuffs, closed suctioning system, VAP bundle), sometimes temporarily. Although the increase in VAP incidence in hospital F was duly investigated at that time, no cause was found. Hospital A did not introduce any interventions. Hospital B, with zero VAPs, used a complete SDD regimen during the surveillance. Therefore, while hospitals appear to be able to reduce VAP, either by introducing an intervention or by surveillance alone, success does not seem to be guaranteed. The observed effects of patient and treatment characteristics vary among studies [9,[11][12][13][14][15][16][17][18], resulting from differences in the other measured covariates and case mix, and, frequently, low statistical power. In our results, Apache II score, specialty, intubation site, length of hospital/ ICU stay before ventilation, postsurgical ventilation, inhalation trauma, stress ulcer The curve shows the estimated risk attributed to exposures on each day prior to the last day of follow-up (i.e., the event date or the censoring date) and the grey ribbon shows the 95% confidence interval. A value of one indicates no effect of the exposure at that time. At times where the grey ribbon includes one, the effect is considered to be statistically insignificant. https://doi.org/10.1371/journal.pone.0218372.g003 Risk factors for VAP with flexible methods prophylaxis, corticosteroids, neutropenia, intestinal prophylaxis, and feeding method were not significantly associated with VAP. Although the overall study population in our study was relatively large, low numbers of patients in certain categories could explain failures to detect associations. In our data, a higher Ramsay score was associated with an increased VAP risk, which corresponds with the increased risk associated with coma or increased Glasgow coma scale [9-12, 16, 17]. Systemic antibiotics (not for ivSDD) appeared to lower the VAP risk, as identified by others [9,10]. In most analyses where antibiotic use was not analysed as a time-dependent variable, 'prior antibiotic use' (or variations) was not found to be associated with VAP [12,16,17]. Our results demonstrated that both ivSDD and SOD were associated with a VAP reduction. De Smet et al, performing a multicenter randomized clinical trial, concluded that both SOD and SDD (full regimen) led to a reduction in respiratory tract colonisation with highly resistant microorganisms [21]. As intestinal prophylaxis partly coincided with ivSDD and SOD, we evaluated a model without ivSDD/SOD. Intestinal prophylaxis was not significantly associated with VAP in this model either.
More surprisingly, because unlike other studies where the risk for older patients was similar or increased [3,11,17,18] younger patients (16-40 years) appeared to be at increased risk. This may have resulted from an overrepresentation of young patients in neurosurgery and traumatology, both specialties with high VAP rates. This association remained borderline significant when adjusted for specialty and Apache II score or when these two specialties were excluded from the analysis. When excluding one center, where 11 of the 32 patients aged 16-40 developed pneumonia, the association was still present (HR 3.2 (0.49-21.2), although not significant anymore. While COPD is generally found to be either unassociated with VAP [3,16], or to increase the VAP risk [11,17], here COPD appeared to be associated with a lower risk. Our model did not detect a significant interaction with systemic antibiotics. In previous studies, COPD patients had longer ventilation durations relative to patients without COPD [3,34]. In our study, the first ventilation periods for patients with and without COPD were similar Risk factors for VAP with flexible methods in duration (average of 8.2 and 8.4 days, respectively), as was the total ventilated duration (9.3 and 8.7 days) and the proportion of ventilation periods lasting 28 days(1.4 and 2.8%, respectively), suggesting that the similarities in ventilated duration are not subject to censoring bias. Furthermore, mortality data, available from five hospitals, showed that mortality was comparable (40.4% and 38.5%, respectively). The lower risk in our study may reflect differences in criteria for COPD diagnosis [35], but also improved care. For example, Funk analysed the outcomes of COPD patients admitted in ICU, for the period 1998-2008, and found that riskadjusted mortality had improved [36].
The use of jet nebulizers demonstrated a slightly reduced risk compared to no inhalation therapy, whereas MDI did not. Appropriate delivery of medication is reported to be more challenging with MDI than with jet nebulizers. In one study no difference in VAP risk was found between both methods, but the population size was small [37].
For the daily measured ivAB, ivSDD, and inhalation therapy, the WCE approach explained the associated VAP risk better than the standard Cox regression analysis. Additionally, this method provides insight into the cumulative effect over time and the relevant retrospective timeframe in which exposures are associated with a VAP. IvSDD affected this hazard for 24 days. It must be noted that in one hospital many patients received ivSDD throughout the entire ventilation period instead of the currently recommended first four days only. Although differences in the AIC were usually small (Table C in S1 File) the WCE approach identified ivSDD as an

B) ivAB, (C) Jet nebulizer (compared to no inhalation therapy) and metered dose inhalers (MDI) (compared to no inhalation therapy).
The curve shows the estimated risk attributed to exposures on each day prior to the last day of follow-up (i.e., the event date or the censoring date) and the grey ribbon shows the 95% confidence interval. A value of one indicates no effect of the exposure at that time. At times where the grey ribbon includes one, the effect is considered to be statistically insignificant. https://doi.org/10.1371/journal.pone.0218372.g004 Risk factors for VAP with flexible methods independent factor affecting the risk to develop VAP whereas standard Cox regression of current values or of one or two days earlier did not. The WCE approach is worth considering in future analysis of time-dependent risk factors of VAP, and for other device-associated infections.
Our study has some limitations. First, while the relatively large study population and the availability of daily data on time-varying risk factors represent strengths, our study made use of surveillance data and therefore missed more detailed clinical data, e.g. type of ivAB or the presence of sepsis, that could also be associated with VAP and affect estimated associations for other factors. Second, hospital-level treatment preferences resulted in low variability within individual hospitals, reducing the power to detect associations on these variables. Third, our results may not be generalizable to other settings. VAP in the Netherlands is usually diagnosed and treated based upon clinical features and/or tracheal aspirate cultures. Therefore, some risk factors may differ when VAP is diagnosed in a more invasive way. Top clinical hospitals were overrepresented in the sample of participating hospitals, whereas academic hospitals did not participate. Along with the variation in VAP rates observed this implies the average VAP incidence density may not be representative of all hospitals. Fourth, we could not distinguish between new and recurrent ICU admissions, which may have led to some overrepresentation of complex patients. Fifth, including only the first ventilation episode may have led to some bias in the assessment of risk factors. Lastly, in-house infection control professionals collected the data. Apart from the low inter-observer reliability of diagnosing VAP [31,32], this may have led to differences in the application of the surveillance protocol. We aimed to minimize this possible bias by arranging meetings for the involved professionals to discuss data collection and infection criteria.
Surveillance results are limited for clinical use or pathophysiological insight but the data of Dutch ICUs participating in the VAP surveillance system revealed risk factors on both patient (age, sedation score) and treatment level (SDD, oropharyngeal prophylaxis, other antibiotics, nebulizer type) that can be useful for case mix adjustment and evaluation of VAP prevention strategies. The introduction of SDD or oropharyngeal prophylaxis was associated with low or zero VAP incidences. Surprisingly, COPD was associated with a reduced VAP risk, which merits further evaluation. For some time-dependent covariates, the WCE approach was preferable over standard Cox proportional hazard regression and additionally provided insight into the relevant retrospective timeframe of past exposures.