Long-Term PM10 Exposure and Cause-Specific Mortality in the Latium Region (Italy): A Difference-in-Differences Approach.

BACKGROUND: The link between particulate matter (PM) exposure and adverse health outcomes has been widely evaluated using large cohort studies. However, the possibility of residual confounding and lack of information about the health effects of PM in rural and suburban areas are unsolved issues.
OBJECTIVE: Our aim was to estimate the effect of annual PM≤10µg (PM10) exposure on cause-specific mortality in the Latium region (central Italy, of which Rome is the main city) during 2006-2012 using a difference-in-differences approach.
METHODS: We estimated daily PM10 concentrations for each 1 km2 of the region from 2006 to 2012 by use of satellite data, land-use predictors, and meteorological parameters. For each of the 378 regional municipalities and each year, we averaged daily PM10 values to obtain annual mean PM10 exposures. We applied a variant of the difference-in-differences approach to estimate the association between PM10 and cause-specific mortality by focusing on within-municipality fluctuations of mortality rates and annual PM exposures around municipality means, therefore controlling by design for confounding from all spatial and temporal potential confounders. Analyses were also stratified by population size of the municipalities to obtain effect estimates in rural and suburban areas of the region.
RESULTS: In the period 2006-2012, we observed deaths due to three causes: 347,699 nonaccidental; 92,787 cardiovascular; and 16,509 respiratory causes. The annual average (standard deviation, SD) PM10 concentration was 21.9 (±4.9) µg/km3 in Latium. For each 1-µg/m3 increase in annual PM10 we estimated increases of 0.8% (95% confidence intervals (CIs): 0.2%, 1.3%), 0.9% (0.0%, 1.8%), and 1.4% (-0.4%, 3.3%) in nonaccidental, cardiovascular, and respiratory mortality, respectively. Similar results were found when we excluded the metropolitan area of Rome from the analysis. Higher effects were estimated in the smaller municipalities, e.g., those with population < 5,000 inhabitants.
CONCLUSION: Our study suggests a significant association of annual PM10 exposure with nonaccidental and cardiorespiratory mortality in the Latium region, even outside Rome and in suburban and rural areas. https://doi.org/10.1289/EHP3759.

Introduction

