Short-term effects of air pollution on hospital admissions for cardiovascular diseases and diabetes mellitus in Sofia, Bulgaria (2009–2018)

Data on hospital admissions in Sofia from 1 January 2009 to 31 December 2018 were provided by the National Health Insurance Fund of Bulgaria and we classified those of interest for our study by the WHO International Classification of Diseases (ICD-10) codes, as follows: ischaemic heart diseases (IHD; I20–25), cerebral infarction (CI; I63), and type 2 diabetes mellitus (T2DM; E11). As the main outcomes we used the total daily number of hospital admissions per diagnosis. Then, to identify potential gender and age-related differences, we stratified the number of cases as follows: men <65 years, men ≥65 years, women <65 years, and women ≥65 years, in line with earlier studies (23). We used only the publicly available patient data, and the study was therefore not subject to an approval by a bioethics committee.

Hourly concentrations of particulate matter with a diameter of less than 2.5 (PM_2.5) and 10 µm (PM₁₀), nitrogen dioxide (NO₂), sulphur dioxide (SO₂), ozone (O₃), and carbon monoxide (CO) were obtained from the Executive Environment Agency of Bulgaria. From these, we calculated daily average concentrations for the period of interest to match the health data (after purging the data from outliers). To control for potentially confounding weather conditions, we also calculated mean daily temperature and relative humidity measured at the same monitoring site. Pollutant concentrations and meteorological variables had been measured by an official background urban monitoring station (located at 23.296338, 42.680806; “Hippodrome” Park), which was the only valid source of data on these pollutants (including on PM_2.5) in Sofia for the entire 10-year study period. Automatic air quality monitoring stations such as this are accredited under the BDS EN ISO/IEC 17025 norm “General requirements for the competence of laboratories for testing and calibration” and maintain a management system ensuring the validity and reliability of results. These stations are inspected and calibrated every three months with standard gases according to the procedures set up by the Executive Environment Agency.

The data were pre-processed to create several variants of exposure variables. Pollutant concentrations, except for CO, were re-scaled so that one-unit increase corresponded to 10 µg/m³, while CO was left in its original scale of 1 mg/m³. To account for potential non-linear effects of meteorological variables, we modelled temperature and relative humidity as deciles. Delayed effects of air pollution were explored by constructing variants of each pollutant with one- to seven-day lagged values, as well as by moving average concentrations over three and seven consecutive days. To control for time trend in the outcome data, we constructed flexible splines with 69 knots and indicator variables for month of the year and year. Since hospital admissions are more likely to be registered on certain days of the week owing to administrative reasons, we also considered what day of the week each date in the time series corresponded to.

For general patterns the data were inspected visually, with descriptive statistics, and Spearman’s correlations. Negative binomial models were used because of overdispersion in all outcomes. Presence of zero-inflation and the order of autocorrelations were checked to determine the optimal model for each outcome. For IHD and CI, we fitted negative binomial regressions with Newey-West standard errors to correct for autocorrelations of up to four days. For T2DM, we preferred a zero-inflated negative binomial model. All models were controlled for the time trend, day of the week, temperature, and relative humidity. If inspection of deviance residuals from the model indicated residual autocorrelation, the model was rerun with an additional adjustment for lagged deviance residuals.

For each outcome we used several models. First, we modelled the effect (incidence rate ratio; IRR) of each pollutant on hospital admissions on the same day (lag0) with pollutant concentrations set as continuous and then as categorical variables. From statistically significant effects observed at pollution levels above the exposure thresholds recommended by the World Health Organization (WHO) (9) we calculated population-attributable fractions. Second, we fitted the distributed-lag models, moving average models and cumulative effect models with lagged pollutant concentrations of up to seven days. The distributed-lag models were tested for multicollinearity. Third, the results from the distributed-lag models were disaggregated by gender and age. Finally, the initial lag0 models were stratified by the time of the year (April to September vs. October to March). Instead of using splines to model the time trend for these final stratified models we adjusted for categorical variables indicating the year and month.

Results were considered statistically significant at p<0.05 (two-tailed) and when the 95 % confidence interval of IRRs did not contain 1.00. All analyses were conducted with Stata/MP, version 17, 2021 (StataCorp LLC, College Station, TX, USA).

Over the 10-year period of interest, there were 98,567 hospitalisations for IHD, 41,327 for CI, and 46,643 for T2DM. Plotting daily case numbers against time did not reveal any obvious trend (Figure 1). On the other hand, air pollution levels gradually decreased over time (Figure 2). Expected seasonal patterns in individual pollutants were present, where all pollutants spiked in cold months, except for O₃, which followed the reverse pattern. Descriptive statistics for air pollutants and meteorological variables are shown in Table 1. Overall, there were few missing data on these variables.

