Short term association between ozone and mortality: global two stage time series study in 406 locations in 20 countries

1Ground level ozone is a highly reactive, oxidative gas commonly found in urban and suburban environments, mostly derived from anthropogenic emissions. Numerous epidemiological studies and several reviews from health and environmental agencies worldwide have reported that exposure to this pollutant is associated with adverse health outcomes, including increased short term mortality and morbidity. Evidence on the health impacts related to ozone exposure has important implications in climate change research, as ozone levels are predicted to increase with global warming.

2Short term ozone-mortality associations have been widely assessed in several multi-location time series studies in Europe, the United States, Canada, Latin America, and Asia. The general methodological framework consists of pooling location specific estimated risks, accounting for potential heterogeneity in the magnitude of the effect and uncertainty. In addition, the increased statistical power of multi-location analyses allows for the exploration of potentially complex features of the association (ie, non-linearity, delayed effects and harvesting, or differential risks by season). However, previous multi-location studies included a small number of cities and countries, were generally of limited geographical scope, and applied heterogeneous analytical approaches and modelling choices, making it difficult to draw consistent and comprehensive conclusions across different regions of the world.

3Although ozone-mortality associations have been widely assessed, results are rarely reported in terms of health impacts, such as excess deaths. Available figures are mostly derived from long term exposure metrics and risks estimated in specific subgroups, which are usually extrapolated to the general population. Quantification of health burdens from air pollution can be extremely useful for the design of efficient public health interventions, including the definition, assessment, and review of air quality standards. Current air quality standards vary greatly between countries, and only a few of them meet the stricter World Health Organization recommendation. Comparting the effects on health of ozone levels above different air quality standards can provide valuable insights into potential public health benefits achieved by strengthening current clean air policies. Although a few studies have attempted to tackle this problem, a widespread evaluation across several countries, which would help to identify more affected areas with a greater need for intervention, is still lacking.

4We carried out a multi-location time series analysis of mortality associated with short term exposure to ozone using data from 406 cities in 20 countries from multiple geographical regions. Next, we explored potential complexities of the association—namely, non-linearity, mortality displacement, and seasonality. Finally, we quantified the impacts on ozone associated mortality of specific concentration ranges consistent with the current air quality standards levels and then compared these estimates across countries.

