Open Access
Open Peer Review

This article has Open Peer Review reports available.

How does Open Peer Review work?

Temperature-related mortality estimates after accounting for the cumulative effects of air pollution in an urban area

  • Svetlana Stanišić Stojić1Email author,
  • Nemanja Stanišić2 and
  • Andreja Stojić3
Environmental Health201615:73

Received: 22 February 2016

Accepted: 30 June 2016

Published: 11 July 2016



To propose a new method for including the cumulative mid-term effects of air pollution in the traditional Poisson regression model and compare the temperature-related mortality risk estimates, before and after including air pollution data.


The analysis comprised a total of 56,920 residents aged 65 years or older who died from circulatory and respiratory diseases in Belgrade, Serbia, and daily mean PM10, NO2, SO2 and soot concentrations obtained for the period 2009–2014. After accounting for the cumulative effects of air pollutants, the risk associated with cold temperatures was significantly lower and the overall temperature-attributable risk decreased from 8.80 to 3.00 %. Furthermore, the optimum range of temperature, within which no excess temperature-related mortality is expected to occur, was very broad, between −5 and 21 °C, which differs from the previous findings that most of the attributable deaths were associated with mild temperatures.


These results suggest that, in polluted areas of developing countries, most of the mortality risk, previously attributed to cold temperatures, can be explained by the mid-term effects of air pollution. The results also showed that the estimated relative importance of PM10 was the smallest of four examined pollutant species, and thus, including PM10 data only is clearly not the most effective way to control for the effects of air pollution.


Air pollutionTemperatureMortalityEnvironmental exposureDeveloping countries


Environmental stressors, such as extreme temperature events and air pollution, pose a significant challenge to human societies, particularly to the growing urban population worldwide [1]. While most studies investigating the adverse health effects of air pollution have adjusted for temperature as a confounder, accounting for the effects of air pollution in the studies aimed at assessing the relationship between temperature and mortality has been less common [2]. The relatively few studies considering the interactive effects of meteorological variables and air pollution on daily mortality have reported inconsistent results [3, 4]. Significant inconsistency in the effect estimates across studies is at least partly associated with the methodological differences in exposure assessment and confounder control [5]. Furthermore, the majority of studies assessed air pollution-related mortality by including lag-specific effects of pollutants, despite the fact that detrimental impact of air pollution is not limited to few days preceding adverse health outcome, particularly in highly polluted areas of developing countries.

Better understanding of the temperature-related mortality is urgently needed, when considering the upcoming atmospheric changes and projected global temperature rise of 1.8 to 4.0 °C by 2100 [6]. The relationships between temperature and mortality were reported to be linear or non-linear, with death rates increasing in either direction from an optimum value [7, 8]. However, the relative importance of cold spells and heat waves still remains an issue of scientific debate. The latest study of Gasparrini et al. [9] including more than 74 million deaths recorded in 13 countries and 384 locations showed that cold-related deaths outnumbered heat-related deaths by a factor of nearly 20, and that most of the temperature-related mortality burden was attributable to milder, non-optimum weather, and not to extreme temperature events. Nevertheless, it should be emphasized that factors such as susceptibility or resilience, including air pollutants, in the most of cases have not been included in the analysis [10].

The evidence on the effects of environmental stressors in developing countries is relatively scarce because of poor environmental monitoring and a paucity of health surveillance data [11]. In this study, we compared the relative contributions from heat and cold to circulatory and respiratory mortality burden, before and after accounting for the cumulative effects of air pollution, as estimated by the method proposed herein. The study is based on data from Belgrade, Serbia, which is among the most polluted European cities [12, 13] and has recorded very high cardiovascular death rates over the last 15 years, ranging from 635 to 677 per 100,000 inhabitants [14, 15].


The mortality data were stratified into age categories, and only people aged 65 and older were included in the analysis in order to increase the statistical power. Daily data of circulatory and respiratory mortality (cause-of-death coding I00-I99 and J00-J99 according to the International Classification of Diseases 10th revision, ICD-10 code) for the period from 2009–2014 were obtained from the Institute of Public Health Belgrade. Temperature data was obtained from the Global Data Assimilation System with spatial resolution of 1° and combined over 24-h periods to provide mean values.