Figure 1

Figure 2

Correlations of expected size and direction were found between the variables in the study (Table 2). Most air pollutants, except O₃, positively correlated with each other and inversely with ambient temperature. There were also positive correlations between admission counts for IHD, CI, and T2DM. We also saw some indication that those correlated with some air pollutants.

When air pollution concentrations were modelled as continuous variables, no significant risks for hospital admissions were found at lag0 (Table 3). Dichotomising the exposures at the short-term exposure cut-offs suggested by the WHO (9) yielded 3.90 % (95 % CI: 1.3 %, 6.6 %) higher risk of IHD when NO₂ exceeded 25 µg/m³ (Table 4). That translated into a population-attributable fraction of 2.60 % (95 % CI: 0.90 %, 4.28 %). A counterintuitive decreased risk of T2DM was observed with O₃ ≥60 µg/m³ and CO ≥4 mg/m³. There was some evidence of non-linearity in the risk associated with NO₂ – it became significantly increased when NO₂ exceeded 30 mg/m³ (compared to <15 mg/m³) and was in the ballpark of 5–9 % for IHD, CI, and T2DM (Table 5).

Modelling pollutant concentrations over multiple days revealed an increased risk of IHD associated with the 10 µg/m³ increase in NO₂ at lag0 (0.9 %; 95 % CI: 0.1 %, 1.8 %) and lag 2 (1.1 %; 95 % CI: 0.1 %, 2.1 %), and with 10 µg/m³ increase in O₃ at lag7 (1.5 %; 95 % CI: 0.5 %, 2.5 %) (Figure 3). The risk of CI increased with a 10 µg/m³ increase in PM₁₀ at lag5 (0.6 %; 95 % CI: 0.1 %, 1.0 %), SO₂ at lag4 (4.3%; 95 % CI: 0.8%, 8.0 %), and CO at lag5 (53 %; 95 % CI: 12 %, 110 %) (Figure 4). The risk of T2DM did not increase significantly with any of the exposure variables, although coefficients were borderline significant at some lags (Figure 5). CO exhibited unexpected patterns with decreased cumulative risk over lags 0–7.

Figure 3

Figure 4

Figure 5

Stratification by gender and age group showed that the risk of hospitalisation for IHD increased 1–2 % with increases in PM₁₀, PM_2.5, and NO₂ in people <65 yrs of age. O₃, on the other hand, was associated with an increased risk in those ≥65 yrs. at lags 4 and 7. The 7-day average concentrations and the cumulative effect of SO₂ over 7 days were associated with about 5 % higher risk. At lag3, the risk with CO was very high for women <65 yrs., while at other lags it was lower than 1.00 (Figure 6).

Figure 6

Hospital admissions for CI increased with higher PM₁₀ and PM_2.5 for women ≥65 yrs. at lag1, and women <65 yrs. at lag 2. SO₂ increased the risk for women <65 yrs. at lag4, and NO₂ at lag2. An increased risk in women <65 yrs. was present with O₃ at lag1, as well as with average O₃ concentrations over three and seven days. The cumulative risk with O₃ in men <65 yrs. was also significantly higher. At lag5, the risk with CO was increased for women <65 yrs., yet was below 1.00 at other lags (Figure 7).

The risk of hospital admissions for T2DM in men ≥65 yrs. increased 1–2 % with SO₂ at lag5, NO₂ at lag4, and O₃ at lag2. O₃ at lags 4 and 7 increased the risk in men <65 yrs. and women ≥65 yrs., respectively. At lag5, CO was associated with an increased risk in men ≥65 yrs., and a cumulative effect below 1.00 (Figure 8)

Table 6 shows seasonal differences in the observed effects. In warm months, NO₂ significantly increased the risk of hospital admissions for all three outcomes by 2–4 %, and PM₁₀ increased the risk of IHD admissions. Conversely, in cold months, coefficients consistently suggested lower risk of all admissions with higher air pollution levels.