5We initially extracted data for 434 locations across the 20 countries from the database of the Multi-city Multi-country (MCC) Collaborative Research Network (http://mccstudy.lshtm.ac.uk/) available at the time of the study. These include location specific daily mortality counts and environmental measures (weather and air pollutants) in largely overlapping periods from 1 January 1985 to 31 December 2015. For each location we derived daily time series of ozone (maximum eight hour average), particulate matter with an aerodynamic diameter less than or equal to 10 µm (PM10, per µg/m3, 24 hour average), particulate matter with an aerodynamic diameter less than or equal to 2.5 µm (PM2.5, per µg/m3, 24 hour average), nitrogen dioxide (24 hour average), total mortality, mean temperature (°C), and relative humidity (%). Mortality was represented by all cause deaths in Canada, the Czech Republic, Estonia, France, Germany, Greece, Italy, Japan, Mexico, Portugal, South Africa, South Korea, Sweden, Taiwan, UK, and US, whereas deaths due to non-external causes (eg, excluding self-intentional harm, poisoning) were used in Australia, China, and Spain, and non-external causes other than unintentional injuries in Switzerland (see supplementary eMethods 1 for the specific international classification of diseases codes used in each country). City specific air pollution series were derived from daily measurements of one or more monitors of the national or regional network. When more than one monitor was available, we computed the daily level of each pollutant (24 hour average or eight hour maximum) as the average across monitors of the city, consistent with previous multi-city studies. We excluded 28 cities as a result of poor quality data or limited periods (less than three years), with 406 locations included in the final analysis (see supplementary eMethods 1 for a detailed description of the data, exposure assessment, and exclusion criteria).

6The general statistical framework applied here is an extension of the classic two stage design and incorporates complex multivariable associations, hierarchical pooling methods, and the computation of impact measures. Briefly, we first estimated city specific ozone-mortality risks from separate time series regression models and then pooled these through a meta-analysis in the second stage. In a final step, we derived impact estimates, expressed as excess mortality fractions associated with ozone, from the pooled country specific risks and city specific exposure series. Using this general statistical framework, we performed a set of additional and sensitivity analyses to investigate specific features of the association. The analyses were conducted with R software (version 3.5.2) using the dlnm and mixmeta packages.

7In the first stage, we performed city specific time series analyses using generalised linear models with quasi-Poisson family. In this type of regression model, to properly scale the standard deviation of the coefficients proportionally to the potential overdispersion, a quasi-likelihood is applied. This phenomenon is common in these types of data, when the variability is larger than that expected under the assumption of a Poisson distribution. We assessed short term ozone-mortality associations using unconstrained distributed lag linear models. These models account for delayed effects of time varying exposures and quantify net effects over a predefined lag period. For the main model, we selected lag 0-1, estimating cumulative associations with the same and previous day’s exposures. The regression model included a natural spline of time with seven degrees of freedom each year, selected based on a quasi-likelihood version of the Akaike information criterion for 4, 6, 7, 8, and 10 degrees of freedom, and indicator variables for the day of the week, to control for long term, seasonal, and weekly variations in risk. Unlike in most previous studies on ozone, we applied a stricter control for temperature by using distributed lag non-linear models, an extension of distributed lag linear models for modelling complex non-linear and lagged associations. Following modelling choices applied in published analyses, we modelled the net temperature-mortality association over lag 0-21 (see supplementary eMethods 2).

8In the second stage we pooled city specific estimates through a multilevel meta-analysis. This novel meta-analytical model defines more complex random effects that can account for variations in risk across two nested grouping levels, represented by cities within countries. This approach allowed the derivation of improved estimates of ozone-mortality associations at both city and country level, defined as best linear unbiased predictions. Best linear unbiased predictions borrow information across units within the same hierarchical level and can provide more accurate estimates, especially in locations with small daily mortality counts or short series. We tested the presence of heterogeneity and reported it using multilevel extensions of Cochran Q test and I2 statistic. Association estimates, expressed as relative risk of mortality per 10 µg/m3 increase of ozone and 95% confidence interval were derived for each country from the corresponding best linear unbiased predictions.

9Risk estimates for ozone related mortality were then translated into impact measures, represented by excess mortality, following a method described elsewhere. Briefly, for each city we computed the daily number of deaths attributable to ozone (or daily excess deaths) using the corresponding risk estimate associated with the level of ozone in each day. Regarding the latter, we used country specific best linear unbiased predictions instead of the city specific estimates to avoid imbalances due to selection of cities and periods within each country. City specific estimates were reported as annual average number of excess deaths and 95% confidence intervals, so allowing for a proper comparison between locations with different lengths of study period. Then, country specific impacts were represented by excess mortality fractions (%) computed as the sum of the city specific daily excess deaths divided by the total mortality for each country. We used fractions instead of number of excess deaths, as excess deaths are not comparable across countries given the dependency on the denominator (ie, total mortality), which at the same time depends on the number of locations included. Although no evidence of a “safe” threshold exists, we computed associated deaths only for days with ozone levels above 70 µg/m3, as in previous health impact assessments. We considered this counterfactual scenario of 70 µg/m3 because ozone levels below this threshold could be mostly attributed to non-anthropogenic sources. A counterfactual scenario defined at 0 µg/m3 would not be appropriate either as it is not realistic given the ubiquitous presence of low levels of ozone derived from natural sources. We also disaggregated mortality impacts into contributions for exposure ranges above and between current air quality standards: 100 µg/m3 (WHO), 120 µg/m3 (European Union directive), 140 µg/m3 (National Ambient Air Quality Standard (NAAQS) in the US, about 0.070 parts per million), and 160 µg/m3 (Chinese Ambient Air Quality Standard (CAAQS) level 2).

10We performed a series of additional subanalyses to explore more complex features of the association, such as potential non-linearity, lagged effects, and seasonal differences. Firstly, we modelled exposure-response functions with a non-linear function consisting of a cubic B spline with internal knots at 50 µg/m3 and 60 µg/m3 of ozone. Secondly, we assessed delayed risks and potential mortality displacement by extending the lag dimension of the distributed lag linear model up to 30 days. Lag-response associations were modelled using a natural cubic spline with three internal knots placed at equally spaced lag values in the log scale. Thirdly, we assessed seasonal differences through interaction models between an indicator of season and the distributed lag linear model of ozone, as described elsewhere. We derived the ozone-mortality risk for the warm season (June to August in the northern hemisphere, December, January, and February in the southern hemisphere) and cold seasons (the remaining months).

11Modelling choices in the main model and extensions previously described were assessed and compared through the quasi-likelihood version of the Akaike information criterion and multivariate extensions of the Wald test. For sensitivity analyses, we first assessed changes in control for time trends and the potential confounding from other air pollutants (PM10, PM2.5, and nitrogen dioxide) and relative humidity by including each of these terms separately in the model. We then assessed the exclusion of a subset of US cities with data for summer only, which were included in the main analysis, and then different modelling approaches to control for temperature. See supplementary eMethods 1 and 2 for a description of the modelling details.

12This was a multinational collaboration using aggregated city level mortality and environmental data. Patients and members of the public did not contribute to the steering committee, design, or other areas of the study, which involved complex research methods and analysis.

13Table 1 provides a summary description of the data included for each country. A total of 45 165 171 deaths were analysed in the 406 cities, with an average time series of 13 years. Average annual mean ozone levels were widely heterogeneous across cities both between and within countries (fig 1). For example, lower levels were registered in Australian and northern European cities, whereas higher annual averages were found in some cities in the central area of the US, in Mexico, and in Taiwan. Supplementary eTable1 provides country specific descriptive summaries of the other air pollutants and humidity and eTable 2 reports the corresponding city specific descriptive results.

14On average, each 10 µg/m3 increase in ozone was associated with an overall relative risk of mortality of 1.0018 (95% confidence interval 1.0012 to 1.0024) (fig 2). Some heterogeneity was found across country and city specific risks (I2=29.8%, Cochran Q P<0.001). Larger risk estimates were found in the UK (1.0035 (1.0024 to 1.0046)), South Africa (1.0027 (1.0013 to 1.0042)), Estonia (1.0023 (1.0006 to 1.0040)), and Canada (1.0023 (1.0013 to 1.0032)), whereas Australia, China, the Czech Republic, France, Germany, Italy, Japan, South Korea, Sweden, Switzerland, and the US showed similar risks, ranging between 1.0014 and 1.0020. Lower and imprecise associations were estimated for Greece (1.0011 (0.9995 to 1.0028)), Mexico (1.0008 (1.000 to 1.0015)), Portugal (1.0011 (0.9997 to 1.0026)), Spain (1.0006 (0.9992 to 1.0019)), and Taiwan (1.0010 (0.9999 to 1.0021)). Supplementary eFigure 1 provides the corresponding figures with the relative risks for an increase in 10 parts per billion of ozone.

15Figure 3 depicts the excess mortality fractions above the WHO guideline and their distribution across intervals between the other air quality standards for each country, whereas supplementary eTable 3 and eTable 4 report the corresponding figures for excess fractions for total ozone (>70 µg/m3) and above and between air quality standards. Table 2 shows fractions and annual number of excess deaths associated with ozone for the total range of exposure and above the WHO guideline for a selection of the main cities in each country and overall across the 406 locations (supplementary eTable5 shows the estimates for all cities). Total mortality associated with ozone greater than 70 µg/m3 accounted for 0.26% of deaths (95% confidence interval 0.24% to 0.28%), which translates into 8203 annual excess deaths (95% confidence interval 3525 to 12 840) across the 406 locations studied (table 2). A substantial residual excess mortality of 0.20% (95% confidence interval 0.18% to 0.22%) corresponding to 6262 (95% confidence interval 1413 to 11 065) annual excess deaths remained when restricting to days with levels above the WHO guideline of 100 µg/m3. This proportion varied greatly by country, with considerably larger fractions in Mexico (0.52% (0.14% to 0.92%)) and Taiwan (0.37% (0.08% to 0.64%)) (fig 3, supplementary eTable 3). A mortality excess around 0.20% was estimated in Canada, China, Italy, Japan, South Africa, Switzerland, and the US, whereas France, Germany, South Korea, and the UK reported smaller percentages, ranging between 0.14% and 0.05% (fig 3, supplementary eTable 3). Imprecise or almost null estimates were found in the Czech Republic, Estonia, Greece, Portugal, Spain, and Sweden (supplementary eTable 3). Overall mortality fractions above more lenient air quality standards (ie, the European Union, NAAQS, and CAAQS) decreased progressively to 0.14%, 0.09%, and 0.05%, respectively (supplementary eTable 3). Only Mexico reported a considerably higher fraction, of 0.35% above the highest air quality standards of 160 µg/m3, although this finding was highly uncertain (black bar in fig 3, supplementary eTable 3). Null excess deaths were found in Australia, as daily exposure levels were all below 70 µg/m3. A similar pattern was found across estimates for the main cities in each country (table 2). A substantial number of annual excess deaths were associated with ozone levels above the WHO guideline—namely, 694 (95% confidence interval 22 to 1317) in the Valley of Mexico, 211 (112 to 307) in Los Angeles, 170 (40 to 304) in Tokyo, 128 (59 to 197) in Toronto, 82 (19 to 148) in Johannesburg, 48 (0 to 96) in Paris, and 37 (15 to 57) in London (table 2). Supplementary eTable 5 shows the corresponding estimates for the 406 cities.

16Additional analyses suggested no evidence of non-linearity in the concentration-response association (according to the quasi-likelihood version of the Akaike information criterion) (supplementary eFigure 2). The assessment of the lagged associations confirmed an immediate ozone-mortality association during the first week. However, lag specific estimates below 1 were found after the second week, which resulted in a slightly lower overall cumulative association of 1.0015 (95% confidence interval 0.9991 to 1.0032) when considering the delayed effects over the first 30 days after the exposure. Finally, no evidence of seasonal differences in ozone-mortality association were found (warm season: 1.0012 (95% confidence interval 1.000 to 1.0026); cold season: 1.0015% (1.0006 to 1.0024), Wald test P=0.37).

17Results from sensitivity analyses suggest that risk estimates of the main analysis were robust to the different modelling choices related to the control for time trends and adjustment by the three air pollutants and humidity (supplementary eTable 6). However, ozone-mortality risk estimates seemed to be sensitive to the approach to control for temperature (supplementary eFigure 3). We found larger ozone-mortality association estimates using less stringent control, although quasi-likelihood Akaike information criterion values suggested that the model with distributed lag non-linear model of temperature (main model) provided the best fit.

18On average, this study found an overall short term ozone-mortality association of 1.0018 (95% confidence interval 1.0012 to 1.0024) per 10 µg/m3 increase in ozone. This evidence is supported by previous epidemiological and experimental studies suggesting several pathophysiological mechanisms (e.g. systemic inflammation, haemostatic alterations). Larger associations were found in previous multi-country studies, including a subset of countries investigated here (eg, relative risk of 1.0022 in APHEA (Air Pollution and Health: A European Approach), 1.0026 in APHENA (Air Pollution and Health: A Combined European and North American Approach), per 10 µg/m3 increase), or single country studies (eg, relative risk of 1.0025 in the US (originally 1.0052 per 10 parts per billion increase), and China 0.55% per 10 µg/m3 increase, and 1.015 in Italy). Differences in the definition of the exposure variable (eg, moving average, single lag) and modelling approach could explain these discrepancies in the magnitude of the association. For example, compared with previous studies, we applied a stronger control for temperature (ie, distributed lag non-linear models), fully accounting for non-linearity and lagged temperature-mortality associations. In fact, results from sensitivity analyses are consistent with previous findings showing that ozone-mortality risk estimates were sensitive to the modelling strategy to control for temperature, reporting larger risks when using simpler approaches (supplementary eFigure 3). Moreover, one of the novelties of the applied statistical framework is the use of multilevel meta-analytical models in the second stage, properly accounting for heterogeneity across cities and countries.

19Our results showed important differences in the ozone-mortality association across countries. For example, while some areas such as UK, South Africa, Canada, and Estonia reported the largest risk estimates above 1.0020, smaller or imprecise estimates below 1.0011 were found in Greece, Mexico, Spain, and Taiwan. This unclear pattern would suggest that although several community level factors have been proposed as potential modifiers in single country studies (eg, population characteristics), these might not fully characterise differences between countries. Future multi-country studies are needed to provide further evidence on the factors defining the level of vulnerability of a population to air pollution.

20This study also provides evidence on the potential public health benefits of stricter clean air policies. In particular, we found that 0.20% excess mortality, which translates into more than 6000 deaths each year, related to short term exposure to ozone could have been avoided if ambient levels were below the WHO guideline of 100 µg/m3 in the 406 cities included in the study. Recent reviews found that most of the current air quality standards do not comply with the WHO air quality guideline, and that 80% of the world’s population in urban areas are exposed to air pollution levels above this threshold. Moreover, an additional 0.06% of excess deaths is associated with ozone levels between 70 µg/m3 and 100 µg/m3. These findings support the WHO initiative of encouraging countries to reconsider current air quality standards and enforce stronger emission restrictions and other public health interventions to meet its recommendations. Additionally, our results have important implications for healthcare practice. Apart from the implementation of clean air policies, individual strategies to reduce personal exposure to air pollutants are also desirable. In this regard, clinicians play an important role in counselling patients with potentially a higher susceptibility to adverse health outcomes related to air pollution. For example, professionals can advise sensitive individuals to stay indoors or avoid doing exercise during episodes of high ambient ozone.

21Previous studies showed that important health benefits could be achieved if reductions of ozone levels are reached. However, in this multi-country study we compared excess mortality estimates across air quality standards and countries, providing additional insights on specific areas with more urgent need of further interventions. For example, we found that 0.52% of total mortality in Mexico was associated with ozone above the WHO limit, the largest mortality fraction among the studied countries. This was associated with the highest ozone levels registered in the Mexican cities, especially above the 160 µg/m3 limit, which is close to its current air quality standards of 156 µg/m3. This means that attaining the current lenient standards would prevent a substantial proportion of ozone related deaths in this country. In contrast, results for the UK show a lower mortality fraction, despite the strongest ozone-mortality association, owing to the lower ozone levels registered in this country.

22This large epidemiological investigation on short term ozone-mortality associations included almost 50 million deaths from 406 cities in 20 countries from different regions across the world. Given its large sample size and wide geographical coverage, we were able to obtain consistent evidence of an association between short term exposure to ozone and total mortality. In addition, we provided ozone related impact estimates, quantified as excess mortality, across different air quality standards, countries and cities, providing evidence with important public health implications.

23We were able to explore additional complexities of the association by taking advantage of the large statistical power and advanced statistical techniques. Firstly, our results support the conclusions of previous studies on a generally linear concentration-response functions, with no indication of threshold. Secondly, we found evidence of a potential mortality displacement in the third and fourth week after the exposure. A similar lag pattern has been previously observed. However, potential mechanisms explaining this delayed and sustained pattern remain unclear. Finally, we found no evidence of seasonal differences in the ozone-mortality association. Previous multi-site studies have provided conflicting results, with larger risks in cold seasons in Asia and in warm seasons in the US and Europe. Further analyses are warranted to characterise different patterns across regions.

24This study has some limitations. Firstly, our results should not be considered truly global estimates, because several areas of the world such as South America, Africa, and the Middle East are unrepresented or were not assessed. In addition, the reported nationwide results might not be representative of the true impacts for some countries with a limited number of cities included in the study (eg, Sweden, Czech Republic, China). In particular, the estimated number of total excess deaths attributed to ozone should be interpreted as the sum of impacts in the 406 observed locations and not as total estimates across the 20 countries. Although excess fractions could be considered proper representations of the impacts for each country, the total excess number of deaths for each country is highly dependent on the total mortality considered in the study—that is, the number of locations included in each country. Systematic differences could also exist between countries in the characteristics of monitors (type, proximity to study area), study area boundaries, temporal coverage, data processing before data collection, and the collection of mortality data (eg, case ascertainment, codification). However, we ensured that the data fulfilled a minimum set of requirements for quality, a similar definition for the eight hour maximum metric, and location of the monitor (ie, within the study area or close enough to ensure its representativeness). Risks and impact estimates were only reported for total mortality (ie, deaths due to all or non-external causes) and we did not seek to identify the sources of heterogeneity of the results across countries. We acknowledge that the applied approach prevents us from understanding the potential mechanisms or differential susceptibility of the population, together with contextual differences across locations. Further studies are warranted to clarify this complex research question, including, for example, cause specific mortality and morbidity, and more complex two stage analyses. Finally, although the risk estimates were small they apply to the whole population, thus translating into substantial mortality impacts as shown in our estimates of excess mortality. By the same token, owing to the nature of the study design (time series analysis) the obtained excess mortality estimates refer to transient impact measures and not to the mortality burden or person years of life lost attributed to chronic exposure to ozone.

25This large multi-country study provided evidence on the short term association between ozone and mortality. We also show that clean air policies with the enactment of air quality standards can constitute essential public health tools to minimise the health burden. In particular, our results suggest that ozone related health impacts can be largely preventable by attaining effective air quality standards in line with the WHO guideline. Moreover, interventions to further reduce ozone pollution would provide additional health benefits, even in regions that meet current regulatory standards and guidelines. These findings have important implications for the design of future public health actions; particularly, for example, in relation to the implementation of mitigation strategies to reduce the impacts of climate change.

26Studies on the short term association between ground level ozone and mortality have been mostly performed in a few locations, in limited geographical areas, and using various designs and modelling approaches

27Although most of the studies found positive associations, results are heterogeneous, and a critical comparison across different countries and regions is made difficult by the limited statistical power and differences across studies

28Estimates of the association are usually reported as relative risks, a summary measure that does not quantify the actual health impact and makes it difficult to evaluate comparative health benefits of different regulatory limits

29This large multi-country study found increased mortality risks associated with exposure to ozone across locations and countries, with an average 0.18% per 10 µg/m3, reinforcing the evidence of a potential causal association

30Risk estimates were translated in measures of excess mortality, and it was found that more than 6000 deaths each year, corresponding to 0.20% of the total mortality, would have been avoided in the 406 cities studied if countries had implemented stricter air quality standards compliant with the WHO guideline

31Moreover, smaller but still substantial mortality impacts were found below WHO guideline, supporting the WHO initiative of encouraging countries to revisit current air quality guidelines and enforcing stronger emission restrictions to meet these recommendations