In order to examine the associations between environmental stressors and mortality, two separate models were specified. Namely, the first model included temperature and other factors affecting mortality, such as season and day of the week, while the second, along with the above-mentioned factors, included air pollutants: PM10, NO2, SO2 and soot. The daily mean pollutant concentrations of PM10, NO2, SO2 and soot were obtained from 7, 21, 22 and 16 monitoring stations, respectively, uniformly distributed throughout the city area, within the monitoring network of the Institute of Public Health Belgrade. Measurements were conducted at automatic stations using a beta-ray attenuation sampler (Thermo FH 62 I-R) for PM10 and referent sampling devices, Horiba APNA 360 and APSA 360 analysers, for NO2 and SO2, respectively. At semi-automatic measuring stations, PM10 concentrations were determined using a referent Sven Leckel sampler, whereas determination of NO2 and SO2 concentrations was performed based on ISO 6768:1998 (Modified Griess-Saltzman method) and ISO 6767:1990 (tetrachloromercurate (TCM)/pararosaniline method) standards, respectively. Soot concentrations were obtained by the use of a Pro-Ekos device based on ISO 9835:1993 standard.

The effect of temperature on mortality is shown to be both delayed in time and shifted by a so-called harvesting effect or mortality displacement – a temporal advance in deaths occurring in the frailest population, typically followed by a period of reduced mortality [6]. In order to simultaneously capture the non-linear exposure-response dependencies and determine the time period between the exposure and health outcome, we used distributed lag non-linear models (DLNM). This methodology is based on the concept of cross-basis, a bidimensional space of functions that captures simultaneously the shape of the relationship between a predictor and a dependent variable along the space of the predictor and the lag dimension of its occurrence [16]. Implementation of the DLNM modeling framework [7] is available in the statistical environment R [17]. Based on the methodology presented in Gasparrini et al. [9], we modelled the exposure-response curve using a quadratic B-spline with three internal knots placed at the 10th, 75th, and 90th percentiles of the observed temperature distribution, and the lag-response curve (with maximum lag up to 21 days) using a natural cubic B-spline with an intercept and three internal knots equally spaced on the log scale. In order to control for seasonal and long-term effects, we included a natural cubic B-spline of time with 8° of freedom per year. Finally, the model also included a categorical variable indicating the day of the week.

The quasi-Poisson regression model with no pollutants included is specified as follows:
$$ E\left({Y}_t\right)= \exp \left\{\alpha +\beta {T}_{t,l}+NS\left( time,\ df\right)+\lambda DO{W}_t\right\} $$
where E(Y t ) is expected daily death count on day t; α is the intercept; T t,l is a matrix of variables obtained by the transformation of mean daily temperature, β is a vector of coefficients for T t,l , and l is the lag day; NS(time,df) is the natural cubic spline of time; and DOW t and λ are dummy variables representing the day of the week and the corresponding vector of coefficients, respectively. Quasi-Poisson was used to compensate for overdispersion since the residual deviance was larger than the residual degrees of freedom (the dispersion parameter was 1.06). We shall refer to this model in further text as the base model.

Results and discussion

The analysis included a total of 56,920 residents aged 65 years or older who died from circulatory (93.27 %) and respiratory diseases (6.73 %). As shown in Fig. 1, mortality annual variations for the entire period exhibited a strong seasonal pattern, with peak values in the middle of February and minimum values around September, 1st. Two minor increases occurred at the beginning of July and the end of October, suggesting the potential impact of sudden temperature changes on human health. The figures presenting circulatory and respiratory daily mortality segmented by age and gender, and their seasonal variations are shown in Additional files 1 and 2.
Fig. 1

Smoothed empirical mortality annual variations for the period 2009–2014. Mortality variations exhibit seasonal pattern, with peak values in the middle of February and minimum values around the start of September

The climate in Belgrade is moderate continental with all four seasons. During the study period, mean daily temperatures were in the range from −17.4 to 30.1 °C, with the average value of 11.9 °C. The average number of days with temperatures below 0 °C was 40 per year, whereas the coldest and warmest months were January and July, with 0.4 and 22.2 °C, respectively. The average seasonal temperatures for the entire period were 12.1, 21.4, 12.7 and 1.1 °C for the spring, summer, autumn and winter, respectively.

The daily mean concentrations of PM10, NO2, SO2 and soot for the entire period were 48.3, 32.7, 15.5 and 20.7 μg m−3, respectively. The number of days with average PM10 concentrations exceeding 50 μg m−3 was in the range from 75 to 155 per year, which is considerably higher than the Air Quality Standard margin (35 exceedances per year), whereas the mean annual NO2 concentrations did not exceed the value of 40 μg m−3, nor any mean daily SO2 levels higher than the recommended limit of 135 μg m−3 were observed.