The association between particulate matter (PM) with aerodynamic diameter ( exposure and human health has been widely investigated. The Lancet Commission on Pollution and Health estimated that ambient PM causes 4.2 million of deaths each year worldwide (Landrigan et al. 2017) on the basis of data produced for the Global Burden of Disease study (Cohen et al. 2017).

An extensive body of epidemiological literature showed adverse effects of long-term exposure to on mortality and morbidity (e.g., Beelen et al. 2014a; Cesaroni et al. 2014; Di et al. 2017; Hoek et al. 2013) in population-based cohort studies. These studies are characterized by recruitment and follow-up of subjects until the occurrence of the study outcome or censoring. Exposure to air pollutants is usually estimated for each subject at the residential address, either at baseline or as a time-varying annual exposure, and the association with the study outcome is investigated by contrasting individual exposures over space (or space and year), while controlling for potential individual-level or area-level risk factors. One of the major criticisms usually directed to cohort studies is the possibility of residual confounding due to unmeasured or mismeasured confounders (Moolgavkar et al. 2017). Moreover, some other issues are related to standard observational studies. Cohorts are population samples that are not always representative of the whole population. For example, both the American Cancer Society (ACS) cohort and the Nurses’ Health Study (NHS) cohort examined populations with considerably higher levels of education than average (Pope et al. 1995; Puett et al. 2009) or considered only city dwellers (Lepeule et al. 2012) (Jerrett et al. 2013; Krewski 2009).

The need to develop studies that are able to account for measured or unmeasured confounders in the design phase, rather than in the analysis phase, has become increasingly apparent in the last few years. Estimation of a causal effect would require, in principle, the definition of a counterfactual framework, where the study outcomes under alternative exposure scenarios are compared (Hernán and Robins 2019). However, in observational studies, only one exposure distribution is measured (or estimated) for a population for a given time, and causal inference methods try to find the best surrogates for alternative (unobserved) exposure distributions (Rubin 1991). The application of propensity score models suffers from the same concerns about unmeasured confounders (Stürmer et al. 2005), suggesting the importance of alternative study designs that address omitted confounders. The difference-in-differences approach could be a suitable methods to control confounding by design (Card and Krueger 1994). A study conducted in the United States assessed the long-term effects of fine particulate matter with aerodynamic diameter () on mortality in New Jersey (Wang et al. 2016) using a modified version of the difference-in-differences approach by comparing annual changes in mortality rates around census tract–specific means with concurrent annual changes in average exposures within the same census tract. A similar approach was used by Kioumourtzoglou et al. in 207 cities of the northern United States (Kioumourtzoglou et al. 2016).

Another limitation common to most studies investigating air pollution health effects is their focus on urban areas, as the availability of environmental data has been usually confined to the main cities. However, rural and suburban areas represent important settings where large portions of the population live. These areas, despite their lower air pollution concentrations, are influenced by different emission sources (Viana et al. 2008), such as biomass combustion, agricultural emissions, and dust from natural sources (Chen et al. 2017; Pey et al. 2013). Established evidence now shows that low concentrations of PM are related to adverse health outcomes (Makar et al. 2017), that fine and coarse particles (PM between 2.5 and ) display comparable short-term health effects (Zanobetti and Schwartz 2009), and that PM from emission sources other than vehicular traffic (such as desert dust advections and biomass burning) are also related to increased mortality and morbidity (Kim et al. 2013; Stafoggia et al. 2015). In the last few years, data from satellite observations have been made available providing relevant information to supplement ground-level measurements (Kloog et al. 2011, 2014). These data have been used recently in Italy to estimate daily concentration at spatial resolution for the period 2006–2012 (Stafoggia et al. 2017).

The objective of this study was to assess the association between long-term exposure to and cause-specific mortality (nonaccidental, cardiovascular, and respiratory) in the Latium region (central Italy), in the 2006–2012 period by replicating the difference-in-differences approach developed by Wang et al. (Wang et al. 2016). In addition, we aimed to evaluate differential effects of PM on cause-specific mortality in urban, suburban, and rural areas of the region.

Methods

Study Area

The Latium region is in central Italy; it is wide, is divided into 378 municipalities, and had a total population of 5,304,778 inhabitants in 2006 (Figure 1). Rome is the major conurbation of the region, with 2,547,677 residents and an area of . For the purposes of this study, Rome was further divided into 155 urban zones (Figure 2). For each municipality of Latium and the urban zone of Rome (a total of 532 units of observations), we collected annual population data from the National Institute of Statistics (ISTAT) (www.demo.istat.it), which we used as denominators in the statistical analyses described below.

Figure 1.

Population size and concentration in 378 municipalities of the Latium Region during the study period. The population size is reported for the year 2006, and the concentration is the average in the whole period.

Figure 2.

Population size and concentration in 155 urban zones of Rome during the study period. The population size is reported for the year 2006, and the concentration is the average in the whole period.

Health Data

We collected information on mortality data from the Regional Register of Causes of Deaths. For each municipality of the Latium region and each urban zone of Rome, we calculated annual counts of deaths (age ) from nonaccidental (WHO 1978) (International Classification of Diseases Ninth Revision (ICD-9; all codes: 0–799), cardiovascular (CVD; codes 390–459), and respiratory causes (Resp; codes 460–519) during the study period (2006–2012).

Exposure Assessment

Daily concentrations were estimated for each municipality and urban zone and each year on the basis of a spatiotemporal model developed for the entire nation of Italy (Stafoggia et al. 2017). Briefly, we collected multiple sources of data on spatial and spatiotemporal parameters, such as land use and meteorological for each day in 2006–2012 and each square kilometer of Italy, and developed a three-stage model to predict daily values for each grid cell and day by calibrating the satellite aerosol optical depth, land use, and meteorological terms to monitors. We then computed annual for each spatial unit of interest (municipality for the Latium region, urban zone for Rome) by averaging the predicted daily estimates over all grid cells intersecting the spatial units and weighting for the area of intersection (Figure S1 shows the grid cells and the units of the study).

We obtained daily estimates of air temperature by the ERA-Interim reanalysis project (Dee et al. 2011), released by the European Centre for Medium-Range Weather Forecasts (ECMWF). We collected data about air temperature estimated at a height at 0000 hours (12:00 P.M.) and 12000 hours (12:00 A.M.) with spatial resolution of 0.125°×0.125° (approximately ). We calculated a daily mean value of temperature by averaging the two daily retrievals. For the purposes of the study, we considered mean summer temperature (from June to August), mean winter temperature (from December to February), and their standard deviations (SDs). For each spatial unit, mean (and SDs) temperatures were computed by averaging the daily values over all grid cells within the spatial units.

Statistical Analysis

We applied a variant of the difference-in-differences approach to assess the relationship between annual and cause-specific mortality.

The essence of the design is that differences in concentrations across time in a given location (municipality or zone) are related to differences in rates of death in the same municipality or zone so that the role of temporally stable individual and behavioral factors (including unmeasured ones) are canceled out because the comparisons are occurring within the same populations (Card and Krueger 1994). Temporally variable confounders are controlled by contrasting this relation of difference across municipalities or zones with the relation of differences in time in other municipalities (with different exposure differences), thereby controlling for time trends in confounders that are similar across the study region. The standard version of the difference-in-differences approach was based on the comparison between two different areas in two different years (Card and Krueger 1994). Recently, Wang et al. (2016) developed a variant of that study design. The method has been generalized to the case of multiple spatial units and time periods (Wang et al. 2016), by introducing in the statistical model indicator variables for each spatial unit and each time period to remove confounding by all known and unknown factors varying across areas (but considered fixed in time) and over years (but homogenous across space), respectively. Potential confounders, thus, are those factors that: a) display different time variations across the study areas, and b) are correlated with the exposure (annual fluctuations of around the spatial unit average in this study). In the main analysis, we assumed that the only variable with such characteristics might be air temperature, as it seems implausible that, for example, variations from year to year in smoking rates in a municipality or zone around the overall mean for that area and the annual fluctuations in Latium as a whole are correlated with variations from year to year in particle concentrations in that area, around its area specific mean and the overall year-to-year difference in Latium. Therefore, we fit the following model:

