Effect modification by maximum temperature of the association between PM2.5 and short-term cardiorespiratory mortality and emergency room visits in Lima, Peru, 2010-2016.

BACKGROUND: The health effects of fine particulate matter (PM2.5) may be worse at higher temperatures.
OBJECTIVE: To investigate temperature's effect on PM2.5-mortality/morbidity associations in Lima, Peru.
METHODS: Time-series regressions relating PM2.5 and temperature to mortality and emergency room (ER) visits during 2010-2016. Daily PM2.5 levels (assigned to 40 Lima districts) and daily maximum temperature (Lima-wide) were estimated based on ground monitors, remote sensing, and modeling. We analyzed all-cause, cardiovascular (ICD codes I00-I99), and respiratory (ICD codes J00-J99) mortality, and cardiovascular and respiratory causes for ER visits.
RESULTS: The average PM2.5 concentration was 20.9 µg/m3 (IQR 17.5-23.5). The mean daily maximum temperature was 23.8 °C (IQR 20.8-26.9). PM2.5's effect on all-cause, respiratory, and circulatory disease mortality was significantly (p < 0.05) stronger at temperatures above the maximum temperature median. The rate ratios per increase of 10 µg/m3 of PM2.5 for all cause, respiratory, and circulatory mortality respectively were 1.03 (1.00-1.06), 1.04 (0.98-1.10), and 1.04 (0.98-1.10) at temperatures below the median, vs. 1.08 (1.04-1.12), 1.11 (1.03-1.19), and 1.14 (1.05-1.25) when temperatures were above the median. Results were analogous for ER visits for respiratory but not circulatory disease. SIGNIFICANCE: Results strengthen the evidence that air pollution may be more dangerous when temperatures are higher. IMPACT: Our data contribute to a growing body of literature which indicates that the damaging effects of PM2.5 may be worse at higher temperature, adding new evidence from Lima, Peru.

INTRODUCTION

There is evidence that temperature modifies the short-term effect of particulate matter on human health. In a recent systematic review and meta-analysis, Chen et al. [1] explored this issue with respect to mortality in a meta-analysis, reporting stronger effects of PM10 with warmer temperatures, particularly for respiratory causes. The meta-analysis results were based on 16 studies, none of which were from Latin America; most were from developed countries or China, which are likely to differ from Latin America in pollutant sources and levels, built environment, population age structure, and/or underlying health status (amongst other factors){2,3}.

Several relevant studies have been published since the review by Chen et al. [1], again primarily from high-income regions [4-5]. For example, Chen et al. [4] found that high air temperature enhanced the mortality effects of PM2.5 in 8 European cities. Liu et al. [5] studied the mortality effects of PM2.5 in 642 cities, and found that associations were stronger in locations with higher annual mean temperatures. A more recent review by Anenberg et al. [6] covered both mortality and morbidity studies and included 27 studies of particulate matter air pollution and heat. These authors concluded that “consistent exposure-response across studies led us to conclude that there is overall moderate quality and sufficient evidence for synergistic effects of heat and air pollution”.

We have previously studied the effect of PM2.5 on mortality and morbidity in Lima, Peru, over the period 2010-2016 [7-8]. Here we extend that work to consider the possible modification of the PM2.5 effects by maximum daily temperature. Understanding how temperature modifies the relationship between PM2.5 and adverse health outcomes can improve the design of interventions to protect vulnerable populations, including in the context of climate change when average temperatures will increase.

METHODS

Study area

Lima is located on the central coast of Peru at an average of 150 meters above sea level and covering an area of 2819 km2. The city has a population of ~10.7 million, representing roughly 30% of the national population. The city is comprised of 44 districts (including the bordering Callao district) (Figure 1). The contiguous province of Callao that adjoins Lima was included in the study of hospital admissions, but we did not have mortality data for Callao (about 10% of Lima’s population), which was hence not included in mortality analyses. PM2.5 data was lacking for 4 districts at high altitude (see below), and hence they were omitted for both mortality and hospital admissions analyses, but they represented only about 4% of Lima’s population.

Fig. 1

Map of Lima showing the city’s 44 districts (including bordering Callao), and location of the 9 hospitals (points) that contributed data to the study

Mortality data