The levels of PM10 and SO2 exhibited a clear seasonal pattern, with highest concentrations observed during cold season (Additional file 3), whereas the seasonality of soot was less pronounced. Furthermore, the moving average PM10 and SO2 concentrations were highly correlated, suggesting that the variations of two pollutants, particularly those long-term, were driven by the same anthropogenic and meteorological factors. Conversely, NO2 concentrations exhibited pronounced weekly, but not seasonal variations, which can be partly attributed to the traffic-related emissions that are present throughout the year.

The lag-specific effects of temperature on mortality obtained from the base model are provided in Fig. 2. As can be seen, the heat-related mortality appears to be an acute event followed by a reduction in death rates, with maximum effects observed on the same day or lags of 1–3 days. Nevertheless, the effects of cold spells were observed to be more evenly distributed across 3–6 days of exposure with less evidence of subsequent harvesting, which complies with previous studies reporting sustained health effects of low temperatures [18].
Fig. 2

Plot of the exposure-lag-response relationship between temperature and all-cause mortality, with reference at 21 °C. The effects of heat are observed on the same day or after 1–3 days, whereas the effects of cold were distributed across few days of exposure

As regards the cumulative effects of temperature over the entire lag period, the U-shaped exposure-response curve was identified as the most representative. The minimum relative risk corresponding to the optimum temperature of 21 °C is presented in Fig. 3 as a dotted line, whereas the cut-off values at the 1st and 99th percentile of the temperature distribution, which stand for extreme cold and heat, are displayed as dashed vertical lines.
Fig. 3

Overall cumulative exposure-response curve and temperature distribution in Belgrade 2009–2014. Dotted and dashed lines represent optimum and extreme temperatures, respectively

According to the findings of Gasparrini et al. [9], cold temperatures were responsible for advancing a substantial fraction of deaths (7.29 %), compared to hot days (0.42 %), whereas most deaths were attributable to exposure to moderate temperatures. As shown in Fig. 3, the highest relative risk was observed for a small proportion of extremely cold days, whereas the estimated effect of the days with moderate temperatures, corresponding to the smooth response, was also noticeable.

Upon estimating the effects of temperature on mortality, we proceeded to examine the effects of air pollutants. The strong positive correlations between PM10, SO2 and soot (Table 1) suggest that these species share the same origin, and thus, the effects of some pollutants may distort the observable patterns of others.
Table 1

Pollutant correlations for the entire period, and for cold and warm season











































Pollutant correlations for the entire period (left), and for cold and warm season (right, lower and upper triangular, respectively)

Given that the extreme cold events and increased levels of pollutants often coincide, it is necessary to understand the combined effects of these two factors on mortality. We have developed our modeling strategy relying on the following logical sequence:
  1. 1.

    Increased concentrations of air pollutants (at least within their normal range of concentrations) require a certain amount of time to reach their cumulative effect on mortality, the time being not less than 24 h, and thus, daily average is the level of aggregation of the measurement data used in this and numerous other studies. This cumulative effect is expected to be more stable, and thus, statistically robust in comparison to short-term lag-specific effects (referring to lags of up to 7 days).

  2. 2.

    Up to the nth period (day), the longer the high values of pollutants persist, the greater the cumulative effect on mortality. As the cumulative effect on mortality reaches its peak (period n), it largely depletes the pool of those susceptible in the observed population, as their death is being brought forward by a matter of a few days.

  3. 3.

    Owing to the aforementioned harvesting effect, the cumulative primary effect of the pollutant starts to decrease afterward, reaching its minimum value on day n1.

Based on these suppositions, two variables were used to describe the temporal distribution of the cumulative effect of each pollutant. The first variable represented the average pollutant concentrations during the specific timeframe at the end of which (day n) the primary effect of the pollutant on mortality reaches its local maximum. The second variable is defined to capture the temporal positioning of the minimum rate of mortality (day n1) attributable to the post-harvesting reduction in death rates. Accordingly, the cumulative effects of each pollutant are included as two simple moving averages (SMA) of its observed concentrations in 10 μg m−3 (P) over n consecutive days, ending on-and inclusive of-day t.
$$ SMA{(P)}_{t,n}=\frac{1}{n}{\displaystyle \sum_{i=1}^n}{P}_{\left(t-i\right)+1} $$

