Time-series analysis of ambient PM2.5 and cardiorespiratory emergency room visits in Lima, Peru during 2010-2016.

INTRODUCTION: There have been no time-series studies of air pollution in Peru. Here we evaluate the effect of ambient PM2.5 on emergency room (ER) visits in Lima.
METHODS: We estimated daily PM2.5 levels at a 1 km2 resolution during 2010-2016 using ground measurements, satellite data, and chemical transport model simulations. Population-weighted average daily PM2.5 levels were calculated for each district in Lima (n = 40), and assigned to patients based on residence. ER visits for respiratory and circulatory diseases were gathered from nine large public hospitals. Poisson regression was used to estimate the rate ratio for daily ER visits with change in daily PM2.5, controlling for meteorology, time trends, and district.
RESULTS: For each interquartile range (IQR) increase in PM2.5, respiratory disease ER visits increased 4% (95% CI: 0-5%), stroke visits 10% (3-18%), and ischemic heart disease visits (adults, 18-64 years) 11% (-1, 24%). Districts with higher poverty showed significantly stronger associations of PM2.5 and respiratory disease ER visits than districts with lower poverty. Effects were diminished 24-42% using Lima-wide instead of district-specific PM2.5 levels.
CONCLUSIONS: Short-term exposure to ambient PM2.5 is associated with increases in ER visits in Lima for respiratory diseases and stroke, and among middle-aged adults, ischemic heart disease.

INTRODUCTION

Air pollution is considered a major problem for environmental health, being one of the main causes of morbidity and mortality around the world (Achilleos et al., 2017). Of all air pollutants, fine particulate matter (PM2.5) has been most consistently associated with health effects (Pope et al., 2015). Many previous epidemiological reviews have shown an association between adverse pulmonary and cardiovascular diseases and PM2.5 exposure in both adults and children (Apte et al., 2015; Kim et al., 2015). The World Health Organization (WHO) estimated that air pollution is associated with premature deaths related to ischemic heart disease, strokes, chronic obstructive pulmonary disease, acute lower respiratory infections and lung cancer (WHO, 2018). Time-series studies are commonly used to estimate short-term associations between daily concentrations of air pollution and daily reports of health outcomes (Bashkaran et al., 2013). These types of studies have also identified ambient PM2.5 as a risk factor of cardiopulmonary disease morbidity and mortality (Requia et al, 2018; Achilleos et al, 2017; Pope et al, 2015).

A WHO report regarding global outdoor air pollution in 2014 showed that Lima, the capital of Peru, was one of the top three cities in the Americas with highest mean annual PM2.5 values (WHO, 2014). Average PM2.5 levels during 2010–2015 were 26 μg/m3 (Silva et al, 2017). Using air pollution data for Lima during 2001–2011, Gonzales and Steenland (2014) estimated that air pollution was responsible for 2,300 premature deaths related to cardiorespiratory disease in adults per year. Although, these data suggested an important public health problem in Lima and its surroundings, there have been few direct studies on the association between PM2.5 and health risk in Peru (Underhill et al, 2015; Robinson et al, 2011, Carbajal-Arroyo et al. 2007).

Here we conduct a daily time-series analysis to determine the association between daily district-level ambient PM2.5 concentrations and ER visits for cardiopulmonary diseases in Lima, Peru. We used a recently developed model of daily PM2.5 levels in Lima, with specific estimates for its 40 districts during the period 2010–2016 to assign exposures (Vu et al., 2019). Air pollution is heterogeneous in Lima, due to different population densities and wind patterns. Daily district-specific levels, assigned to patients living in these districts, were used to more accurately assign air pollution levels with the goal of obtaining more accurate and precise PM2.5 health effect estimates compared to assigning Lima-wide daily PM2.5 levels.

MATERIALS AND METHODS

Study area

Lima is located on the central coast of Peru, at an average of 150 meters above sea level. Lima is comprised of 44 districts (including the bordering Callao district) and divided into four zones: North Lima, Central Lima, East Lima, and South Lima. The contiguous province of Callao that adjoins Lima was included in the study, and considered as a fifth zone, the West zone (Figure 1). We excluded 4 districts at high altitude (all above 570 meters average altitude) on the eastern edge of the study area due to uncertainty about the PM2.5 model predictions in these districts. The uncertainty was largely driven by the fact that ground monitoring stations providing inputs to the PM2.5 model were all located below 375 meters, requiring a large extrapolation to these four high districts, using the model’s prediction of the altitude effect. These districts (Carabayllo, Chaclacayo, Cienguilla, Lurigancho) represented only 4% of the total population of Lima and Callao.

