Research Volume 123  2015
CommutingAdjusted ShortTerm Health Impact Assessment of Airborne Fine Particles with Uncertainty Quantification via Monte Carlo Simulation
Michela Baccini,^{1}^{,}^{2} Laura Grisotto,^{1} Dolores Catelan,^{1}^{,}^{2} Dario Consonni,^{3} Pier Alberto Bertazzi,^{3}^{,}^{4} and Annibale Biggeri^{1}^{,}^{2}

Background: Exposure to air pollution is associated with a shortterm increase in mortality, and this field has begun to focus on health impact assessment.
Objectives: Our aim was to estimate the impact of PM_{10} on mortality within 2 days from the exposure in the Italian region of Lombardy for the year 2007, at the municipality level, examining exposure entailed by daily intermunicipality commuting and accounting for uncertainty propagation.
Methods: We combined data from different sources to derive probabilistic distributions for all input quantities used to calculate attributable deaths (mortality rates, PM_{10} concentrations, estimated PM_{10} effects, and commuting flows) and applied a Monte Carlo procedure to propagate uncertainty and sample the distribution of attributable deaths for each municipality.
Results: We estimated that annual average PM_{10} concentrations above the World Health Organizationrecommended threshold of 20 μg/m^{3} were responsible for 865 shortterm deaths (80% credibility interval: 475, 1,401), 26% of which were attributable to PM_{10} above the European Union limit of 40 μg/m^{3}. Reducing annual average PM_{10} concentrations > 20 μg/m^{3} by 20% would have reduced the number of attributable deaths by 36%. The largest estimated impacts were along the basin of the Po River and in the largest cities. Commuting contributed to the spatial distribution of the estimated impact.
Conclusions: Our estimates, which incorporated uncertainty quantification, indicate that the shortterm impact of PM_{10} on mortality in Lombardy in 2007 was notable, and that reduction in air pollution would have had a substantial beneficial effect on population health. Using commuting data helped to identify critical areas for prioritizing intervention.