Our analysis of 10 years worth of data on daily hospital admissions and air pollution in the city of Sofia has shown that, in certain scenarios, the risk of hospital admissions for IHD, CI, and T2DM increased between 1 and 5 % on average. With non-linear models this risk climbed to 5–9 % with higher NO₂. These effects took place within a week of high pollution days and were, in general, more common with gaseous pollutants (SO₂, NO₂, and O₃) than particulate matter. When we stratified the results by gender, we found that each pollutant had an effect on hospital admissions in some combination of gender and age subgroup. However, in some of these subgroups, the direction of the effect did not support the clinically based assumption that air pollution should increase the risk of hospital admissions. Unexpectedly, the risk of hospital admissions was higher in the warmer rather than colder months of the year, when air pollution is higher and when its effect on disease exacerbation may be expected to be stronger (24). On the individual level, the reason may be different seasonal behaviour: less time spent outdoors in the winter and more in the summer, which may add co-exposure to intense solar radiation and high temperatures. In addition, seasonal patterns in time spent outdoors, time spent in the city, and medication intake may differ across age groups (25), which we did not explore here. However, in an ecological study such as this it is not possible to rule out higher-level processes affecting data time patterns, such as exacerbation of multiple concomitant diseases in winter months affecting the primary diagnosis registered as a reason for hospitalisation.

Figure 7

Figure 8

Our findings generally agree with the wealth of evidence on air pollution and hospital or emergency department visits for cardiovascular diseases (10–13). However, we have not found a consistent pattern that would support our hypothesis, which too is not uncommon in literature. Although most studies confirm a correlation between hospitalisation for cardiovascular diseases and short-term peaks in air pollution, this finding is not universally true and studies vary in that respect (26). Some report that the risk increases on the same day (lag0 (23) and many other report a delay of several days (10–13).

As for T2DM, we observed a less clear association than for IHD and CI, present only in some age and gender subgroups. In contrast, a recent systematic review reported an elevated risk of T2DM with higher air pollution (27).

We also observed some counterintuitive behaviour for CO at several lags, in line with studies reporting strange patterns with CO at specific lags (28). However, our multicollinearity test suggests that there is no reason for concern. Considering that a large meta-analysis reported that the risk of myocardial infarction associated with a 1 mg/m³ increase in CO is 1.052 (95 % CI: 1.017, 1.089) (29), we believe that the large shifts in our findings for CO, regardless of their direction, are influenced by unaccounted confounding factors or correlations and are therefore not trustworthy.

Demographic stratification also did not reveal a clear-cut pattern, as the hospitalisation increased in both older and younger individuals for different combinations of a pollutant, outcome, and lag. Other authors suggest that the elderly are more susceptible to the effects of air pollution (23). However, recent WHO reviews on short-term effects of different air pollutants on hospital admissions, emergency room visits, and mortality reported much higher diversity for this population group (29–31).

To our knowledge, this is the first study in Bulgaria to model short-term associations between air pollution and IHD, CI, and T2DM hospitalisations that covers such a long period of ten years. In their reports of short-term associations between air pollution and cardiovascular diseases, earlier studies in the country relied on bivariate analyses, which are not ideally suited for modelling count data and which ignore time trends and autocorrelations in time series (19–21).

This study too is not without its limitations. First, as an ecological study it is not particularly informative about the risk of disease exacerbation on the patient level. Aggregated time-series data do offer some insight into the processes and drivers of the use of healthcare services on a population level, but this approach does not address the issue of individual prevention.

Second, because we covered a long period of time, we had to rely on air quality data from only one air quality monitoring station in Sofia, which may not fully reflect local differences in pollutant levels. Other government-owned stations did not measure PM_2.5 or focused only on concentrations near traffic sites. On the other hand, municipal stations in Sofia were installed only recently and are less precise and do not make part of the official national air quality network managed by the Executive Environment Agency under the administration of the Minister of Environment and Water. Moreover, their sensors become less accurate when air humidity level is high, which is a limitation of that sensor technology.

Our study also does not include outdoor measurements in the residential areas of patients, which would allow future research to extrapolate their exposure more accurately. This would require ethical approval and installation of monitoring stations whose air pollution measurements would then be adjusted for the distance to each residence address or modelled in another more sophisticated fashion. However, this was beyond the scope of the present study.

Third, we only had access to data on the total number of hospital admissions per ICD code and lacked information whether the hospitalisation was a repeated or first registered clinical event. The use of a given ICD code does not always reflect the primary cause that has led to hospitalisation. Moreover, we could not distinguish between fatal and non-fatal cases, which is important, since short-term exposure to air pollution increases the risk of mortality (31).

To overcome some of these shortcomings, we intend to extend the analyses, provided the data can be obtained, to other major cities in Bulgaria, and pool the results across the cities. Future research should attempt to extract information on patient co-morbidities, medical history, as well as the outcome of their hospitalisation. It would also be informative to include data of emergency department visits and calls, not hospital admissions alone.