Rather than assume any particular temporal distribution of the effect, as is typically done in the literature, we visually inspected the values of regression coefficient estimates that relate SMA concentrations of each of the four air pollutants to mortality rates within a reasonable timeframe of 60 days, seeking for “peak and trough” temporal patterns that can be attributed to the corresponding harvesting effects of the specific pollutant. The temporal pattern was expected to look as follows: the relative risk should start off at a value close to 1 and stay relatively close to 1 for short time periods (no major effect is expected within the first few days), then steadily increase to a finite value within a reasonable timeframe, and lastly, display a steady decline to a value below 1, with the decline roughly matching the timing and magnitude of the preceding increase. The duration period of the effect depends on the specific system characteristics, such as population susceptibility, pollutant concentrations and their spatio-temporal dynamics.

As certain pattern deviations were expected due to strong correlations (referring to SMA values and not daily observations), we included the pollutants in sequence iteratively, giving priority to those with stronger effects (as estimated by the assessed relative risk) and clearer temporal patterns of these effects. Figure 4 shows the procedure for determining the order of pollutant inclusion in the final model. Lags 0–60 refer to timeframes for which the relative risk was estimated by including one by one timeframe-specific SMA pollutant concentration in the model. As can be seen, once the clear harvesting effect pattern had been identified and attributed to a specific pollutant, the pollutant was added to the final model, and the estimation of the relative risk was repeated for the remaining pollutants.
Fig. 4

The procedure for including air pollutants in the final model. The estimation of the timeframe-specific relative risk was repeated each time after exclusion of the pollutant with the most clear harvesting effect pattern. Dotted lines refer to timeframes that are associated with highest positive or negative relative risk

As can be seen, the effects of pollutants were relatively stable between iterations due to the fact that SMA were used instead of lag values. The final model is specified as follows:
$$ \begin{array}{l}E\left({Y}_t\right)= \exp \left\{\alpha +\beta {T}_{t,l}+NS\left( time,\ df\right)+\lambda DO{W}_t+{\gamma}_1SMA{\left({\mathrm{SO}}_2\right)}_{t,19}+{\gamma}_2SMA{\left({\mathrm{SO}}_2\right)}_{t,47}\right.\hfill \\ {}\kern10.5em +{\gamma}_3SMA{\left(\mathrm{Soot}\right)}_{t,28}+{\gamma}_4SMA{\left(\mathrm{Soot}\right)}_{t,51}+{\gamma}_5SMA{\left({\mathrm{NO}}_2\right)}_{t,53}+{\gamma}_6SMA{\left({\mathrm{NO}}_2\right)}_{t,21}\hfill \\ {}\left.\kern10.5em +{\gamma}_7SMA{\left({\mathrm{PM}}_{10}\right)}_{t,3}+{\gamma}_8SMA{\left({\mathrm{PM}}_{10}\right)}_{t,10}\right\}\hfill \end{array} $$

The detailed model estimates are presented in Additional file 4, in the columns labeled “Saturated model”.

The risk of excess death with short-term exposure to elevated concentrations of PM10 was observed, whereas for SO2 and soot, the marked effects were found within 2–3 weeks of exposure. A 10 μg m−3 increase in PM10 concentrations is followed by a negligible increase in mortality (0.04 %), whereas a moderate increase in death rates in the range from 0.7 to 1.3 % was associated with the same increment in NO2, SO2 and soot concentrations. Comparable short-term health effects of exposure to PM were also reported in the extensive review of Brook et al. [19]. They concluded that the elevated PM2.5 concentrations over 5 days lead to the increase in circulatory mortality risk by 0.4 % to 1.0 %. Thereby, it is important to mention that the health effects of PM2.5 are more severe compared to PM10, because the particles of smaller diameter penetrate deeper in the respiratory tract. The air pollution-related mortality estimates also comply with the findings of other extensive studies [20, 21], reporting moderate increases in mortality associated with the exposure to PM, SO2 and NO2. As can be seen in Fig. 4, the effect pattern was inverted in the case of NO2, which could indicate that the exposure to high NO2 levels reaches their maximum effect (day n) beyond the examined time period. This complies with the results of our previous study [13], aimed at investigating the association between short- and medium-term (up to 90 days) exposure to nitrogen dioxide (NO2) and mortality within the several timeframes. As shown, the short-term exposure to NO2 exhibited negative associations with death rates, whereas the medium-term exposure to NO2 was associated with significant increase in mortality.

Furthermore, it could be assumed that the effects attributed to NO2 exposure are not realistic, but rather a result of omitted variable bias, which is even more likely taking into account its considerable relative importance (Fig. 6), which was highest of all species. Likewise, Brook et al. [22] reported surprisingly strong relationship between NO2 and non-accidental mortality across 10 cities in Canada. They assumed that NO2 concentrations were acting as an indicator of some other variable or process, because NO2 levels exhibit high correlations in time and space with a range of toxic traffic-related species, including volatile organic compounds and polycyclic aromatic hydrocarbons, which could also contribute to the observed effects.