where:

represents the number of deaths in spatial-unit (378 municipalities of Latium and 155 urban zones of Rome), year (7-y period, 2006–2012);

is the annual mean concentration of in unit , time ;

is an offset term which represents person-years at risk in unit , time ;

is a dummy variable for each spatial unit ;

is a dummy variable for each year.

Furthermore:

is the intercept term;

are regression coefficients adjusting for confounding induced by factors varying across spatial units (and considered fixed in time) (), and over time (homogeneously across the study area) ();

, , , and are the regression coefficients for the effects of mean summer and winter temperatures and their standard deviations, respectively; and

represents the effect of .

We used the conditional Poisson regression models to analyze the association above described, using the “gnm” package (Turner and Firth 2015). Briefly, we took advantage of conditional Poisson models parameters conditioning on spatial units and “eliminating” (from “eliminate” option of gnm models) the estimates of the variables that did not contribute to the likelihood. Moreover, conditional Poisson models were perfectly comparable with unconditional ones to account for overdispersion and autocorrelation of time-series data (Armstrong et al. 2014).

Because the second objective of our study was disentangling PM effects across different spatial units, and particularly outside of large urban zones cores, we analyzed the whole Latium region and Rome in different ways. First, we analyzed all municipalities of Latium including and excluding Rome as a single unit; then, we focused only on Rome by considering its 155 urban zones; third, we classified municipalities according to their population [1,000–5,000 (194 municipalities); 5,000–14,000 (123); 14,000–40,000 (108); inhabitants (31)], and provided PM effects for each group of municipalities. The choice of the cutoff points was made considering the distribution of city-specific populations in the whole Latium region. We decided to apply cutoff points to have a homogeneous representation of the municipalities of Latium in each class.