Data on daily mortality in Lima from 2010-2016 were obtained from the country’s Ministry of Health. The dataset included information on the age, gender, district of residence, district of death occurrence, and cause of death coded by the International Classification of Disease 10th revision (ICD-10). Here we consider all-cause mortality, as well as mortality from respiratory (J00-J99) and circulatory (I00-I99) disease. We aggregated individual-level data by residential district (see Figure 1 for district boundaries) for the analysis (daily counts of deaths by district). These data have been used in prior epidemiological analyses [7].

Hospital data

We obtained electronic patient-level data from the emergency room (ER) visits of nine large public hospitals in Lima for the same time period (2010 – 2016), for the same ICD codes as the mortality data (J00-J99 and I00-I99). These hospitals belong to the Ministry of Health (see Figure 1). Ministry of Health public hospitals cover about 60% of the Peruvian population (https://www.who.int/workforcealliance/countries/per/en/). Public hospitals receive patients of low resources who do not have private health insurance, which in general is most of the population. In Lima there are 12 large full-service public hospitals, of which we studied 9 (three lacked electronic data). There are also private hospitals in Lima, but these are all relatively small. ERs in the larger hospitals receive outpatients and can also admit them when necessary .

For each visit, variables included the patient’s primary ICD-10 diagnosis code, district of residence, age, and gender. Visits were aggregated by day and district (based on patients’ district of residence) to obtain daily district-level visit counts. Visits by patients not residing in Lima were excluded.

To evaluate the accuracy of the electronic data, a team of two physicians compared the digital data received from each hospital with the hard copy medical history for a random sample of 100 records at each hospital. Overall, the two data sources matched 86-94% of the time for the key variables, and the electronic data (which we analyze here) was generally found to be more accurate. This process has been described in more detail in a previous publication [8].

PM2.5 and temperature data

Ground-monitoring PM2.5 data in Lima were available for March 2010 through December 2016, from ten stations from the Servicio Nacional de Meteorología e Hidrología del Perú (SENAMHI, Ministry of the Environment), and 6 stations operated during 2011-2012 by Johns Hopkins University (JHU) (Supplemental Figure 1). The SENAMHI data were not consistently available on a daily basis during our study period, covering only about 10% of days. Hence, the ground-monitoring network was considered too sparse to adequately capture the spatiotemporal variability in PM2.5 levels that occurs in Lima. Thus, we based our PM2.5 exposure data from a model developed by Vu et al. [9]; further details can be found there. Briefly, daily PM2.5 concentrations at a 1 km2 spatial resolution for 2010-2016 were estimated using a combination of the available ground measurements, aerosol optical depth (AOD) from satellites, meteorological and land use data, and outputs from a chemical transport model. AOD was obtained from NASA, using the MAIAC (Multi-Angle Implementation of Atmospheric Correction) algorithm. Meteorological fields (temperature, wind, barometric pressure, etc.) were obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF) and the Weather Research and Forecasting model coupled with Chemistry (WRF-Chem). A random forest model was used to regress the available ground measurements with 16 variables, including predicted MAIAC AOD from the gap-filling method, NDVI (vegetation index), percent urbanization, road category 3 distance, elevation, population density, relative humidity, surface solar radiation downwards, interpolated WRF-Chem simulated PM2.5, wind speed and direction, temperature, surface pressure, albedo (proportion of reflected light), cloud cover, and PBL (planetary boundary layer, the lowest part of the atmosphere). The overall cross-validation R2 value (and root mean square prediction error) was 0.70 (5.97 μg/m3), comparing predicted to observed ground level data. The mean difference between ground and predicted measurements was −0.09 μg/m3. This regression model was then used to predict population-weighted daily PM2.5 levels for each km2 grid across Lima. As mentioned above, due to concerns about model accuracy at high altitudes since all ground monitors were located below 150 meters, we excluded the four highest-altitude districts from the analysis; these represented only 4% of Lima’s population.

On every 16th day throughout the study period, we were unable to estimate PM2.5 due to lack of satellite coverage. Furthermore, PM2.5 estimates for October 15 to December 31, 2015 could not be made because the WRF-Chem model failed to estimate data within reasonable bounds for that period. Hence, we had PM2.5 estimates for 2236 (91%) of the 2465 day study period. Days missing PM2.5 were excluded from the analysis.