In Fig. 5, the overall effects of temperature on mortality before and after accounting for air pollutant concentrations are compared.
Fig. 5

Overall cumulative exposure-response curve after accounting for the effects of air pollution and temperature distribution. After adjusting for air pollution, the optimum temperature range was relatively broad

After accounting for the effects of air pollutants, estimates of backward attributable risk of non-optimal temperature decreased from 8.80 % (5,009 people) to 3.00 % (1,706 people) (for detailed information on calculation of attributable risk see [23]). More specifically, the base model results showed that the exposure to cold temperatures was associated with significantly higher risk (7.35 %), compared to hot days (1.50 %). However, after adjusting for air pollutants, the relative risk associated with cold temperatures was lowered to 1.06 %, whereas the risk of heat-related mortality slightly increased to 1.96 %. When adopting a backward perspective, the risks are simultaneously computed for the same time t, and thus, negligible difference (\( b-A{F}_{x,t}^{r_{1+}{r}_2}\le\ b-{A}_{x,t}^{r_1} + b-{A}_{x,t}^{r_2} \)) could be observed between the sum of components (1.06 % and 1.96 %) and the overall attributable risk (3.00 %) [23]. It is also worth noting that after including air pollutants in the model, the optimum temperature range was relatively broad, from −5 to 21 °C, which differs from previous findings that most of related deaths were caused by moderate temperatures [9]. Correspondingly, the relative contributions of extreme temperature events have also changed (Table 2).
Table 2

Relative risk for the 1 and 99th percentile of the observed temperature distribution (95 % CI)


Before accounting for the effects of air pollutants

After accounting for the effects of air pollutants

−10.0 °C

1.40 (1.12, 1.76)

1.14 (0.86, 1.51)

27.0 °C

1.32 (1.10, 1.60)

1.45 (1.20, 1.76)

Anderson and Bell [24] showed that the temperature-related mortality risk was slightly lowered (approx. 0.4 % on average) after controlling for the effects of O3 and PM10. As shown in Additional file 5, the change in mortality risk after accounting for the short-term lag-effects of air pollutants was negligible. However, after accounting for the cumulative effects of air pollutants, as proposed herein, the significant change in estimated risk was observed.

As part of a supplementary analysis, and with the aim of arriving at a more parsimonious model specification, we performed optimization with respect to information criteria, and investigated whether a model with a double-threshold parameterization can be used as a reasonable approximation to the saturated model. Firstly, we examined whether some of the model terms could be left out without substantial loss to the model fit. In order to gauge the relative quality of the model in terms of the goodness-of-fit/complexity trade-off, we used a quasi-likelihood (we used the C-hat [the dispersion parameter] value obtained from the saturated model, which is 1.06) counterpart of Akaike’s information criterion (QAICc), specifically its implementation available in the GLMULTI package in R [25]. An exhaustive screening, which included each of the 2,048 (211) possible models for the specified set of variables, showed that the best model, in terms of the information criteria, is the following:
$$ \begin{array}{l}E\left({Y}_t\right)= \exp \left\{\alpha +\beta {T}_{t,l}+NS\left( time,\ df\right)+{\gamma}_1SMA{\left({\mathrm{SO}}_2\right)}_{t,19}+{\gamma}_2SMA{\left({\mathrm{SO}}_2\right)}_{t,47}+{\gamma}_3SMA{\left(\mathrm{Soot}\right)}_{t,28}\right.\hfill \\ {}\left.\kern10.5em +{\gamma}_4SMA{\left(\mathrm{Soot}\right)}_{t,51}+{\gamma}_5SMA{\left({\mathrm{NO}}_2\right)}_{t,53}+{\gamma}_6SMA{\left({\mathrm{NO}}_2\right)}_{t,21}\right\}\hfill \end{array} $$

Detailed model estimates are presented in Additional file 4, “Optimal model” columns.

The screening process also enabled an insight into the estimated importance of each of the variables. The estimated importance was calculated as the sum of the relative evidence weights of all models in which the term appears. Relative evidence weights are computed as exp(ΔIC/2), where ΔIC is the difference in IC between a model and the best model, and they are normalized so that they sum up to one. They can be interpreted as probabilities for each model to be the best in the set. Given the fact that PM10 concentrations have the smallest estimated relative importance of all examined species (Fig. 6), including only PM10 in the form of a 0–1 lag variable, as performed in research [26] is clearly not the most effective way to control for the effects of air pollution.
Fig. 6