As sensitivity analysis, we tested the effect of the same year and lagged exposure (lag 0–1) averaging the concentrations in the current and 1 y before. We modeled the summer and winter temperatures with natural splines with 3 and 2 degrees of freedom, respectively, to check the nonlinearity of the association between temperatures and mortality (data not shown).

Finally, we provided alternative approaches for confounding adjustment. Specifically, we assumed that variables might exist that displayed different linear time trends across areas, possibly co-varying with air pollution levels. To adjust for them, we applied mixed models by adding a random intercept for municipality and a random slope by year. Similarly, we ran a fixed effects model where we inserted an interaction term between year (as linear term) and municipality, aimed at capturing residual confounding from other spatiotemporal covariates not accounted for in the main approach.

All results are expressed as percent increases of risk (IR%), and relative 95% confidence intervals (95% CI) per increase in annual . All statistical analyses were conducted using R software (version 3.1.2; R Development Core Team), specifically the package “gnm” (Turner and Firth 2015) for the main analysis and “lme4” (Bates et al. 2015) for the mixed models.

Results

In the period 2006–2012, we observed deaths due to three causes: 347,699 nonaccidental; 92,787 cardiovascular; and 16,509 respiratory () in the entire Latium region. Corresponding figures for Rome were 178,557; 47,437; and 8,772 deaths from those three causes. Table 1 shows the distribution of deaths, and temperature variables in the whole period both in the municipalities of Latium and in the urban zones of Rome. The average concentration during the study period across all the study units was in Latium and in Rome, with interquartile ranges (IQRs) of 6.9 and respectively. The spatial distribution of population data (year 2006) and concentrations (mean of 2006–2012) are reported in Figure 1 (for Latium) and Figure 2 (for Rome). Population size was represented as total number of residents in the 2006 because no significant changes occurred during the subsequent six years.

Table 1

Environmental and mortality data in the Latium region (378 municipalities) and in the urban zones of Rome (155 units) over the period 2006–2012.

Latium region	Mean	SD	Min	Percentiles			Max	IQR*
				25th	50th	75th
Cause-specific mortality
Nonaccidental	131.41	1317.54	0	5	15	38	26,987	33
Cardiovascular	35.07	349.64	0	1	5	12	7,032	11
Respiratory	6.24	64.84	0	0	0	2	1,354	2
Environmental data
PM_10 (μg/m^3)	21.91	4.87	9.11	18.62	22.25	25.50	33.81	6.88
Mean winter Temperature (°C)	11.91	3.25	5.43	8.63	13.69	14.73	17.19	6.11
Mean Summer Temperature (°C)	17.81	3.67	12.10	14.62	15.91	21.59	25.20	6.97
Standard deviation of winter T (°C)	2.73	0.72	1.07	2.12	2.79	3.24	4.25	1.12
Standard deviation of summer T (°C)	3.86	1.12	1.02	2.83	4.33	4.82	5.51	1.99
Urban zones of Rome
Cause-specific mortality
Nonaccidental	134.75	132.67	0	27	97	194	650	167
Cardiovascular	53.79	54.19	0	11	39	77	264	66
Respiratory	8.61	9.09	0	2	6	13	49	4
Environmental data
PM_10 (μg/m^3)	31.66	3.81	22.84	28.67	31.20	34.44	43.09	5.77
Mean winter Temperature (°C)	12.01	2.59	6.99	9.93	13.13	14.30	15.92	4.37
Mean Summer Temperature (°C)	19.10	3.78	13.83	15.76	16.93	23.53	24.87	7.76
Standard deviation of winter T (°C)	0.63	0.45	0	0.09	0.65	1.01	1.73	0.92
Standard deviation of summer T (°C)	0.88	0.40	0	0.58	0.89	1.24	1.82	0.66