We estimated PM2.5 by district of residence, as district of residence was the level of detail available in mortality and morbidity data. We calculated a population-weighted average of 1 km2 grid estimates within each district to get the district average.

The Peruvian National Meteorology and Hydrology Service (SENAMHI) provided daily maximum temperatures (tmax) at 0.1° gridded spatial resolution (each grid approximately 10 km x 10 km) for the country of Peru. Temperature estimates were constructed by SENAMHI utilizing data from 178 quality-controlled, gap-filled, and homogenized air temperature stations, remote-sensed data, and a set of topographic predictors ( http://iridl.ldeo.columbia.edu/documentation/.pisco/.PISCOt_report.pdf). Estimates were constructed for daily maximum temperature for Lima by averaging grids that fall within the borders of Lima.

Statistical analysis, mortality

We used generalized linear models with quasi-Poisson regression (correcting for overdispersion) to estimate associations between daily district-level PM2.5 levels and daily counts of deaths for the outcomes of interest. PM2.5 effects were assessed using same day (lag 0), previous day (lag 1), two day (lag 2), and three days (lag 3), as well as the average of two (0-2) and three days (0-3), and finally the average for the prior month. Lag 1 for PM2.5 was eventually chosen based on superior fit based on Akaike’s Information Criteria (AIC). To control for spatially varying factors and allow the analysis to be based on temporal contrasts only, the models included indicator variables for district to represent the geographical area over which death counts were spatially aggregated [10]. The models also included variables for day of week, daily relative humidity (a linear term), and maximum daily temperature lagged 1 day (categorized into deciles). We compared two methods for controlling for long-term trends in mortality, either via parametric cubic splines with monthly knots, or with variables for month, year and an interaction between these variables (month*year). The month/year combination fit the data better as per the AIC, and was used. We also considered two methods of modeling temperature, either via deciles or a restricted cubic spline with 3 knots (25, 50, 75th percentile). Deciles fit the data better for the ER analyses, while the two models were equivalent for the mortality analyses. We used the decile model for both outcomes. The continuous PM2.5 variable was considered as a linear term, or log transformed, or considered as polynomials or cubic splines [11]. Linear models fit as well or better than other exposure metrics, again using the AIC to assess fit; therefore, we present the linear model results for mortality.

To assess the effect modification of temperature, we ran models separately for days when the maximum temperature (lagged 1 day) was either below or above the median maximum temperature across all study days. In both strata, low and high temperature, temperature itself was further controlled in the model via inclusion of temperature quintiles. Interaction between maximum temperature and PM2.5 was formally assessed with an interaction term between a dichotomous variable for above or below median maximum temperature and a linear continuous terms for PM2.5. To determine whether effect modification by temperature affected specific age groups differently, in supplemental analyses we stratified the data by two age groups (19-64, 65+; there were few events in those under 19). We also stratified the data by socioeconomic (SES) group using district-level census data, above and below the median poverty index as determined by the Peruvian census, to see whether effect modification by temperature differed by SES.

Standard errors of coefficients were adjusted for over-dispersion, which generally was very modest. Analyses were conducted using SAS v9.4 PROC GENMOD (SAS Institute Inc., Cary, NC, USA).

Statistical analysis, ER visits

Analyses of ER visit data were analogous to the analyses described above for mortality: quasi-Poisson regression adjusted for over-dispersion and including the same covariates. While we found here that splines controlling for long term time trends performed slightly better than the month/year/month*year interaction terms that we used in the mortality data, we also found that, when divided into high and low temperature, the models using splines for long-term trends sometimes did not converge; hence we again used month, year, and month-by-year interactions. We also found that same day PM2.5 measures out-performed all other lags, as judged by the AIC. For the shape of the exposure response function, we found that neither linear nor log transformed PM2.5 nor polynomial models performed as well as spline models, so we used restricted cubic splines with 3 knots at the 25th, 50th, and 75th percentiles [11](details of the spline curves are provided in supplemental material). We also used same day maximum temperature in these models to correspond to the same day PM2.5.

Results

The population-weighted average PM2.5 concentrations estimated for Lima, across all districts and years, was 20.9 ug/m3 (s.d. 4.91). The highest average concentrations (29.3 ± 4.5 ug/m3) were observed in East Lima and the lowest in Central Lima (17.7 ± 2.3 ug/m3) (Figure 1). Average PM2.5 levels in North and South Lima, were 22.7 (± 4.2) ug/m3, 20.5 (± 3.5) ug/m3 and 18.6 (± 1.6) ug/m3, respectively. All district averages exceeded the WHO long-term annual air quality standard of 10 ug/m3.

The average maximum temperature across Lima during the study period was 23.8°C, with a standard deviation of 3.7. The 25th, 50th, and 75th percentiles were 20.8, 23.6, and 26.9, with a minimum of 14.4 and a maximum of 34.6.

There were 183,893, 53,304, and 37,522 deaths from all causes, respiratory disease, and circulatory disease during the study. There were 70,453 and 582,443 ER visits for circulatory and respiratory disease, respectively, during the study period. The daily mean, standard deviation, maximum, and minimum for all deaths were 72, 11, 27 and 155 respectively. The same statistics were 21, 5, 6, 41 for respiratory deaths, and 15, 4, 1, and 31 for circulatory deaths. Of the respiratory visits,76% were among children, 18% were among those aged 19-64, and 6% among those aged 65 and older. Among the children, 58% of visits were for upper respiratory infections and another 23% were for asthma. Among circulatory disease, 4% were among children, 52% among adults aged 19-64, and 44% among adults aged 65+

There was little correlation between daily PM2.5 data and daily maximum temperature (Spearman rho =0.05). Spearman correlation between daily maximum temperature and monthly average maximum temperature was high (rho=0.93)

Figure 2 reports the rate ratios from the linear model for each mortality cause, for a 10 ug/m3 increase in PM2.5. For all three causes, risks from PM2.5 were significantly higher when daily temperatures were above the median maximum temperature compared to below it, and were significant (at the 5% level) only on those warmer days. The effect modification by temperature was most pronounced for circulatory mortality (RR 1.14 (1.05-1.25) for high temperature vs 1.04 (0.98-1.14) for low temperature) followed by respiratory (RR 1.11 (1.03-1.19) for high, 1.04 (0.98-1.10) for low) and then all-cause (RR 1.08 (1.04-1.12) for high, 1.03 (1.00-1.06) for low). Formal tests for interaction between PM2.5 and temperature were significant in all cases (all-cause p=0.0002, circulatory causes p=0.002, respiratory causes p = 0.02). Rate ratios corresponding to Figure 2 are reported in Supplementary Table 1. In supplemental analyses of the mortality data (Supplemental Tables 2 and 3) we found that the effect modification by temperature was generally stronger in the age 19-64 group compared to the age 65+ group. Regarding SES, those with higher SES showed more effect modification by temperature than those with lower SES. However, in many of these analyses, confidence intervals were wide and overlapping.

Figure 2.

Rate ratios per 10 ug/m3 increase in exposure from linear models of the PM2.5-mortality associations

Figures 3 and 4 present analogous results for ER visits; the figure shows the spline models for respiratory and circulatory visits to the ER, respectively, in relation to PM2.5, for temperatures below and above the median maximum temperature. Figure 3 shows results consistent with mortality results for respiratory causes, with increasing PM2.5 associated with increased respiratory ER visits only for temperatures above the median maximum temperature. The spline curves for low and high temperature are clearly separated for most of the range of PM2.5, although the confidence intervals cross at higher PM2.5 levels where the spline is less precise. On the other hand, Figure 4 shows that for circulatory ER visits, there is no association with PM2.5 for either low or high temperature (nor is there any effect of PM2.5 when the data are combined across temperature).

Figure 3.

Spline curve for the association of PM2.5 with respiratory hospital (ER) visits, Lima 2010-2016, for days below or above the median maximum Lima temperature.

Figure 4.

Spline curve for the association of PM2.5 with circulatory hospital (ER) visits, Lima 2010-2016, for days below or above the median maximum Lima temperature.

DISCUSSION

We have explored whether the short-term association of PM2.5 with adverse health outcomes is modified by ambient temperature in Lima, Peru. For mortality, our results suggest that risks are higher when temperatures are higher. This result also holds for emergency room visits for respiratory disease, but not for circulatory diseases.

In sub-group analyses (mortality only), there was some indication that the modification by temperature was stronger in younger populations and in those with higher SES. It is possible that younger people are more likely to be physically active and higher temperature exacerbates the PM2.5 effect more among those more active. Similarly, one can speculate that those with lower SES are more accustomed to heat and less susceptible to an exacerbation of the PM2.5 effect by heat. But there is very little data in the literature on these sub-group differences.

To our knowledge this is the only study of this issue in Spanish-speaking Latin America, and one of few studies conducted in a developing country (outside of China), where the PM2.5 constituents, demographic characteristics (e.g. age structure) and health profile of the population may differ from those in higher income countries. The majority of Lima’s population is low-income and our results may be particularly relevant to other low income populations.

Our results are largely consistent with the existing literature on effect modification by temperature on the association of PM2.5 and mortality. Of the 29 studies included in the systematic review of mortality by Chen et al. [1], the majority reported significant interactions with temperature for the causes explored here, where the effects of particulate matter were stronger when temperature was higher. Anenberg et al. [6], in a more recent review which included both mortality and morbidity studies, also concluded that there was evidence that the effects of particulate matter were greater when temperatures were higher.

There is also a literature on how temperature modifies the PM-morbidity relationship. For example, a recent study from New England, USA reported higher effects of PM2.5 on warmer days for respiratory hospital admissions, but the opposite for cardiac admissions [12]. The latter result is consistent with a different study from New York [13], but contrasts with recent findings from Cape Town, South Africa [14] that reported higher cardiovascular admissions on hotter days (neither study explored respiratory causes). A recent systematic review on a mix of causes concluded that the effects of PM2.5 on cause-specific morbidity seem not to be significantly modified by season [15]. In our morbidity data the effect of temperature on the association of PM2.5 and ER visits was only observed for respiratory causes, in contrast to our mortality findings, where we found an effect for both respiratory and circulatory causes. It is possible that the mechanism by which temperature modifies PM2.5 on mortality may affect many causes, while the effect modification for morbidity is restricted primarily to respiratory disease. Perhaps temperature variability is relatively more important for cardiovascular morbidity versus mortality10, and PM2.5 might not affect variability as much as underlying cardiovascular pathology.

Aside from the modifying effect of temperature, it is surprising that we find an overall effect of PM2.5 on circulatory mortality but not circulatory ER visits. It is possible that ER visits for circulatory disease are a broad category, and that only a subset of these reflect disease which might be affected by PM2.5. For example we found that ischemic heart disease (IHD) (14% of ER visits) did show a relatively strong effect of PM2.5, when restricted to those 19-64 (RR 1.14 for a 10 ug/m3 increase). IHD is a more common subset in mortality, representing 25% of circulatory disease mortality in our data. In the IHD aged 19-64 subset of ER visits, we did find a strong modifying effect of temperature (RR per 10 ug/m3 increase in PM2.5, 1.25 (0.97-1.61) when temperature>=median, and 1.01 (0.80-1.26 when temperature < median). Our results should be interpreted in light of the study’s limitations. Due to poor spatial and temporal coverage of the ground monitoring system, we relied on a remote sensing-based PM2.5 exposure model, which required us to exclude a small number (n=4) of high-altitude districts and days when satellite AOD was not available. However, the excluded districts represented only 4% of Lima’s population. In addition, our morbidity data did not include every hospital in the city, and so may be a non-random selection of Lima’s population. However, we believe our hospital data represent about half the population of Lima, and that the relationship between air pollution and ER visits for this part of the population is unlikely to differ substantially from the rest of Lima, although we have no data to confirm this. A further limitation is the lack of data on occupation and the location of occupation, which could affect PM2.5 exposure, which we based on location of residence.

Both Annenberg et al. [6] and Chen et al. [1] discuss possible mechanisms by which higher temperature may increase the damage of air pollution on the lung and heart. Heat increases systemic and pulmonary inflammation as part of thermoregulation. Heat also increases cardiovascular problems via inflammation and blood clotting. In turn, air pollution increases oxidative stress leading to pulmonary inflammation, and also may affect cardiovascular disease by inflammation and altered coagulation. These shared mechanisms may act synergistically.

This study contributes to the literature on how temperature modifies the short-term effect of fine particulate matter on both mortality and morbidity, adding to the sparse evidence from developing countries. Understanding this relationship, including by cause, is important in the design of air pollution warning systems in order to identify dangerous air pollution days and to protect the most vulnerable populations. The results also suggest that, all else equal, air pollution effects may become more problematic in the future given the increasing temperatures expected with climate change.