Model-averaged relevance of terms. Dotted line represents threshold importance

Since the saturated model indicates the existence of an optimum temperature range, it is valuable to know whether a model with the double-threshold parameterization (with −5 to 21 °C being the thresholds) could be used as a simple approximation to it, without significant loss to the model fit. The analysis showed, however, that the loss in the model fit, from 0.2832 for saturated to 0.2574 for double-threshold model as measured by correlation between the predicted and the observed values (pseudo R squared), was not compensated for in terms of information criteria—QAICc increased from 12518.25 to 12553.38. Visual inspection indicated that the linear components of the double-threshold model are incapable of accurately describing the progressive nature of the effects of extreme temperatures.

In comparison to the herein proposed SMA method, DLNM framework is technically more sophisticated and more appropriate for the exploration of the particular lag-specific and exposure-level-specific pollutant effects. It is also superior when long-time series are available. However, given that DLNM framework does not employ the concept of a stable cumulative effects between pollution and mortality, but rather tries to describe the entire exposure-lag-response association, it is demanding regarding the number of parameters that need to be estimated. For instance, using DLNM for describing the effect of temperature on mortality would involve estimating a total of 5 × 5 = 25 parameters.1 Furthermore, describing the effects of four air pollutants would require an additional 100 parameters on top of the 80 parameters already included in the model (25 parameters are estimated for describing the effect of temperature, 8 × 6 years = 48 parameters for capturing seasonal and long-term trend, 6 parameters for day of the week, and 1 overall intercept). Conversely, the herein proposed procedure requires 8 + 8 = 16 parameters.

Furthermore, the herein proposed SMA method is the most parsimonious method that is consistent with the assumption that significant cumulative, medium-term effects of air pollutant levels on mortality exist. If this assumption is valid, which has been confirmed by previous studies, then it suffices to include only two properly selected SMA terms for each pollutant that captures the majority of their respective effects. The choice of SMA lags describes the temporal positioning of the cumulative effect’s peak, while the regression estimates of the two corresponding regression coefficients describe the magnitude of the effect for each particular air pollutant. Hence, in the case of short-time series, when air pollutants are included into the research model solely to account for their effect on mortality, the SMA approach may be comparatively more efficient than DLNM.

The results of using the DLNM to model the effects of four pollutants are presented in Additional file 6. As can be seen, the observed shapes of the effects confirm the central hypothesis of our study, which is that, besides the well-evidenced short-term effects of air pollution on mortality (most visible in the case of soot), significant cumulative (perhaps even seasonal) medium-term effects of air pollutant levels on mortality (most visibly in the case of SO2) are observed. The harvesting effect pattern is also clearly visible in both short- and medium-term effects. In Additional file 7, temperature-related mortality risk is presented, after accounting for air pollutants using SMA and DLNM.

Some general limitations of retrospective epidemiological studies aimed at health risk estimation are also applicable to this study. By analyzing the effects of air pollution and temperature, the results obtained from a number of central monitoring stations are assumed to be representative of individual exposure, although people have different activity patterns and spend a significant period of time indoors. Furthermore, certain species are selected as indicators of air pollution, although pollution refers to a complex mixture of gaseous species and particles that permanently interact and undergo physico-chemical transformations [27]. However, regardless of the potential limitations, the results of retrospective epidemiological studies based on large-sample statistics remain an important source of information for designing environmental regulations [28]. The advantage of the herein proposed method is that it includes the cumulative mid-term effects of air pollution in the traditional Poisson regression model in a parsimonious way, which complies with the nature of their impact. The limitation of the proposed method is possible overfitting due to the limited number of observations, which needs to be eliminated through further studies.


The herein proposed method is based on the theoretical premise that there are significant cumulative mid-term effects of air pollution on mortality, which are more stable, and thus more statistically robust, than the lag-specific effects commonly included in regression models. The proposed method requires the researcher to visually identify a “peak and trough” temporal pattern for the effects of each of the pollutants, which is then included in the model using only two parameters, yielding a rather parsimonious model specification.

The presented results suggest that the inference that cold-attributable mortality currently accounts for one order of magnitude more deaths than mortality associated with heatwaves might not be globally applicable. Accounted for in an appropriate way, the effects of air pollutants explained most of the mortality previously associated with cold temperatures, and a greater temperature-attributable mortality burden was associated with the contribution of heat.