Citation: Baccini M, Grisotto L, Catelan D, Consonni D, Bertazzi PA, Biggeri A. 2015. Commutingadjusted shortterm health impact assessment of airborne fine particles with uncertainty quantification via Monte Carlo simulation. Environ Health Perspect 123:27–33; http://dx.doi.org/10.1289/ehp.1408218
Address correspondence to M. Baccini, Department of Statistics, Informatics and Applications “G. Parenti,” viale Morgagni 5950134 Firenze, Italy. Telephone: 39 320 4209085. Email: baccini@disia.unifi.it
We thank the Air Quality Unit of the Lombardy Regional Agency for the Environment and the Health Directorate of the Region of Lombardy Government for providing the data that made this work possible.
This work was supported by the Regional Government of Lombardy (grant VIII/10462) and the Ministry of University and Scientiﬁc Research (grant PRIN 20072S2HT8).
The authors declare they have no actual or potential competing financial interests.
Received: 4 February 2014
Accepted: 15 October 2014
Advance Publication: 17 October 2014
Final Publication: 1 January 2015  Supplemental Material (435 KB) PDF
Introduction
The role of air pollution in shortterm and longterm disease causation is widely recognized [International Agency for Research on Cancer 2013; World Health Organization (WHO) 2013]; many health impact assessments have been published, and others are ongoing (Baccini et al. 2011; Ballester et al. 2008; Heal et al. 2013; Künzli et al. 2000; Li et al. 2010; Martuzzi et al. 2006; Medina et al. 2004). Estimates of the shortterm impact of air pollution on mortality (i.e., the impact within a limited number of days from the exposure, usually less than a week) have advantages over estimates of longterm impact for assessing the effectiveness of emissions reduction policies, because they are less influenced by latency time, cumulative exposure, and demographic dynamics (Flachs et al. 2013). Moreover, because the prevalence of exposure is high, populationlevel effects on mortality may be large, despite the small relative risks measured by the epidemiological models.
The treatment of uncertainty arising from different sources is a major point of concern when evaluating the shortterm impact of air pollution (Knol et al. 2009). If air pollution level and total number of events are known, only sampling variability around the effect estimate has to be considered (Baccini et al. 2011; Orru et al. 2009). However, exposure and mortality levels are often predicted rather than observed—for example, when information for a specific area and/or time period is lacking. In these situations, approaches that deal simultaneously with different sources of uncertainty are needed.
Uncertainty can be accounted for by using Monte Carlo (MC) techniques that generate samples from unknown output probability distributions based on repeated random sampling from independent known input distributions, particularly when a closed form of the output distribution is impossible or difficult to obtain, as in complex nonlinear or multidimensional problems (Kumamoto and Henley 1996). This makes it possible to incorporate in the output the uncertainty affecting the inputs. MC methods for propagating uncertainty have not been widely applied in the context of environmental health impact assessment. MesaFrias et al. (2013) conducted a systematic review of methods used to quantify uncertainty in this field: Only 19 studies between 2000 and 2012 addressed uncertainty; of these, only 14 adopted probabilistic approaches.
A second concern is the assumption of a static (noncommuting) population. Few studies have investigated the influence of daily commuting on population exposure to air pollution (e.g., Berrocal et al. 2011; Calder et al. 2008). Recently, the European Topic Centre on Air Pollution and Climate Change Mitigation evaluated the contribution of population commuting on exposure to particulate matter in European urban areas and stated that exposure estimates that include commuting may be larger than those based on a static assumption because people usually move from lower to higherpolluted areas (Larssen et al. 2012).
Our aim in this paper was to estimate the shortterm impact of high concentrations of particles up to 10 μm in diameter (PM_{10}) on mortality in the Italian region of Lombardy in 2007 at the municipality level, considering commuting and accounting for uncertainties that arise in the calculation of attributable deaths (AD) via MC simulation. We focused on PM_{10} because the evidence for the causal mechanism between exposure and health damage is more consolidated for this air pollutant than for others (Anderson et al. 2012).
The region of Lombardy. Lombardy is a 23,865km^{2} region in northwestern Italy, which can be geographically and economically divided into three zones: the mountain range of the Alps, the sloping foothills, and the highly industrialized and populated basin of the Po River. The latter is characterized by a high level of air pollution due to frequent thermal inversion, with pollution being trapped close to the ground level. In 2007, Lombardy had 9.8 million inhabitants living in 1,546 municipalities located among 11 provinces (Figure 1).
Figure 1 – Subdivision of Lombardy by province in 2007. The provincial capitals and the cities of Busto Arsizio and Vigevano are indicated in white. Grayscale expresses altitude, with darker tones indicating mountain areas. Reproduced from Baccini et al. (2011), with permission of Oxford University Press.
Mobility in Lombardy is substantial, with the highest percentage in Italy (53%) of residents commuting daily to workplaces or schools, half of them out of their residence municipality (Regione Lombardia 2004).
Levels of air pollution vary throughout the region. Between 2003 and 2006 the annual average concentration of PM_{10} in Milan, the capital of the region (1,299,633 inhabitants in 2007), was 52.5 μg/m^{3}, with 95% of daily concentrations > 20 μg/m^{3}, whereas the average across the other most densely populated areas of Lombardy was 45.4 μg/m^{3} (Baccini et al. 2011).
Methods
Data and input distributions. To derive the probabilistic distributions for input data used in the MC procedure, we used information on total mortality, PM_{10} concentrations, PM_{10} effects, and intermunicipality commuting flows, arising from different sources.
Smoothed crude mortality rates. For each municipality, the average number of deaths during the period 2000–2004 from the Regional Mortality Register and the number of inhabitants at the 2001 Italian population census were available [Italian National Institute of Statistics (Istat) 2014].
Because risk estimates from small areas suffer from substantial variability, we applied a method, widely used in disease mapping, to smooth crude mortality rates at the municipality level (Besag and Kooperberg 1995; Besag et al. 1991). Specifically, we specified a Bayesian model that accounted for structured and unstructured spatial variability, according to the Besag, York, and Mollie’s (BYM) proposal (Besag et al. 1991). Let deaths_{i} be the total number of deaths during the period 2000–2004 in municipality i (i = 1, 2,…,1,546). We assumed that deaths_{i} followed a Poisson distribution with mean given by the product of the death rate r_{i} and the population size D_{i}, which we estimated to be five times the population size in 2001. Moreover, we specified a randomeffect loglinear model on r_{i}:
log(r_{i}) = u_{i} + v_{i}, [1]
where u_{i} and v_{i} were independent random terms representing the unstructured spatial variability component and the structured spatial variability component, respectively. We assumed that u_{i} were independent and normally distributed with mean 0 and variance τ_{u}^{2}, and that v_{i} followed an intrinsic conditional autoregressive (ICAR) model. In other words, for each S_{i}, the set of the n_{i} municipalities adjacent to municipality i, we assumed the following conditional distribution for v_{i}:
v_{i}  v_{j} Є S_{i} ~ N(v̄_{i}, τ_{v}^{2}/n_{i}), [2]
where v̄_{i} was the mean of v_{i} for the municipalities belonging to S_{i}, and τ_{v}^{2}/n_{i} was their conditional variance (Besag and Kooperberg 1995). Through the random terms u_{i} and v_{i}, the BYM model shrinks the crude rate estimates toward both the local and the general mean.
After defining noninformative inverse gamma priors on the variance parameters, we used Monte Carlo Markov Chain (MCMC) methods to obtain a sample from the joint posterior distribution of the smoothed crude rates r_{i}. The smoothed crude rates reflect spatial variability of mortality, but at the same time are more stable than the observed crude rates, which could be very unstable for the smallest cities. This analysis was performed using WinBUGS 1.4.3 (Lunn et al. 2000).
Predicted annual average concentrations of PM_{10.} Data about annual concentrations of PM_{10} were available from two sources: a Eulerian photochemical model developed by the Regional Agency of Environmental Protection (Milan, Italy) that accounted for transport, chemical conversion, and deposition of atmospheric pollutants (Silibello et al. 2008), and the regional monitoring network for air quality control. The Eulerian model provided predictions at the level of 4 × 4 km grid cells that covered the region; predictions were in the form of cell averages and, because of the deterministic nature of the model, did not convey any information about the inherent uncertainty around the prediction process. Data from monitors were sparsely collected from 58 monitors in the region. To obtain an estimate of the uncertainty around the predictions while retaining coverage of the entire region, we used data from both sources to specify a Bayesian geostatistical model, in the form of universal kriging, and used MCMC methods to obtain a sample from the joint posterior predictive distribution of the annual averages at the municipality level.
We assumed that the vector of the log annual average concentrations from the 58 monitors, log(x), was a realization from a multivariate normal distribution with a mean vector μ depending on z _{s}, the vector of the annual average concentrations predicted by the Eulerian model for the grid cells where the 58 monitors were located, and variance covariance matrix Σ = {Σ_{ij}} expressing the spatial correlation structure as a function of the distance between pairs of monitors. Specifically, μ = α + βz _{s}, where α and β were regression coefficients, and Σ_{ij} = σ^{2}f(d_{ij};κ,φ), where f(d_{ij};κ,φ) = exp(–φ d_{ij} κ), σ^{2} was a variance parameter, d_{ij} was the known distance between monitor i and monitor j, and φ and κ were positive parameters (κ ≤ 2) that controlled the rate of decline of the correlation by distance and the amount of spatial smoothing, respectively. Priors for φ and κ were chosen to produce zero correlation between any pair of points at the maximum distance (250 km) and one correlation at the minimum distance (3 km); uninformative priors were specified for all other parameters in the model.
We applied MCMC methods to obtain a sample from the joint posterior distribution of the model parameters and, using as covariate values z the predictions from the Eulerian model in each 4 × 4 km cell, we obtained a sample from the joint posterior predictive distribution of the annual average levels of PM_{10} at the cell centroids of the entire grid. Finally, we derived the joint posterior predictive distribution for the annual average concentration of PM_{10} for each municipality by integrating over these predicted cell values.
We confirmed the validity of the geostatistical model through leaveoneout crossvalidation. Root mean square error (RMSE) and fractional bias (FB) values, evaluated on the log scale, were small (RMSE = 0.010; FB = 0.035%), indicating good predictive performance of the model. This analysis was performed using WinBUGS 1.4.3 (Lunn et al. 2000).
Estimates of PM_{10} effect. Estimates of PM_{10} effect were derived from the Bayesian random effects metaanalysis in Baccini et al. (2011), which combined estimates of the effect of PM_{10} during the same day and previous day (lag 0–1) on mortality in 13 areas in Lombardy during 2003–2006. The 13 areas were the entire administrative agricultural district of Lodi; the provincial capitals of Bergamo, Brescia, Como, Cremona, Lecco, Mantova, Milano, Pavia, Sondrio, and Varese; and the large municipalities of Busto Arsizio and Vigevano. For the present impact evaluation we used the joint posterior distribution of the cityspecific effects for the 13 areas included in the metaanalysis, and the posterior predictive distribution of a generic cityspecific effect for all other municipalities (Riley et al. 2011). A sample from these distributions was obtained using MCMC methods.
Probability of commuting. Data on regular commuting flows over the region were derived from the 2001 Italian population census (see Supplemental Material, Figure S1). For each municipality i, we knew the absolute number of inhabitants that regularly commuted to school or work in municipality j (c_{ij}, j ≠ i). Specifying a vague a priori beta distribution on the commuting probability from municipality i to j (p_{ij}), and a binomial likelihood for c_{ij} with binomial denominator equal to the population of i at the 2001 census (pop_{i}^{2001}), the resulting posterior distribution for p_{ij} was a beta distribution with parameters 1 + c_{ij} and 1 + pop_{i}^{2001} – c_{ij}. We assumed that commuting probabilities were mutually independent. Because of the small size of the commuting probabilities, no constraint was needed in practice to assure that, once i was fixed, the sum of all p_{ij} (where j ≠ i) was < 1.
AD calculation. Let us assume that the exposure is homogeneous within the municipality; in the absence of intermunicipality commuting, the fraction of deaths among residents of municipality i that are attributable to PM_{10} depends only on the exposure–response curve and the concentration of PM_{10} in municipality i. Then, if one assumes linearity on a log scale, the deaths attributable to exposures exceeding the threshold x_{0} (AD _{i}) can be calculated for each municipality as
AD_{i} = y_{i} – y_{i}/exp[β _{i} × (x_{i} – x_{0}) × I(x_{i} > x_{0})], [3]
where y_{i} is the total number of deaths in municipality i, β _{i} is the estimated PM_{10} effect in i, x_{i} is the annual average PM_{10} concentration in i, and I(x_{i} > x_{0}) is a indicator function with I = 1 for x_{i} > x_{0} and I = 0 otherwise. Usually y_{i} is assumed to follow a Poisson distribution with mean μ_{i} = pop_{i} × r_{i}, where r_{i} is the crude mortality rate and pop_{i} represents the personyears at risk in municipality i, which roughly correspond to the total number of inhabitants in i during the year of interest.
In commuting, individuals can be exposed to different levels of PM_{10} during the day, so that Equation 3 is no longer valid. Therefore, for each municipality i we calculated three different quantities: A_{i}, the number of deaths among residents of i attributable to exposure in i; B_{i}, the number of deaths among residents of i attributable to exposure in j (where j ≠ i); and C_{i}, the number of deaths among nonresidents of i attributable to exposure in i.
Assuming that, on average, regular commuters spend onethird of their time in municipality j (where they work or study), and twothirds of their time in municipality i (where they live),
A_{i} = y_{i}^{S} – y_{i}^{S}/exp[β_{i} × (x_{i} – x_{0}) × I(x_{i} > x_{0})], [4]
where y_{i}^{S} is the total number of deaths among the residents in municipality i associated with the personyears at risk in i [with the superscript S indicating the static (noncommuting) portion of the population]. To account for intrinsic variability related to random variation of deaths count, we assumed that y_{i}^{S} followed a Poisson distribution, with mean μ_{i}^{S} equal to the product between personyears at risk and the crude mortality rate r_{i}:
μ_{i}^{S} = [pop_{i} – 1/3 × Σ_{j} _{≠}_{i} exit(i,j)] × r_{i}, [5]
where pop_{i} was the intercensual estimate of the population of i in 2007 (Istat 2014), and exit(i,j) was the number of individuals who regularly commuted from i to work or school in j (where j ≠ i). In turn, exit(i,j) was assumed to follow a binomial distribution with a probability of success equal to the probability of commuting from i to j, and the number of trials equal to pop_{i}. Because the probabilities of commuting were small relative to the probability of not commuting, no constraint was needed to avoid negative values of μ_{i}^{S}.
B_{i} is the sum over j, with j ≠ i, of deaths among residents of i attributable to the exposure in j, which were estimated based on the personyears at risk spent in municipality j by residents of i, the crude mortality rate for municipality i (r_{i}), and the attributable risk associated with exposure in j:
B_{i} = Σ _{j} _{≠}_{i} {y_{ij}^{C} – y_{ij}^{C}/exp[β_{j} × (x_{j} – x_{0}) × I(x_{j} > x_{0})]}, [6]
where y_{ij}^{C} is the total number of deaths among commuters from i to j (with the superscript C denoting the commuting portion of the population), assumed to follow a Poisson distribution with mean μ_{ij}^{C} = 1/3 × exit(i,j) × r_{i}.
Finally, C_{i} is the sum over j, with j ≠ i, of deaths among commuters from j to i attributable to exposure in i, which were estimated based on the personyears at risk spent in municipality i by residents of j, the crude mortality rate for municipality j (r_{j}), and the attributable risk associated with exposure in i:
C_{i} = Σ_{j} _{≠}_{i} {y_{ji}^{C} – y_{ji}^{C}/exp[β_{i} × (x_{i}–x_{0}) × I(x_{i} > x_{0})]}, [7]
where y_{ji}^{C}, the total number of deaths among commuters from j to i, was assumed to follow a Poisson distribution with mean μ_{ji}^{C} = 1/3 × exit(j,i) × r_{j}. In turn, the number of individuals who regularly commuted from j to work or school in i, exit(j,i), was assumed to follow a binomial distribution with the probability of success equal to p_{ji} and the number of trials equal to pop_{j}.
As the formulas show, we assumed that the PM_{10} effect depends on the municipality where the individual is exposed, whereas the mortality rate is that of the municipality where the individual resides.
By combining A_{i}, B_{i}, and C_{i}, we estimated two impact measures:
 AD_{i}^{A} ^{+}^{B} = A_{i} + B_{i}, the number of deaths among residents of i attributable to exposures in the region (i.e., the exposure in i or in other municipalities of the region); and
 AD_{i}^{A} ^{+}^{C} = A_{i} + C_{i}, the number of deaths among residents of the region (i.e., the residents in i or in other municipalities of the region) attributable to exposure in i.