Figure 1.

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

Hospital data collection

The study was performed using electronic patient visit-level data obtained from the emergency rooms of 9 large public hospitals in Lima for the period March 1, 2010 to December 30th, 2016 (Figure 1). These hospitals belong to the Ministry of Health (MoH). They are located in the four zones of Lima: North (Carlos Lanfranco and Cayetano Heredia hospitals), Center (Instituto Nacional de Salud del Niño, A. Loayza, and Santa Rosa hospitals), South (Maria Auxiliadora hospital), East (Vitarte and San Juan de Lurigancho hospitals) and Callao (Daniel A. Carrión hospital).

For each visit, variables included the patient’s primary International Classification of Diseases 10th Revision (ICD10) 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 for selected disease categories of interest. Respiratory diseases (RD) were analyzed as a group and included visits with ICD10 codes J00 - J45; two RD subcategories of infectious RD (ICD10 codes J00-J06, J09-J22) and non-infectious RD (ICD10 codes J30-J45) were also analyzed. We also analyzed ischemic heart disease (I20 - I25), and stroke (G45, I63 - I67).

To evaluate the electronic patient visit data, we compared the digital data received from each hospital with the hard copy medical history, for a random sample of 100 records at each hospital. Our validation team included two physicians. Among the variables collected for this study, date of emergency care had the highest matching rate (94%). In cases with a mismatch in dates, the disagreement typically was one or two days off, with the medical history date being earlier than the date in the digital record; this discrepancy, extrapolated to our full dataset, may have led to some mismeasurement of exposure in 6% of our data. ICD10 diagnosis code had the lowest matching rate at 86%. The disagreement observed was generally for cases where the patient presented diffuse symptoms in which there was no precise diagnosis, then he or she was hospitalized for a few days. This lack of a definitive diagnosis was recorded in the hard copy emergency room record. However upon further searching the hospital database we found that the final diagnosis upon hospital discharge matched the electronic visit data that we received. Hence, we believe that in most cases of diagnosis discrepancy, the electronic data were correct. We found some mismatch among the electron and hard copy records for the variables age, sex, and district; for these variables, the disagreement was more due to missing data in the hard copy than any true discrepancy.