*Interquartile range = 75th–25th percentiles.

Table 2 describes the fluctuations around the area-level mean for annual concentrations and cause-specific mortality rates in Latium (378 municipalities) and Rome (155 urban zones) over the period 2006–2012. Figure 3 describes the SD of the annual concentrations over the study period across the units in Latium.

Table 2

Annual variation in area-level concentrations and cause-specific mortality rates in Latium region (378 municipalities) and in Rome (155 urban zones): Absolute changes are reported for , percent changes are reported for mortality rates.

Latium region	Mean	SD	Percentiles			IQR
			25th	50th	75th
Cause-specific mortality
Nonaccidental	0	48.9	−19.1	−1.1	16.7	35.8
Cardiovascular	0	78.1	−39.7	−3.3	24.4	1.29
Respiratory	0	140	−100	−37.51	41.1	141.1
Particulate matter
PM_10 (μg/m^3)	0.00	2.44	−1.91	−0.39	1.99	3.90
Urban zones of Rome
Cause-specific mortality
Nonaccidental	1.33	29.23	−8.11	−0.041	7.80	15.91
Cardiovascular	2.50	37.42	−12.11	−0.17	11.67	23.78
Respiratory	14.23	70.36	−19.96	1.83	29.85	49.81
Particulate matter
PM_10 (μg/m^3)	0	2.96	−2.51	−0.88	2.59	5.1

Figure 3.

Standard deviation of the annual concentrations for each municipality in the whole region over the period 2006–2012.

The results of the main analysis are reported in Table 3. We found positive associations between long-term exposure to and cause-specific mortality in the whole Latium region, with IR% (95% CI) of 0.75% (0.17, 1.34), 0.93% (0.03, 1.83), and 1.42% (, 3.25) per increase in , for nonaccidental, cardiovascular, and respiratory causes, respectively. In Rome, we found a significant effect on nonaccidental mortality [0.55% (0.04, 1.07)], whereas the effect estimates for cardiovascular and respiratory mortality had large confidence intervals. The overall effects in the Latium region were not substantially affected by inclusion or exclusion of Rome.

Table 3

Associations between long-term exposures to environmental variables and cause-specific mortality. Results are expressed as percent increase of risk and relative 95% confidence intervals (CI) per increase of .

Area/cause-specific mortality	Mortality
	IR%	95% CI
Latium Region
Nonaccidental	0.75	0.17	1.34
Cardiovascular	0.93	0.03	1.83
Respiratory	1.42	−0.38	3.25
Latium region without Rome
Nonaccidental	0.57	−0.07	1.22
Cardiovascular	0.59	−0.38	1.57
Respiratory	2.02	0.05	4.04
Rome (155 urbanistic zones)
Nonaccidental	0.53	−0.05	1.12
Cardiovascular	0.22	−0.64	1.08
Respiratory	0.57	−1.43	2.62

Table 4 displays effect estimates according to population size. We observed a tendency of higher risks on nonaccidental mortality in smaller spatial units in comparison with larger ones, ranging from 0.76% (, 1.88) in municipalities with fewer than 5,000 inhabitants to (, 0.90) in cities with more than 40,000 residents (except Rome) (p-value for ). In contrast, there was a tendency toward stronger effects on respiratory mortality in larger municipalities (p-value for ).

Table 4

Associations between long-term exposures to and cause-specific mortality in different type of municipalities (the number of municipalities in each class are reported in brackets) in the Latium region. Results are expressed as percent increase of risk and relative 95% confidence intervals (CIs) per increase of .

Effect modifiers	Mortality
	Nonaccidental				Cardiovascular				Respiratory