It should be noted that Σ_{i}B_{i} = Σ_{i}C_{i}, which is the total estimated number of deaths due to commutingrelated exposure in the region, and Σ_{i}AD_{i}^{A + B} = Σ_{i}AD_{i}^{A + C}.
Air pollutant reduction scenarios. We estimated AD under different reduction scenarios (RS) corresponding to different definitions of the threshold x_{0}:
 RS0: x_{0} = 20 μg/m^{3}, the WHO Air Quality Guideline threshold for PM_{10} annual average (WHO 2005);
 RS1: x_{0} = 40 μg/m^{3}, the European Union (EU) limit for PM_{10} annual average (EU 2008);
 RS2: x_{0} equal to a reduction of 20% in the observed annual concentration of PM_{10}, provided it is > 20 μg/m^{3};
 RS3: x_{0} equal to a reduction of 20% in the observed annual concentration of PM_{10}, provided it is > 40 μg/m^{3}.
The RS2 and RS3 scenarios define intermediate targets coherent with realistic policies of progressive reduction in PM_{10} concentration until the limits of 40 μg/m^{3} or 20 μg/m^{3} are reached (EU 2008).
Under RS0 and RS1, for the 11 provinces and their capitals, we also calculated the attributable community rates (ACRs), which allow impact comparison among different populations (Wacholder 2005). ACR is defined as the difference between overall crude risk in the population and risk in the unexposed individuals; in our context, ACR was estimated as the ratio between AD and population size.
Monte Carlo method. After obtaining 1,000 independent realizations from the distribution of each input quantity (“Data and input distributions” section; see also Supplemental Material, Table S1, for a summary description of the input distributions), we combined them accordingly to Equations 4–7 to obtain 1,000 values from the joint posterior distribution of all quantities of interest (A_{i}, B_{i}, C_{i}, AD_{i}^{A + B}, AD_{i}^{A + C}; i = 1,2,…,1,546).
MC simulation was performed using R statistical software (R Core Team; http://www.Rproject.org/). The code is available at http://www.biostatistica.net/sites/MC_mobility/prog_data.zip.
Evaluating the role of commuting. For each municipality we calculated the logarithm of the ratio between the posterior median of C_{i} and the posterior median of B_{i}, as a measure of the balance between AD “exported” (deaths among individuals residing elsewhere in the region attributable to the air pollution level in the municipality i) and AD “imported” (deaths among the residents in the municipality i, attributable to their commutingrelated exposure elsewhere in the region). Negative values of the ratio indicated that the imported AD exceeded the exported AD. Positive values indicated that the exported AD exceeded the imported AD.
The problem of sampling negative values of air pollutant effect. Impact measures are appropriate when dealing with risk factors that are causally related to the outcome, and under these circumstances, negative estimates are structurally impossible. Therefore, we set AD = 0 whenever a negative value for the PM_{10} effect was sampled during the MC procedure, because we were interested in the AD distribution conditional on rejection of the null hypothesis of no PM_{10} effect in favor of the unilateral alternative hypothesis that exposure increased the risk of death.
To obtain conservative estimates avoiding an overestimation of impact, we summarized AD distributions using the median instead of the mean. We also provided the 80% credibility interval (CrI), defined as the 10th and 90th percentiles of the posterior distribution, and the posterior probability of a positive number of AD.
Results
The highest annual average concentrations of PM_{10} were in the central area of the region (Figure 2A). The air pollution level in the municipalities around Milano and Brescia was clearly > 40 μg/m^{3}. Milan and the neighboring municipalities were characterized by annual average concentrations > 50 μg/m^{3}.
Figure 2 – Posterior means of annual average PM_{10} concentrations (μg/m^{3}) (A) and smoothed crude mortality rates (per 100,000 residents) (B); posterior medians of attributable deaths among residents (C) and posterior probabilities (%) of a nonnull impact (D) under the RS0 scenario (for annual average PM_{10} > 20 μg/m^{3}), by municipality.
The smoothed crude mortality rates reflected the age distribution of the resident population, with higher estimated rates in municipalities with a larger percentage of elderly people (Figure 2B). Mortality estimates were lower in the central part of the region, except for the three largest cities: Milan, Bergamo, and Brescia.
Estimates of the shortterm effect of PM_{10} were similar among cities except for Milan, where the estimated effect was twice the overall metaanalytic estimate (Baccini et al. 2011). The probability of a negative value for the percent variation was close to 0 for Milan, but > 30% for Brescia and Lecco, with the other municipalities in between (see Supplemental Material, Figure S2).
The average percentage of commuters by residence municipality was around 33%, with lower values in the provincial capitals (only 7% in Milan), which, on the contrary, catalyzed entry flows (see Supplemental Material, Figure S1).
Estimated impacts by province of residence are reported in Table 1 as the posterior medians of AD^{A + B} and 80% CrIs. During 2007, we estimated that under the RS0 scenario there were 865.3 (80% CrI: 475.3, 1400.6) deaths attributable to annual levels of PM_{10} > 20 μg/m^{3} in the region as a whole, and that under the RS1 scenario, approximately 26% of these deaths (AD = 224.9; 80% CrI: 110.6, 367.8) could have been avoided by reducing the annual average of PM_{10} to the EU limit of 40 μg/m^{3}. The largest number of AD was estimated for the province of Milan [495.1 for the 20 μg/m^{3} limit (RS0) and 157.6 for the 40 μg/m^{3} limit (RS1)], followed by the provinces of Bergamo and Brescia. The smallest estimated impact was for the province of Sondrio, which is located in the mountain area, with only 3.1 deaths (80% CrI: 0.8, 5.9) under RS0, and a nearly null impact under RS1.
Table 1 – Attributable deaths (AD) among residents that were attributable to local and regional exposure (AD^{A + B}) under different scenarios by province: posterior medians and 80% CrIs.
We estimated that applying a 20% reduction of the annual average PM_{10} concentration in areas where the annual average was > 20 μg/m^{3}_{,} or a lower reduction if sufficient to reach the limit of 20 μg/m^{3} (RS2 scenario), would have prevented 311.4 (80% CrI: 167.3, 506.2) deaths in the region (Table 1). Applying a 20% reduction to reach up to 40 μg/m^{3} in areas where the annual average was > 40 μg/m^{3} (scenario RS3) would have prevented 189.4 (80% CrI: 98.1, 304.1) deaths, corresponding to 84% of the total estimated burden of mortality attributable to PM_{10} above the 40 μg/m^{3} EU limit under scenario RS1.
The province with the smallest estimated percentage of municipalities characterized by a nonnull impact was the province of Sondrio (59.6% of municipalities with a positive estimated value for AD_{i}^{A + B} under RS0, 6.3% under RS1), whereas the province of Milan had the largest percentage (98.7% and 97.8% under RS0 and RS1, respectively) (see Supplemental Material, Table S2).
In all provinces except Sondrio, we estimated ACRs of > 5 AD per 100,000 inhabitants due to PM_{10} > 20 μg/m^{3} (RS0), with a maximum of 12.7 in the province of Milano (Table 2). The estimated ACR for the entire Lombardy region was 9.1 per 100,000 inhabitants, but in the urban context of the capital cities the ACR reached 15.4/100,000. In fact, with few exceptions, ACRs were higher in provincial capitals than in the provinces as a whole. We estimated that PM_{10} > 40 μg/m^{3} (RS1) resulted in 2.4 deaths per 100,000 inhabitants in the Lombardy region, and 4.4 in the provincial capitals.
Table 2 – Attributable community rate (ACR) per 100,000 inhabitants under two scenarios by province of residence and in the capital of each province: posterior medians.
The larger cities belonging to the basin of the Po River were characterized by the largest estimated impacts in the region (Figure 2C). The posterior probability of a nonnull impact was large in the Milan and Brescia areas and in other urban centers of the Po valley (Figure 2D). The estimated impacts in Milano and Brescia stood out even when the less stringent limit of 40 μg/m^{3} was considered (see Supplemental Material, Figure S3C).
The probability that the number of exported AD (i.e., the deaths among nonresidents of i attributable to PM_{10} exposure in i; C_{i}) was greater than the number of imported AD (deaths among residents of i attributable to exposure outside of i; B_{i}), was > 60% for all capital cities, except for Lecco and Lodi, where it was around 50% (Table 3). However, the absolute differences between AD_{i}^{A + C} and AD_{i}^{A + B} were always < 1, except for Milan. We estimated 265.6 (80% CrI: 143.8, 414.9) deaths among residents of Milan attributable to PM_{10} exposure anywhere in the region, and 283.8 (80% CrI: 152.1, 443.7) deaths among residents of the region attributable to PM_{10} exposure in Milan, 22.3 (80% CrI: 12.1, 34.7) of which were among individuals residing outside of the capital (C_{i}).
Table 3 – Estimated values of AD^{A + B} and AD^{A + C} in the provincial capitals under the RS0 scenario: posterior medians, 80% CrIs, and probabilities that the “exported” attributable deaths exceed the “imported” ones.
Figure 3 shows the log ratio between the posterior median of C_{i} and the posterior median of B_{i}. The basin of the Po River was characterized by an overall balance between the two estimated flows, except for the cities of Milan, Brescia, Pavia, Cremona, and Mantova, where the ratio was large and positive, indicating that these cities tended to “export” impact elsewhere in the region rather than to “import” it. In contrast, our estimates suggest that the southern part of the Pavia province and the municipalities of the Alpine area, except Sondrio, tended mainly to “import” impact from the rest of the region.
Figure 3 – Log ratio between the posterior median of “exported” attributable deaths (C_{i}, attributable deaths among nonresidents due to exposure in city i) and the posterior median of the “imported” attributable deaths (B_{i}, attributable deaths among residents of city i due to exposure outside of city i), by municipality.
Discussion
Our analysis indicates that a 20% reduction in annual average PM_{10} concentrations > 40 μg/m^{3} would have substantially reduced the shortterm impact of PM_{10} exposures on population mortality in the Lombardy region in 2007. The largest estimated impacts were in the municipalities along the basin of the Po River and in the largest cities, particularly Milan and Brescia. Around 57% of the total estimated number of deaths attributable to annual average PM_{10} concentrations > 20 μg/m^{3} was in the province of Milan. For PM_{10} > 40 μg/m^{3}, this percentage rose to 70%. This is not surprising, because the Po basin is one of the less windy areas in Europe, with poor rainfall and frequent thermal inversions that trap pollutants close to the ground. In addition, the areas around Milano and Brescia are the most densely populated and the most industrialized, with tens of thousands of enterprises connected through road transportation.
In estimating the impact of PM_{10} exposures, we relaxed the strong hypothesis of a static population by considering commuting. This enabled us to identify areas where local PM_{10} exposure spread its impact throughout the region (Milan and its hinterland, and other main cities in the basin of the Po River) and municipalities whose residents were primarily impacted by exposure to PM_{10} in the other areas where they worked or attended school (the mountain municipalities). Not considering commuting would have resulted in underestimating the impact of pollution in those areas. Regarding the size of the phenomenon, the estimated number of deaths due to the exposure outside the municipality of residence was sometimes not negligible.
The MC approach allowed us to account for sampling variability resulting from the use of estimated values for PM_{10} effects, PM_{10} concentrations, mortality rates, and commuting probabilities, and for intrinsic variability related to random variation of deaths counts and commuting flows. Moreover, by sampling from the joint posterior distributions of the input quantities, we preserved betweenmunicipality correlations in the output so that reliable CrIs for the total number of AD in the region, and by province, could be derived by summing the estimated impacts over municipalities.
As a consequence of having considered many sources of uncertainty, CrIs were quite large. For example, the total estimated impact for the region ranged from 475.3 to 1400.6 AD with an 80% probability. These large CrIs indicate that one should focus on the whole AD distribution and not only on the point estimate. We reported also the posterior probabilities of AD > 0. The existence of a relevant impact of PM_{10} over the region is strongly supported by the fact that for the most municipalities the probability of a nonnull impact was > 90%; restricting the attention to the basin of the Po River, for most municipalities this probability was > 95%.
Individuals who commute from their residence area to work or attend school may be less frail and less susceptible to adverse effects of PM_{10} than the general population; if so, applying mortality rates and effect estimates derived for the general population to commuters could overestimate the commutingrelated impact of exposure. We think that, although some degree of overestimation is possible, it may be offset by underestimation of PM_{10} exposures among commuters, who are probably exposed to higher levels of air pollution than noncommuters because they travel along highly polluted “corridors,” and because PM_{10} concentrations during working hours tend to be higher than nighttime concentrations (Larssen et al. 2012). However, it is difficult to quantify these potential biases.
We addressed the issue of intermunicipality commuting, but did not consider intraurban commuting, which would require finer data at the suburban level and imply removing the assumption of homogeneous exposure within municipality. Also, although we did not account for commuting from neighboring regions, we accounted for commuting to them. This could explain unexpected low or high values of the log ratio between C_{i} and B_{i} for municipalities situated on the border of the region.
We did not consider “structural uncertainty” related to the model assumptions (Smith 2006). Some of these assumptions are common to analyses of the shortterm effect of air pollution (homogeneity of exposure within municipality, loglinear exposure–response relationship), and others are related to AD calculation in the presence of commuting. In particular, we assumed that the mortality rate depends on the municipality where the individual resides, whereas the PM_{10} effect depends on the municipality where the individual is exposed. The first assumption implies that the mortality rate depends on the socioeconomic, demographic, and environmental context where the individual resides. The second assumption implies that the exposure effect depends on local environmental factors that modulate the action of PM_{10} on health (for example, meteorological conditions, PM_{10} composition, or presence of special emission sources). This second assumption does not account for discrepancies among effect estimates that could be related to demographic and socioeconomic factors.
Finally we assumed that commuters spend a third of their time in the municipality where they work or study, supposing an overall balance between time spent traveling and weekend break. Changing this proportion might alter results at the municipality level.
Because our impact assessment refers to 1 year in the past, the epistemic component of the uncertainty related to the incomplete knowledge of future exposures and future sociodemographic conditions did not play a role in this analysis.
Conclusions
Using data on population daily commuting, we estimated for each municipality the shortterm impact of PM_{10} exposure not only on the mortality of residents, but also on the mortality of commuters residing elsewhere in the region, and we identified critical areas where PM_{10} pollution was likely to have the greatest impacts on population health. Our estimates accounted for several sources of uncertainty. Overall, our findings suggest an important impact of PM_{10} exposure on mortality, and also that different scenarios of PM_{10} reduction would have had a substantial beneficial effect, especially considering that we did not consider morbidity attributable to air pollution, which could be considerable. The results of this study can help policy makers prioritize interventions.
References
Anderson JO, Thundiyil JG, Stolbach A. 2012. Clearing the air: a review of the effects of particulate matter air pollution on human health. J Med Toxicol 8:166–175.
Baccini M, Biggeri A, Grillo P, Consonni D, Bertazzi PA. 2011. Health impact assessment of fine particle pollution at the regional level. Am J Epidemiol 174:1396–1405.
Ballester F, Medina S, Boldo E, Goodman P, Neuberger M, Iñiguez C, et al. 2008. Reducing ambient levels of fine particulates could substantially improve health: a mortality impact assessment for 26 European cities. J Epidemiol Community Health 62:98–105.
Berrocal VJ, Gelfand AE, Holland DM, Burke J, Miranda ML. 2011. On the use of a PM_{2.5} exposure simulator to explain birthweight. Environmetrics 22:553–571.
Besag J, Kooperberg C. 1995. On conditional and intrinsic autoregressions. Biometrika 82:733–746.
Besag J, York J, Mollié A. 1991. Bayesian image restoration, with two applications in spatial statistics. Ann Inst Stat Math 43:1–20.
Calder CA, Holloman CH, Bortnick SM, Strauss W, Morara M. 2008. Relating ambient particulate matter concentration levels to mortality using an exposure simulator. J Am Stat Assoc 103:137–148.
EU (European Union). 2008. Directive 2008/50/EC of the European Parliament and of the Council of 21 May 2008 on Ambient Air Quality and Cleaner Air for Europe. 2008. Offic J Eur Union 51(6):L152/1–L152/44. Available: http://eurlex.europa.eu/LexUriServ/LexUriServ.do?uri=OJ:L:2008:152:0001:0044:EN:PDF [accessed 1 October 2014].
Flachs EM, Sørensen J, Bønløkke J, BrønnumHansen H. 2013. Population dynamics and air pollution: the impact of demographics on health impact assessment of air pollution. J Environ Public Health 2013:760259; doi: 10.1155/2013/760259.
Heal MR, Heaviside C, Doherty RM, Vieno M, Stevenson DS, Vardoulakis S. 2013. Health burdens of surface ozone in the UK for a range of future scenarios. Environ Int 61:36–44.
International Agency for Research on Cancer. 2013. Air Pollution and Cancer. IARC Scientific Publication 161 (Straif K, Cohen A, Samet J, eds). Lyon:IARC.
Istat (Italian National Institute of Statistics). 2014. I.Stat: Your Direct Access to the Italian Statistics [in Italian]. Available: http://dati.istat.it [accessed 1 October 2014].
Knol AB, Petersen AC, van der Sluijs JP, Lebret E. 2009. Dealing with uncertainties in environmental burden of disease assessment. Environ Health 8:21; doi: 10.1186/1476069X821.
Kumamoto H, Henley EJ. 1996. Probabilistic Risk Assessment and Management for Engineers and Scientists. 2nd ed. New York, NY:IEEE Press.
Künzli N, Kaiser R, Medina S, Studnicka M, Chanel O, Filliger P, et al. 2000. Publichealth impact of outdoor and trafficrelated air pollution: a European assessment. Lancet 356(9232):795–801.
Larssen S, de Leeuw F, Kurfürst P. 2012. Estimating the contribution of commuting on exposure to particulate matter in European urban areas. ETC/ACC Technical Paper 2012/2. Available: http://acm.eionet.europa.eu/reports/docs/ETCACM_TP_2012_2_PM_exposure_urban_commuting.pdf [accessed 2 December 2014].
Li Y, Gibson JM, Jat P, Puggioni G, Hasan M, West JJ, et al. 2010. Burden of disease attributed to anthropogenic air pollution in the United Arab Emirates: estimates based on observed air quality data. Sci Total Environ 408:5784–5793.
Lunn DJ, Thomas A, Best N, Spiegelhalter D. 2000. WinBUGS—a Bayesian modelling framework: concepts, structure, and extensibility. Stat Comput 10:325–337.
Martuzzi M, Mitis F, Iavarone I, Serinelli M. 2006. Health Impact of PM_{10} and Ozone in 13 Italian Cities. Copenhagen:World Health Organization Regional Office for Europe. Available: http://www.euro.who.int/__data/assets/pdf_file/0012/91110/E88700.pdf [accessed 1 October 2014].
Medina S, Plasencia A, Ballester F, Mücke HG, Schwartz J, Apheis Group. 2004. Apheis: public health impact of PM_{10} in 19 European cities. J Epidemiol Community Health 58:831–836.
MesaFrias M, Chalabi Z, Vanni T, Foss AM. 2013. Uncertainty in environmental health impact assessment: quantitative methods and perspectives. Int J Environ Health Res 23:16–30.
Orru H, Teinemaa E, Lai T, Tamm T, Kaasik M, Kimmel V, et al. 2009. Health impact assessment of particulate pollution in Tallinn using fine spatial resolution and modeling techniques. Environ Health 8:7; doi: 10.1186/1476069X87.
Regione Lombardia. 2004. Pendolari in Lombardia. Censimento 2011 [in Italian]. Available: http://www.statistica.regione.lombardia.it/DocumentiNews/Pendolari_Lomb_cens_2001.pdf [accessed 1 October 2014].
Riley RD, Higgins JPT, Deeks JJ. 2011. Interpretation of random effects metaanalyses. BMJ 342:d549; doi: 10.1136/bmj.d549.
Silibello C, Calori G, Brusasca G, Giudici A, Angelino E, Fossati G, et al. 2008. Modelling of PM_{10} concentrations over Milano urban area using two aerosol modules. Environ Model Softw 23:333–343.
Smith EP. 2006. Uncertainty Analysis. In: Encyclopedia of Environmetrics. 4 (ElShaarawi AH, Piegorsch WW, eds). Hoboken, NJ:John Wiley & Sons Ltd.; doi:10.1002/9780470057339.vau001.pub2.
Wacholder S. 2005. The impact of a prevention effort in a community. Epidemiology 16:1–3.
WHO (World Health Organization). 2005. Air Quality Guidelines: Global Update 2005. Particulate Matter, Ozone, Nitrogen Dioxide and Sulfur Dioxide. Copenhagen:WHO. Available: http://www.euro.who.int/__data/assets/pdf_file/0005/78638/E90038.pdf [accessed 1 October 2014].
WHO (World Health Organization). 2013. Review of Evidence on Health Aspects of Air Pollution—REVIHAAP Project. Technical Report. Copenhagen:WHO Regional Office for Europe. Available: http://www.euro.who.int/__data/assets/pdf_file/0004/193108/REVIHAAPFinaltechnicalreportfinalversion.pdf [accessed 1 October 2014].
Announcements
EHP has joined the bioRxiv B2J program! Effective immediately, authors who have already posted their paper to bioRxiv can choose to submit their work to EHP through the bioRxiv platform. (more…)