Ambient PM2.5 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)(Bose et al. 2018)(Figure 2). 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. (2019). Briefly, daily PM2.5 concentrations at a 1 km2 spatial resolution for 2010–2016 were estimated using a combination of the available ground measurements plus aerosol optical depth (AOD) data from satellites, and meteorological and land use data and chemical transport model-WRF-Chem (Grell et al., 2005). AOD was obtained from NASA, using the MAIAC (Multi-Angle Implementation of Atmospheric Correction) algorithm. Meteorological fields (temperature, wind, and barometric pressure) were obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF) and the Weather Research and Forecasting model coupled with Chemistry model (WRF-Chem) (Sánchez-Ccoyllo et al. (2018). A random forest model was used to regress the available ground measurements with 14 variables, including MAIAC AOD, meteorological variables from WRF-Chem and ECMWF, and land use variables. The overall cross-validation (CV) R2 value (and root mean square prediction error) was 0.70 (5.97 μg/m3), comparing predicted to observed ground level data. The observed mean PM2.5 for ground measurements was 24.7 μg/m3 while the prediction model estimated mean PM2.5 was 24.9 μg/m3 in the cross-validation dataset. The mean difference between ground and predicted measurements was −0.09 μg/m3. This regression model was then used to predict daily PM2.5 levels for the full 1 km2 grid across Lima. For use in epidemiologic analyses, for which daily ER visit counts were aggregated by district of patient residence, daily population-weighted average PM2.5 levels were calculated for each district from the 1 km2 gridded data

Figure 2.

Map of Lima showing the locations of ten automatic PM2.5 monitoring stations operated by the Ministry of the Environment (MINAM/SENAMHI) and six PM2.5 stations operated by a John’s Hopkins University (JHU) research group.

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, 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 days (91%) out of the 2465 days during the study period.

Statistical analysis

We used Poisson generalized linear models to estimate associations between daily district-level PM2.5 levels and counts of ER visits for the outcomes of interest. PM2.5 effects were assessed using same day (lag 0), previous day (lag 1), and day before previous (lag 2), as well as the average across all 3 days, in separate models. 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 ER visit counts were spatially aggregated; this also controlled for spatial autocorrelation in the baseline ER visits across the districts (Sarnat et al., 2013). The models included indicator variables for day of week. Long-term trends in ER visit rates were controlled with parametric cubic splines with monthly knots each year. Meteorology was controlled using terms for daily lag 0 mean temperature (linear, quadratic, cubic), and daily lag 0 mean relative humidity. Indicator variables for each hospital for each day, indicating whether or not the hospital contributed any cases for that day, were included to account for missing data for relatively short periods at some hospitals.

To assess effect modification by socioeconomic status (SES), we obtained estimates of the percent of each district’s households above the poverty line from the Peruvian Institute of Statistics and Computer Science (INEI), which collects census data. Poverty is defined by INEI as a household’s per capita spending insufficient to purchase a basic food and nonfood basket (clothes, housing, health, etc). We then divided districts above and below the median poverty percentage (which was 12.9%), added this dichotomous variable to our model, and tested an interaction term with PM2.5.

We also assessed effect modification by age. Using information on patient age, we separated our outcome visit counts into three age strata, young (<18), middle-aged (19–64), and old (65+), to obtain daily district-level outcome counts for each age strata. We then ran our model separately for all age-specific outcomes. Differences in the observed association of PM2.5 between age strata was considered indicative of effect modification by age.

Standard errors of coefficients were adjusted for over-dispersion, which generally was very modest. Best model fit was determined via the Akaike Information Criterion (AIC). Analyses were conducted using SAS v9.4 (SAS Institute Inc., Cary, NC, USA).

In developing our final epidemiologic model structure, we began by controlling for temporal trends. We considered cubic spline terms with 6 knots per year as a priori superior to other measures such as month and year (and their interaction), based on our prior experience. We then tested the splines vs other methods of controlling for temporal trends, and indeed the splines outperformed other measures as evidenced by a better fit to the data (lower AIC). We then added each other variable chosen a priori as a candidate to be in the model, and we explored functional forms separately for each, and chose the best form via lowest AIC. For example, for temperature we determined that the inclusion of linear, quadratic, and cubic terms significantly improved the model over inclusion of a linear term, or a linear and quadratic term. For relative humidity we found a simple linear term fit best. Other terms such as variables for district, hospital, day of week, and hospital were inherently categorical, and we included these again based on significant lowering of the AIC.

We also tested whether those days missing PM2.5 differed with regard to the number of daily health events for all our major outcomes. To do so, we created a dichotomous variable for days missing PM2.5, and ran our health outcomes models using this dichotomous variable in place of PM2.5. ‘Days missing PM2.5’ was not associated with any outcome, with all p-values being >0.80.

Computer code for analyses is available upon request to authors. This study was approved by the Institutional Review Board (IRB) of the Universidad Peruana Cayetano Heredia.

RESULTS

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

Table 1.

Average concentrations of PM2.5 (μg/m3) by zone and district in Lima during 2010–2016.

District	X	SD	District	X	SD
North Lima
Ancon	22.3	2.43	Puente Piedra	27.2	2.93
Comas	27.3	3.32	San Martin de Porras	18.3	2.19
Independencia	23.4	2.43	Santa Rosa	21.1	1.97
Los Olivos	19.3	2.3
Center Lima
Cercado de Lima	18.4	2.13	Miraflores	16.9	1.57
Barranco	16.8	1.46	Rimac	20.4	2.35
Breña	17.6	2.21	San Borja	19.6	2.27
Jesus María	16.5	2.3	San Isidro	17.0	1.97
La Victoria	19.2	2.25	San Luis	20.3	2.22
Lince	17.3	2.24	San Miguel	17.1	1.31
Magdalena	16.3	1.51	Santiago de Surco	20.3	1.79
Pueblo Libre	16.8	1.68	Surquillo	17.4	1.93
South Lima
Chorrillos	17.9	1.25	San Bartolo	21.3	2.79
Lurin	18.6	1.39	San Juan de Miraflores	20.3	1.87
Pachacamac	27.7	1.55	Santa Maria	18.0	1.11
Pucusana	17.9	1.03	Villa El Salvador	19.4	1.95
Punta Hermosa	19.6	1.82	Villa Maria del Triunfo	24.7	2.14
Punta Negra	19.8	2.1
East Lima
Ate	29.0	4.12	San Juan de Lurigancho	32.1	4.91
El Agustino	27.4	3.77	Santa Anita	28.7	4.86
La Molina	29.1	3.38
West Lima
Callao	18.6	1.58

Table 2 provides descriptive statistics on the number of cases of respiratory and circulatory ER admissions by age., gender, and zone in Lima.

Table 2.

Summary of the respiratory and circulatory ER visits at 9 participating hospitals in Lima during 2010 to 2016.

Characteristic	Respiratory (n=595,174 visits total)		Circulatory (n=71,984 visits total)
	N	%	n	%
Patient Age (years):
0 – 18	450,488	75.7	3,418	4.8
19 – 64	109,066	18.3	37,465	52.1
>65	35,620	6.0	31,101	43.2
Patient Sex:
Female	291,635	49.0	39,994	55.6
Male	303,539	51.0	31,990	44.4
Patient Zone of Residence in Lima:
North	134,490	22.6	16,921	23.5
Center	173,430	29.1	19,910	27.7
South	96,490	16.2	15,332	21.3
East	128,038	21.5	9,039	12.6
West	62,726	10.5	10,782	15.0

Same day (lag 0) PM2.5 levels generally produced the best fit models between daily air pollution and ER visits, as judged by the AIC, compared to models estimating lag 1 and lag 2 effects. Lag 0 model results are presented in Table 3, for each disease category of interest (lag 1 and lag 2 model results are presented in Supplementary Tables 1–2.). Results are presented as rate ratios (RRs) and 95% confidence intervals (CIs) calculated for an interquartile range (IQR) increase in PM2.5 across all districts and years (6.1 μg/m3). PM2.5 was significantly positively associated with all respiratory disease ER visits, and the sub-categories of infectious and non-infectious respiratory disease, in age groups under 65. We found significant positive RRs for stroke ER visits for ages >18 years, and a borderline significant RR for ischemic heart disease ER visits for persons aged 18–64 years.

Table 3.

Associations of same day (lag 0) district-level PM2.5 and ER visits for respiratory and circulatory diseases. Effect estimates presented as rate ratios (RR) and 95% confidence intervals (CI) per interquartile range (IQR) increase in PM2.5.*

Disease	Age Group	n	RR	LCL	UCL	p-value
Respiratory	all	595,174	1.04	1.03	1.05	<0.0001
	<18	261,750	1.03	1.02	1.04	<0.0001
	18–64	109,666	1.09	1.06	1.11	<0.0001
	65+	35,620	1.02	0.98	1.06	0.3440
Infectious Respiratory	All	376,333	1.05	1.04	1.06	<0.0001
	<18	304,075	1.04	1.03	1.05	<0.0001
	18–64	57,055	1.10	1.07	1.13	<0.0001
	65+	15,203	1.03	0.98	1.09	0.1832
Non-infectious respiratory	All	218,841	1.03	1.02	1.05	<0.0001
	<18	146,413	1.03	1.01	1.04	0.002
	18–64	52,011	1.08	1.05	1.11	<0.0001
	65+	20,417	1.01	0.96	1.06	0.591
Stroke	All	10,239	1.10	1.03	1.17	0.0034
	<18	195	Did not converge
	18–64	4,262	1.11	1.01	1.02	0.03
	65+	5,872	1.10	1.01	1.20	0.02
Ischemic heart disease	all	5,134	1.02	0.96	1.15	0.27
	<18	83	Did not converge
	18–64	3,059	1.11	0.99	1.25	0.07
	65+	1,992	0.96	0.82	1.12	0.59

IQR for district level PM2.5 was 6.1 μg/m3. Respiratory diseases (RD) J00 – J45, infectious respiratory disease (codes J00–J06, J09 – J22), non-infectious respiratory disease (codes J30–J45). Ischemic heart disease (I20 – I25), and stroke (G45, I63 – I67). Models adjusted for district, temperature, relative humidity (RH), day of week (DOW) and hospitals.

We found effect modification by SES for respiratory disease ER visits, including both visits for infectious and non-infectious respiratory disease. Rate ratios per IQR for districts above and below the median percentage of population living below the poverty line are presented in Table 4. Effect modification was present for all respiratory disease and respiratory disease in those under 18. In all models, a significant effect of PM2.5 in increasing ER visits was found only for those living in poorer districts. SES did not modify associations for stroke or ischemic heart disease ER visits (results not shown).

Table 4.

Associations of same day (lag 0) district-level PM2.5 and ER visits for respiratory by age group, for outcomes where there was significant (<0.05) effect modification by socioeconomic status (SES). Effect estimates presented as rate ratios (RR) and 95% confidence intervals (CI) per interquartile range (IQR) increase in PM2.5.

Disease	Age Group	Poverty Level*	RR	LCL	UCL	p-value
Respiratory	all	richer	1.01	0.99	1.03	0.29
	all	poorer	1.06	1.02	1.10	<0.0001
	0–18	richer	0.99	0.97	1.01	0.37
	0–18	poorer	1.05	1.01	1.09	<0.0001
Infectious Respiratory	all	richer	1.02	1.00	1.05	0.04
	all	poorer	1.07	1.02	1.12	0.0001
	0–18	richer	1.00	0.98	1.02	0.87
	0–18	poorer	1.06	1.01	1.11	<0.0001
Non-infectious respiratory	all	richer	0.99	0.97	1.02	0.56
	all	poorer	1.05	1.00	1.11	<0.0001
	0–18	richer	0.98	0.95	1.01	0.20
	0–18	poorer	1.04	0.97	1.11	0.0003

Richer districts were those where the average household poverty level, as defined by the Peruvian census, was lower than the median percentage of households living in poverty across all districts, which was 12.4%. Poorer districts had poverty levels of 12.4% or greater.

As a sensitivity analysis, we estimated the same associations using daily PM2.5 levels for Lima as a whole. Depending on average time-activity-location and mobility patterns of the patient population, Lima-wide PM2.5 levels are less likely to reflect patients’ exposures to ambient PM2.5 compared to PM2.5 levels in the districts where they live. Thus, models using Lima-wide PM2.5 may be subject to more measurement error and be biased to the null compared to our primary analyses using district-specific PM2.5. Results showed that in general models using PM2.5 in Lima as a whole had lower exposure-response coefficients than did models using district-specific PM2.5. For example, the exposure-response coefficient for respiratory disease using Lima-wide PM2.5 was 24% lower than that using district-specific PM2.5; the coefficient for stroke was 42% lower.

DISCUSSION

Environmental epidemiological studies of daily air pollution have mostly been conducted in developed countries, which often have good daily data on air pollutants. In low and middle income countries (LMIC), the scarcity of monitoring stations and lack of daily measurements make time-series studies of air pollution difficult. Researchers have used a variety of techniques to develop predictive models that can estimate daily levels absent adequate ground monitoring data (Liu et al, 2009; Miri at al, 2019; Paciorek and Liu, 2012, Levy et al 2010). In our case, we used a satellite-driven PM2.5 exposure model developed by Vu et al. (2019) that incorporated ground-based measurements of PM2.5, satellite data, and chemical transport model simulations, and ultimately provided daily population-weighted average PM2.5 concentrations for all districts of Lima for the 2010 to 2016 study period. The present study is thus the first that analyzes the association between PM2.5 exposure and respiratory and circulatory disease in Lima using a time-series approach.

We found positive associations of ambient PM2.5 with respiratory disease ER visits, with an increase of 4% per IQR (6.1 μg/m3) increase in PM2.5 driven largely by those under 65. These results are reasonably congruent with other studies. For example, a recent systematic review of 16 time-series studies of hospital admissions found arespiratory disease risk of 2.7% per 10 μg/m3 increase in PM2.5. (Requia et al. 2018).

We found a significant association between PM2.5 and ER visits for stroke in adults over 18 years, and a borderline significant association with ER visits for ischemic heart disease in adults ages 18–64 years. The literature on ER visits for stroke in relation to PM2.5 is relatively sparse; Chen et al. (2014) found a positive association similar to ours, while another study in China did not find an effect (Ge et al. 2018). A study in Arkansas found an effect only in winter (Rodopoulou et al. 2015). For ischemic heart disase, literature supports such an association between either ischemic heart disease or heart attacks identified using ER visit records (Zhang et al. 2016, Xu et al. 2016, Weber et al. 2016, Weichenthal et al. 2016).

We found that a model with no lag fit better than a model with a lag of 1 day or two days, or an average of 3 days. Previous research on ambient air pollution and cardiorespiratory ED visits has shown that relevant lags range from lag 0 up to 7+ days, sometimes depending on the outcome and pollutant (e.g., Metzger et al 2004, Peel et al. 2005, Strickland et al. 2010).

The association between PM2.5 concentrations and ER visits for respiratory and circulatory diseases in our study was evident at concentrations which are in the range of the annual permissible level in Lima (25 μg/m3), and considerably below the permitted 24 hour level (50 μg/m3). (https://busquedas.elperuano.pe/normaslegales/aprueban-estandares-de-calidad-ambiental-eca-para-aire-y-e-decreto-supremo-n-003-2017-minam-1529835-1/). Average population-weighted PM2.5 was 20.9 μg/m3 for Lima, but there is considerable variation in PM2.5 levels across the city, with highest concentrations in the east and north of Lima, due to the pattern of local winds entering from the coast in a southwest direction. Some government regulations have been implemented to reduce pollution. One of the most successful regulations has been the vehicular restructuring in downtown Lima in 2011 (lowering speed limits, restricting bus stops), where a reduction of 60% in the PM2.5 average was observed (Tapia et al, 2018). Our results suggest that more such measures should be considered.

We also found significant modification of PM2.5 effects by SES, whereby districts with lower SES (poorer) had higher rate ratios due to PM2.5, for respiratory disease ER visits than districts with higher SES (richer); effect modification by SES was not observed for stroke or ischemic heart disease. It is possible that those who are more poor are more likely to use the ERs for respiratory disease for children than those who are richer. But, among adults faced with stroke or a heart attack, both rich and poor end up in the ER. Effect modification like this has been seen in most other studies that have examined the issue, but findings are not consistent. Clougherty et al. (2014), in a review, cited three studies of mortality and hospitalization, two of which found worse air pollution effects with lower SES, and the third found no modification. O’Lenick et al. (2017) found greater air pollution effects on asthma ER visits among children among those with lower SES. These same authors cite 17 prior population-based studies of childhood asthma in relation to air pollution where SES effect modification was studied. In seven studies, there was no modification by SES; in eight studies, lower SES was associated with worse pollution effects; and in two studies, higher SES was associated with worse air pollution effects. In newer studies not cited by O’Lenick et al., stronger air pollution health associations with lower SES was found by Cakmap et al. (2016) for children/respiratory disease, Wang et al. (2017) for mortality, Chi et al. (2016) for cardiovascular disease; in contrast Goodman et al. (2017) found no modification by SES for asthma in children.

The main strength of our study was the use of a predictive model that estimated daily PM2.5 concentrations for each district of Lima. This model allowed the construction of consistent long-term historical measurements to supplement the relatively sparse data from ground monitoring, and the assignment of PM2,5 levels to people based on their district of residence. We found stronger associations with ER visits using district-specific PM2.5 estimates than using Lima-wide estimates, suggesting that our primary measure of exposure had less exposure measurement error than the traditional use of central site or city-wide average pollutant data in air pollution time-series studies. Another strength was our use of ER visits of public hospitals of the MoH. In general MoH 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. There are also private hospitals in Lima, and there are three public hospitals that we did not study due to lack of data, but these are all relatively small. Of the public hospitals, we studied only 9 of the 12 due to lack of available electronic data from 3 hospitals. We believe our 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. Additional limitations of our study include low power to capture some specific sub-categories of disease, or to explore in more detail effect modification by other factors.

Conclusions

The findings from this first time-series study of ambient PM2.5 and morbidity in Lima indicate that short-term exposure to ambient PM2.5 is associated with increases in emergency room visits for respiratory diseases, stroke, and ischemic heart disease. These associations were observed at PM2.5 concentrations similar to the annual permitted level and considerably below the daily permitted level in Lima, which suggests that current PM2.5 standards may not be adequate for the health of the Lima population.