Population*	IR%	95% CI		p-int	IR%	95% CI		p-int	IR%	95% CI		p-int
1,000–5,000 (194)	0.76	−0.35	1.88	—	−0.13	−1.63	1.39	—	−0.96	−4.21	2.40	—
5,000–15,000 (123)	0.38	−0.68	1.44	0.286	0.50	−1.03	2.06	0.251	0.99	−2.30	4.40	0.175
15,000–40,000 (108)	0.44	−0.46	1.35	0.302	0.34	−0.95	1.65	0.285	3.57	0.78	6.44	0.007
>40,000 (31)	−0.01	−0.92	0.90	0.103	0.50	−0.87	1.89	0.230	2.98	0.24	5.80	0.015

Size (inhabitants).

Note: p-Int, p-value for interaction.

We also investigated the role of summer and winter temperatures on mortality in the whole Latium region. There was a suggestion of a harmful effect of summer temperature with an IR% (for 1°C increases) equal to 4.74 (, 14.10), 2.98 (, 17.17), and 17.2 (, 51.87) for nonaccidental, cardiovascular, and respiratory mortality, respectively. In contrast, winter temperatures showed no effect on cause-specific mortality.

When we evaluated lagged effects of (lag 0–1) the results were similar in comparison with the main analysis (lag 0) (see Supplemental Material). Finally, modeling the summer and winter temperatures with natural splines with 3 and 2 degrees of freedom, respectively, did not change the effect of on cause-specific mortality (see Supplemental Material).

Table 5 shows the comparison between the main analysis and the sensitivity approach for confounding adjustment (mixed model with random intercepts and slopes, fixed effects model with interaction terms). These models were run for the whole Latium region only. We observed similar estimates for nonaccidental mortality across the three models, ranging from 1.03% (0.49, 1.57) to 0.69% (0.35, 1.04). Cardiovascular mortality displays no longer positive association with in sensitivity analyses, whereas respiratory mortality shows similar effects with the mixed model and higher estimates with the fixed effects model.

Table 5

Associations between long-term exposures to and mortality for nonaccidental, cardiovascular and respiratory causes over the Latium region by three different modeling approaches: Difference-in-differences (base model), fixed effects model, and mixed model. All results are expressed for increase in .

Approach	IR%	95%CI
Difference in differences
Nonaccidental mortality	0.75	0.17	1.34
Cardiovascular mortality	0.93	0.03	1.83
Respiratory mortality	1.42	−0.38	3.25
Fixed effects model
Nonaccidental mortality	1.03	0.49	1.57
Cardiovascular mortality	0.59	−0.44	1.63
Respiratory mortality	3.34	0.73	6.02
Mixed model
Nonaccidental mortality	0.69	0.35	1.04
Cardiovascular mortality	−0.01	−0.45	0.44
Respiratory mortality	0.94	0.83	1.04

Discussion

We investigated the relationship between long-term exposure to and cause-specific mortality in the Latium region during the period 2006–2012. Using a variant of difference-in-differences approach, our results provided evidence of a link between and mortality in the study area, with percent increases of mortality rate equal to (for increase of ) to 0.75%, 0.93%, and 1.42% for nonaccidental, cardiovascular, and respiratory causes, respectively. If we consider the average ratio for the fixed monitors in Rome equal to 0.65, these estimates correspond to a long-term effect of equal to 1.16%, 1.43%, and 2.19%, respectively. In addition, we tested alternative approaches for confounding adjustment in sensitivity analyses, obtaining results similar to results using the main approach.