The more realistic estimates of the mortality risk associated with extreme temperature events becomes increasingly important for planning of future public health interventions and adaptation measures, with the global warming impact that is already underway. However, further research is needed to identify all environmental risk factors and clarify their complex associations, which have been largely obscured by the current approach.


PM10, particles smaller than 10 μm; NO2, nitrogen dioxide; SO2, sulfur dioxide; CI, confidence interval


We used a quadratic B-spline with three internal knots and no intercept for the exposure-response curve: 3 df for knots + 2 df for the degree of the piecewise polynomial (quadratic) + 0 df for intercept (no intercept) = 5 df, and a natural cubic B-spline with three internal knots and an intercept for the lag-response curve: 3 df for knots + 1 df for intercept + 1 df = 5 df.




The authors would like to express their deepest gratitude to the Institute of Public Health Belgrade for technical support and additional analysis of mortality rates in Belgrade.


This study was conducted as part of the projects No III43007 and No III41011, which were financed by the Ministry of Education, Science and Technological Development of the Republic of Serbia within the framework of integrated and interdisciplinary research for the period 2011–2016.

Availability of data and materials

Raw data will not be shared because the authors are not authorized for distribution of data.

Authors’ contributions

SSS and AS conceived the study design. The analysis was conducted by NS. The paper was jointly written by SSS, AS and NS. All authors approved the final draft.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Ethics approval and consent to participate

Not applicable.

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Authors’ Affiliations