In the last years, many cohort studies investigated the association between long-term exposure to and cause-specific mortality worldwide. Hoek et al. performed a systematic review in 2013 (Hoek et al. 2013) in which they meta-analyzed 11 studies reporting the long-term effects of on all-cause mortality and 10 studies on cardiovascular mortality. The authors concluded that each increase of in was associated to an excess risk of 6% (4, 8%), and 11% (5, 16%) for nonaccidental and cardiovascular mortality, respectively. No effect was found for respiratory mortality. In Europe, the ESCAPE Project analyzed the role of long-term exposure to air pollution on several health outcomes in 22 European cohorts. The authors estimated the exposure to several pollutants using land use regression models and examined the association by Cox regression analyses. They estimated hazard ratios (HRs) of 1.04 (1.00–1.09), 1.02 (0.92–1.14), and 0.86 (0.67–1.04) for each increase of for nonaccidental, cardiovascular, and nonmalignant respiratory mortality, respectively (Beelen et al. 2014b, 2014a; Dimakopoulou et al. 2014). Their effect estimates were adjusted for environmental and behavioral variables such as smoking status, alcohol consumption, fat intake, and BMI, increasing the importance of the evidence published. Recently, Vodonos et al. (2018) in a systematic review found that for all-cause all-age mortality, a increase in was associated with a 1.29% increase in all-age all-cause mortality (95% CI: 1.09, 1.50) at a mean exposure of , which decreased to 1.03% (95% CI: 0.97, 1.11) at a mean exposure of (the mean level across all studies), and to 0.82% (95% CI: 0.52, 1.12) at Vodonos et al. 2018). The estimates for cardiovascular and respiratory mortality were 1.46% (95% CI: 1.25, 1.67), and 1.13% (95% CI: 0.85, 1.41), respectively, at a mean exposure of .

Although quantitative comparisons of the effects of and are not entirely appropriate, all the results indicate a clear role of PM exposure on mortality for nonaccidental, cardiovascular, and respiratory causes. We found results very similar to those reported in a recent systematic review conducted by Vodonos et al. It should be noted that there are differences between the methods used in our analysis and those used in previous cohort studies. Most of the available cohort studies used exposure levels at a fixed time and considered only some of the potential individual confounders. In our case, we controlled for each slow-changing variable in time, such smoking status, BMI, socioeconomic status, etc., by design. We also controlled for variables that changed from year to year similarly across the region by design, thus allowing us to control also for unmeasured confounders. Therefore, we have assurance that most unmeasured potential confounders of the previous studies do not confound our results, under the assumption that variations in other variables around location specific means and region specific time trends are not correlated with similar variations in .

Recent studies applied a similar design to investigate the association between particles and mortality. Wang et al. used a variant of difference-in-differences method to assess the effect of long-term exposure to on natural causes mortality in New Jersey during the 2004–2009 period (Wang et al. 2016). They assigned an average exposure to for each year and for each ZIP code in New Jersey using satellite-based data. They found an IR% of 15.5 (0.8, 32.3) for increase of pollutant in the entire population of New Jersey. Although the methods used in that study were similar to our approach, the effect estimates were slightly higher than those that we found. The higher estimates found by Wang et al. could be explained by the pollutant considered. In fact, it is reasonable to expect that would show a lower effect in comparison with that of exposure. However, some other untestable hypotheses could play a role, such as the different study period and study area, different composition of particles, or the economic crisis that struck Europe in 2008 leading to a progressive weakness of the study population (Karanikolos et al. 2013). Kioumourtzoglou et al. used a similar approach to evaluate the effect of long-term exposure to in 207 cities in the United States (Kioumourtzoglou et al. 2016). They applied a variant of difference-in-differences method analyzing Medicare enrollees (). The authors assigned to each subject of the same city the same annual exposure to obtained by city-specific monitors. The results showed a strong positive association with a HR of 1.19 (95% CI: 1.11, 1.28) per increase of related to total mortality. In this case, the comparison with our study is not as immediate because they combined a city-specific exposure with individual information using the counting process extension of the proportional hazards model by Andersen and Gill (Andersen and Gill 1982) instead of the Poisson regressions at city levels applied in our study. However, the two methods are statistically comparable as demonstrated in several frameworks (Peters et al. 2006). In addition, the results displayed similar trends with effect estimates in line with other epidemiological evidence.