Faculty of Physical Chemistry, University of Belgrade
Singidunum University
Institute of Physics Belgrade, University of Belgrade


  1. Stieb DM, Judek S, Burnett RT. Meta-analysis of time-series studies of air pollution and mortality: effects of gases and particles and the influence of cause of death, age, and season. J Air Waste Manage Assoc. 2002;52:470–84.View ArticleGoogle Scholar
  2. Analitis A, Michelozzi P, D’Ippoliti D, de Donato F, Menne B, Matthies F, Atkinson RW, Iñiguez C, Basagaña X, Schneider A, Lefranc A. Effects of heat waves on mortality: effect modification and confounding by air pollutants. Epidemiology. 2014;25:15–22.View ArticleGoogle Scholar
  3. Pattenden S, Armstrong B, Milojevic A, Heal MR, Chalabi Z, Doherty R, Barratt B, Sari Kovats R, Wilkinson P. Ozone, heat and mortality: acute effects in 15 British conurbations. Occup Environ Med. 2010;67:699–707.View ArticleGoogle Scholar
  4. Burkart K, Canário P, Breitner S, Schneider A, Scherber K, Andrade H, Alcoforado MJ, Endlicher W. Interactive short-term effects of equivalent temperature and air pollution on human mortality in Berlin and Lisbon. Environ Pollut. 2013;183:54–63.View ArticleGoogle Scholar
  5. Hoek G, Krishnan RM, Beelen R, Peters A, Ostro B, Brunekreef B, Kaufman JD. Long-term air pollution exposure and cardio-respiratory mortality: a review. Environ Health. 2013;12:43–59.View ArticleGoogle Scholar
  6. Ye X, Wolff R, Yu W, Vaneckova P, Pan X, Tong S. Ambient temperature and morbidity: a review of epidemiological evidence. Environ Health Persp. 2012;120:19–28.View ArticleGoogle Scholar
  7. Gasparrini A. Distributed lag linear and non-linear models in R: the package dlnm. J Stat Softw. 2011;43:1–20.View ArticleGoogle Scholar
  8. Pinheiro SDLLD, Saldiva PHN, Schwartz J, Zanobetti A. Isolated and synergistic effe17cts of PM10 and average temperature on cardiovascular and respiratory mortality. Rev Saude Publ. 2014;48:881–8.View ArticleGoogle Scholar
  9. Gasparrini A, Guo Y, Hashizume M, Lavigne E, Zanobetti A, Schwartz J, Tobias A, Tong S, Rocklöv J, Forsberg B, Leone M, De Sario M, Bell ML, Guo Y-LL, Wu C-f, Kan H, Yi S-M, Coelho M de SZS, Saldiva PHN, Honda Y, Kim H, Armstrong B. Mortality risk attributable to high and low ambient temperature: a multicountry observational study. Lancet. 2015;386:369–75.View ArticleGoogle Scholar
  10. Dear K, Wang Z. Climate and health: mortality attributable to heat and cold. Lancet. 2015;386:320–2.View ArticleGoogle Scholar
  11. Goldberg MS, Burnett RT, Valois M-F, Flegel K, Bailar III JC, Brook J, Vincent R, Radon K. Associations between ambient air pollution and daily mortality among persons with congestive heart failure. Environ Res. 2003;91:8–20.View ArticleGoogle Scholar
  12. Stojić A, Stojić SS, Reljin I, Čabarkapa M, Šoštarić A, Perišić M, Mijić Z . Comprehensive analysis of PM10 in Belgrade urban area on the basis of long-term measurements. Environ Sci Pollut R. 2016;23:10722–32.View ArticleGoogle Scholar
  13. Stanišić Stojić S, Stanišić N, Stojić A, Šoštarić A. Single and combined effects of air pollutants on circulatory and respiratory system-related mortality in Belgrade, Serbia. J Toxicol Env Heal A. 2015; doi: 10.1080/15287394.2015.1101407
  14. Statistical Office of the Republic of Serbia. Vital statistics - the additionally processed mortality data. 2015.Google Scholar
  15. World Health Organization, Regional Office for Europe. Detailed Mortality Database. 2015. Accessed 27 July 2015.
  16. Gasparrini A, Armstrong B, Kenward MG. Distributed lag non-linear models. Stat Med. 2010;29:2224–34.View ArticleGoogle Scholar
  17. R Core Team. R. A language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2014. Accessed 28 Dec 2015.Google Scholar
  18. Breitner S, Wolf K, Devlin RB, Diaz-Sanchez D, Peters A, Schneider A. Short-term effects of air temperature on mortality and effect modification by air pollution in three cities of Bavaria, Germany: a time-series analysis. Sci Tot Environ. 2014;485:49–61.View ArticleGoogle Scholar
  19. Brook RD, Rajagopalan S, Pope CA, Brook JR, Bhatnagar A, Diez-Roux AV, Holguin F, Hong Y, Luepker RV, Mittleman MA, Peters A. Particulate matter air pollution and cardiovascular disease an update to the scientific statement from the American Heart Association. Circulation. 2010;121:2331–78.View ArticleGoogle Scholar
  20. Goldberg MS, Burnett RT, Stieb DM, Brophy JM, Daskalopoulou SS, Valois MF, Brook JR. Associations between ambient air pollution and daily mortality among elderly persons in Montreal, Quebec. Sci Tot Environ. 2013;463:931–42.View ArticleGoogle Scholar
  21. Shah AS, Langrish JP, Nair H, McAllister DA, Hunter AL, Donaldson K, Newby DE, Mills NL. Global association of air pollution and heart failure: a systematic review and meta-analysis. Lancet. 2013;382:1039–48.View ArticleGoogle Scholar
  22. Brook JR, Burnett RT, Dann TF, Cakmak S, Goldberg MS, Fan X, Wheeler AJ. Further interpretation of the acute effect of nitrogen dioxide observed in Canadian time-series studies. J Expo Sci Environ Epidemiol. 2007;17:S36–44.View ArticleGoogle Scholar
  23. Gasparrini A, Leone M. Attributable risk from distributed lag models. BMC Med Res Methodol. 2014;14:55–63.View ArticleGoogle Scholar
  24. Anderson BG, Bell ML. Weather-related mortality: how heat, cold, and heat waves affect mortality in the United States. Epidemiology. 2009;20:205–13.View ArticleGoogle Scholar
  25. Calcagno V, de Mazancourt C. glmulti: an R package for easy automated model selection with (generalized) linear models. J Stat Softw. 2010;34:1–29.View ArticleGoogle Scholar
  26. Lin YK, Ho TJ, Wang YC. Mortality risk associated with temperature and prolonged temperature extremes in elderly populations in Taiwan. Environ Res. 2011;111:1156–63.View ArticleGoogle Scholar
  27. Brauer M, Amann M, Burnett RT, Cohen A, Dentener F, Ezzati M, Henderson SB, Krzyzanowski M, Martin RV, Van Dingenen R, van Donkelaar A. Exposure assessment for estimation of the global burden of disease attributable to outdoor air pollution. Environ Sci Technol. 2012;46:652–60.View ArticleGoogle Scholar
  28. Solomon PA, Costantini M, Grahame TJ, Gerlofs-Nijland ME, Cassee FR, Russell AG, Brook JR, Hopke PK, Hidy G, Phalen RF, Saldiva P, Sarnet SE, Balmes JR, Tager IB, Özkaynak H, Vedal S, Wierman SSG, Costa DL. Air Qual Atmos Health. 2012;5:9–62.View ArticleGoogle Scholar


© The Author(s). 2016