Our results suggested that in towns with fewer than 5,000 inhabitants the effect of is higher for nonaccidental mortality and null for respiratory mortality. For bigger cities, we observed an opposite trend with higher effect estimates for respiratory mortality in comparison with other causes. A possible explanation of these results might be chance alone or a different mixture composition of ambient particles in rural areas in comparison with urban areas. Some studies indicate that PM composition in rural areas is characterized mostly by crustal [aluminum (Al), silicon (Si), calcium (Ca), iron (Fe)], and sea-salt components [sodium (Na), chlorine (Cl), magnesium (Mg)], whereas in urban settings, vehicular and industrial constituents represent the greater part of the total composition (Götschi et al. 2005; Viana et al. 2008). Epidemiological evidence showed that some components displayed different effects on cause-specific mortality (Atkinson et al. 2015). For example, Si was associated with nonaccidental mortality in the 15 cohorts of the ESCAPE framework (Beelen et al. 2015). Although these pieces of evidence supported our findings, we were not able to separate the component-specific effects in our study; therefore, this lack of information is an untestable hypothesis in our case.

Our main approach assumed that time trends might be captured adequately by dummies for years, under the assumption that these were uniform in the study region. In the sensitivity analysis, we substantially relaxed such an assumption by adding interactions between linear time trends and municipalities. Results were very similar for natural mortality and comparable for cardiovascular and respiratory outcomes. Although the sensitivity analysis points against substantial residual confounding, we still cannot entirely rule out the existence of omitted covariates, changing differently across years from area to area in a nonlinear way, partially biasing our estimates.

To our knowledge, this investigation is one of the few studies that attempted to estimate the link between long-term exposure to and cause-specific mortality using a difference-in-differences approach. Our approach presents some limitations that we must mention. Apart from temperature, we did not adjust for other spatiotemporal predictors potentially confounding the association under investigation. We could not conceive of any variable displaying different time trends across municipalities and co-varying with annual PM concentrations, other than temperature (which we controlled for) and influenza epidemic. Especially relevant for respiratory mortality, we are not able to control for influenza epidemics, which might vary both spatially and temporally during the study period and across the study area. In addition, in the extreme case, socioeconomic trends, smoking rates, etc., could be considered under the (unrealistic) hypothesis that they varied differently across the study areas and were also related to annual fluctuations. To critically evaluate our findings, we computed the E-value proposed by VanderWeele and Ding in 2017 (VanderWeele and Ding 2017) to evaluate the effect of potential unmeasured confounding. Briefly, E-value is defined by the authors as “the minimum strength of association, on the risk ratio scale, that an unmeasured confounder would need to have with both the treatment and outcome, conditional on the measured covariates, to fully explain away a specific treatment–outcome association.” E-value could be easily calculated by the formula: “ ).” In our case, we estimated an E-value for nonaccidental mortality of 1.39 for the effect estimate and 1.27 for the lower limit of the CI related to the relative risk. This approach means that an unmeasured confounder could explain our estimates away if it were associated with both exposure and outcome with a risk ratio of 1.39. Similarly, it should be associated with both exposure and outcome with a risk ratio of 1.27 to move the CI to include the unity. We consider highly unlikely that such an unmeasured (spatiotemporal) confounder exists. Considering exposure assessment, we know that exposure misclassification can occur from assigning the same yearly averaged in each spatial unit for all residents. In addition, we used average exposure for each city introducing a possible exposure misclassification with relative Berksonian error, which should not bias our estimates but would increase the CIs. Finally, it should be considered that our approach may have a more limited statistical power in comparison with traditional cohort analyses as we have only a limited fluctuation in the concentration around its overall mean over the study period; this lack of power is reflected in the large CIs of the effect estimate for respiratory mortality.

Conclusion

Exposure to ambient air particles represent an important risk factor for human health worldwide; however, the epidemiological findings should be reinforced with alternative study designs and modeling approaches. Our analysis suggests an effect of on cause-specific mortality in the Latium region that is unlikely to be explained by confounding. PM affects mortality not only in urban settings but also in suburban and rural areas. Finally, the difference-in-differences approach might be considered as a good alternative method to estimate an association between long-term exposure to air pollution and mortality.

Click here for additional data file.

Click here for additional data file.