Skip to content

Environmental Health Perspectives

Facebook Page EHP Twitter Feed Open Access icon  

Research Volume 123 | 2015

Environ Health Perspect; DOI:10.1289/ehp.1408218

Commuting-Adjusted Short-Term 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 Biggeri1,2

Author Affiliations open
1Department of Statistics, Informatics and Applications “G. Parenti,” University of Florence, Florence, Italy; 2Biostatistics Unit, Cancer Prevention and Research Institute (ISPO), Florence, Italy; 3Epidemiology Unit, Department of Preventive Medicine, Fondazione IRCCS (Istituto Di Ricovero e Cura a Carattere Scientifico) Ca` Granda, Ospedale Maggiore Policlinico, Milan, Italy; 4Department of Clinical Sciences and Community Health, Università degli Studi di Milano, Milan, Italy

PDF icon PDF Version (3 MB)

  • Background: Exposure to air pollution is associated with a short-term increase in mortality, and this field has begun to focus on health impact assessment.

    Objectives: Our aim was to estimate the impact of PM10 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, PM10 concentrations, estimated PM10 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 PM10 concentrations above the World Health Organization-recommended threshold of 20 μg/m3 were responsible for 865 short-term deaths (80% credibility interval: 475, 1,401), 26% of which were attributable to PM10 above the European Union limit of 40 μg/m3. Reducing annual average PM10 concentrations > 20 μg/m3 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 short-term impact of PM10 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. Commuting-adjusted short-term health impact assessment of airborne fine particles with uncertainty quantification via Monte Carlo simulation. Environ Health Perspect 123:27–33;

    Address correspondence to M. Baccini, Department of Statistics, Informatics and Applications “G. Parenti,” viale Morgagni 59-50134 Firenze, Italy. Telephone: 39 320 4209085. E-mail:

    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 Scientific 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

  • PDF icon Supplemental Material (435 KB) PDF


The role of air pollution in short-term and long-term 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 short-term 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 long-term 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, population-level 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 short-term 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. Mesa-Frias 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 higher-polluted areas (Larssen et al. 2012).

Our aim in this paper was to estimate the short-term impact of high concentrations of particles up to 10 μm in diameter (PM10) 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 PM10 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,865-km2 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 - Map of study area.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.

View larger image (TIF File)

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 PM10 in Milan, the capital of the region (1,299,633 inhabitants in 2007), was 52.5 μg/m3, with 95% of daily concentrations > 20 μg/m3, whereas the average across the other most densely populated areas of Lombardy was 45.4 μg/m3 (Baccini et al. 2011).


Data and input distributions. To derive the probabilistic distributions for input data used in the MC procedure, we used information on total mortality, PM10 concentrations, PM10 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 deathsi be the total number of deaths during the period 2000–2004 in municipality i (i = 1, 2,…,1,546). We assumed that deathsi followed a Poisson distribution with mean given by the product of the death rate ri and the population size Di, which we estimated to be five times the population size in 2001. Moreover, we specified a random-effect log-linear model on ri:

log(ri) = ui + vi, [1]

where ui and vi were independent random terms representing the unstructured spatial variability component and the structured spatial variability component, respectively. We assumed that ui were independent and normally distributed with mean 0 and variance τu2, and that vi followed an intrinsic conditional autoregressive (ICAR) model. In other words, for each Si, the set of the ni municipalities adjacent to municipality i, we assumed the following conditional distribution for vi:

vi | vj Є Si ~ N(v̄i, τv2/ni), [2]

where v̄i was the mean of vi for the municipalities belonging to Si, and τv2/ni was their conditional variance (Besag and Kooperberg 1995). Through the random terms ui and vi, 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 ri. 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 PM10. Data about annual concentrations of PM10 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 = σ2f(dij;κ,φ), where f(dij;κ,φ) = exp(–φ dij κ), σ2 was a variance parameter, dij 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 PM10 at the cell centroids of the entire grid. Finally, we derived the joint posterior predictive distribution for the annual average concentration of PM10 for each municipality by integrating over these predicted cell values.

We confirmed the validity of the geostatistical model through leave-one-out cross-validation. 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 PM10 effect. Estimates of PM10 effect were derived from the Bayesian random effects meta-analysis in Baccini et al. (2011), which combined estimates of the effect of PM10 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 city-specific effects for the 13 areas included in the meta-analysis, and the posterior predictive distribution of a generic city-specific 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 (cij, ji). Specifying a vague a priori beta distribution on the commuting probability from municipality i to j (pij), and a binomial likelihood for cij with binomial denominator equal to the population of i at the 2001 census (popi2001), the resulting posterior distribution for pij was a beta distribution with parameters 1 + cij and 1 + popi2001cij. 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 pij (where ji) 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 PM10 depends only on the exposure–response curve and the concentration of PM10 in municipality i. Then, if one assumes linearity on a log scale, the deaths attributable to exposures exceeding the threshold x0 (AD i) can be calculated for each municipality as

ADi = yiyi/exp[β i × (xix0) × I(xi > x0)], [3]

where yi is the total number of deaths in municipality i, β i is the estimated PM10 effect in i, xi is the annual average PM10 concentration in i, and I(xi > x0) is a indicator function with I = 1 for xi > x0 and I = 0 otherwise. Usually yi is assumed to follow a Poisson distribution with mean μi = popi × ri, where ri is the crude mortality rate and popi represents the person-years 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 PM10 during the day, so that Equation 3 is no longer valid. Therefore, for each municipality i we calculated three different quantities: Ai, the number of deaths among residents of i attributable to exposure in i; Bi, the number of deaths among residents of i attributable to exposure in j (where ji); and Ci, the number of deaths among non-residents of i attributable to exposure in i.

Assuming that, on average, regular commuters spend one-third of their time in municipality j (where they work or study), and two-thirds of their time in municipality i (where they live),

Ai = yiSyiS/exp[βi × (xi – x0) × I(xi > x0)], [4]

where yiS is the total number of deaths among the residents in municipality i associated with the person-years 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 yiS followed a Poisson distribution, with mean μiS equal to the product between person-years at risk and the crude mortality rate ri:

μiS = [popi – 1/3 × Σj i exit(i,j)] × ri, [5]

where popi 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 ji). 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 popi. Because the probabilities of commuting were small relative to the probability of not commuting, no constraint was needed to avoid negative values of μiS.

Bi is the sum over j, with ji, of deaths among residents of i attributable to the exposure in j, which were estimated based on the person-years at risk spent in municipality j by residents of i, the crude mortality rate for municipality i (ri), and the attributable risk associated with exposure in j:

Bi = Σ j i {yijCyijC/exp[βj × (xj – x0) × I(xj > x0)]}, [6]

where yijC 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 μijC = 1/3 × exit(i,j) × ri.

Finally, Ci is the sum over j, with ji, of deaths among commuters from j to i attributable to exposure in i, which were estimated based on the person-years at risk spent in municipality i by residents of j, the crude mortality rate for municipality j (rj), and the attributable risk associated with exposure in i:

Ci = Σj i {yjiCyjiC/exp[βi × (xix0) × I(xi > x0)]}, [7]

where yjiC, the total number of deaths among commuters from j to i, was assumed to follow a Poisson distribution with mean μjiC = 1/3 × exit(j,i) × rj. 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 pji and the number of trials equal to popj.

As the formulas show, we assumed that the PM10 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 Ai, Bi, and Ci, we estimated two impact measures:

  • ADiA +B = Ai + Bi, 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
  • ADiA +C = Ai + Ci, 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 ΣiBi = ΣiCi, which is the total estimated number of deaths due to commuting-related exposure in the region, and ΣiADiA + B = ΣiADiA + C.

Air pollutant reduction scenarios. We estimated AD under different reduction scenarios (RS) corresponding to different definitions of the threshold x0:

  • RS0: x0 = 20 μg/m3, the WHO Air Quality Guideline threshold for PM10 annual average (WHO 2005);
  • RS1: x0 = 40 μg/m3, the European Union (EU) limit for PM10 annual average (EU 2008);
  • RS2: x0 equal to a reduction of 20% in the observed annual concentration of PM10, provided it is > 20 μg/m3;
  • RS3: x0 equal to a reduction of 20% in the observed annual concentration of PM10, provided it is > 40 μg/m3.

The RS2 and RS3 scenarios define intermediate targets coherent with realistic policies of progressive reduction in PM10 concentration until the limits of 40 μg/m3 or 20 μg/m3 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 (Ai, Bi, Ci, ADiA + B, ADiA + C; i = 1,2,…,1,546).

MC simulation was performed using R statistical software (R Core Team; The code is available at​bility/

Evaluating the role of commuting. For each municipality we calculated the logarithm of the ratio between the posterior median of Ci and the posterior median of Bi, 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 commuting-related 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 PM10 effect was sampled during the MC procedure, because we were interested in the AD distribution conditional on rejection of the null hypothesis of no PM10 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.


The highest annual average concentrations of PM10 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/m3. Milan and the neighboring municipalities were characterized by annual average concentrations > 50 μg/m3.

Figure 2 - Four maps.Figure 2 – Posterior means of annual average PM10 concentrations (μg/m3) (A) and smoothed crude mortality rates (per 100,000 residents) (B); posterior medians of attributable deaths among residents (C) and posterior probabilities (%) of a non-null impact (D) under the RS0 scenario (for annual average PM10 > 20 μg/m3), by municipality.

View larger image (TIF File)

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 short-term effect of PM10 were similar among cities except for Milan, where the estimated effect was twice the overall meta-analytic 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 ADA + 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 PM10 > 20 μg/m3 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 PM10 to the EU limit of 40 μg/m3. The largest number of AD was estimated for the province of Milan [495.1 for the 20 μg/m3 limit (RS0) and 157.6 for the 40 μg/m3 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 - See HTML for full tableTable 1 – Attributable deaths (AD) among residents that were attributable to local and regional exposure (ADA + B) under different scenarios by province: posterior medians and 80% CrIs.

View Table (HTML Version)
View larger image (TIF File)

We estimated that applying a 20% reduction of the annual average PM10 concentration in areas where the annual average was > 20 μg/m3, or a lower reduction if sufficient to reach the limit of 20 μg/m3 (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/m3 in areas where the annual average was > 40 μg/m3 (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 PM10 above the 40 μg/m3 EU limit under scenario RS1.

The province with the smallest estimated percentage of municipalities characterized by a non-null impact was the province of Sondrio (59.6% of municipalities with a positive estimated value for ADiA + 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 PM10 > 20 μg/m3 (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 PM10 > 40 μg/m3 (RS1) resulted in 2.4 deaths per 100,000 inhabitants in the Lombardy region, and 4.4 in the provincial capitals.

Table 2 - See HTML for full tableTable 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.

View Table (HTML Version)
View larger image (TIF File)

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 non-null 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/m3 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 PM10 exposure in i; Ci) was greater than the number of imported AD (deaths among residents of i attributable to exposure outside of i; Bi), was > 60% for all capital cities, except for Lecco and Lodi, where it was around 50% (Table 3). However, the absolute differences between ADiA + C and ADiA + B were always < 1, except for Milan. We estimated 265.6 (80% CrI: 143.8, 414.9) deaths among residents of Milan attributable to PM10 exposure anywhere in the region, and 283.8 (80% CrI: 152.1, 443.7) deaths among residents of the region attributable to PM10 exposure in Milan, 22.3 (80% CrI: 12.1, 34.7) of which were among individuals residing outside of the capital (Ci).

Table 3 - See HTML for full tableTable 3 – Estimated values of ADA + B and ADA + C in the provincial capitals under the RS0 scenario: posterior medians, 80% CrIs, and probabilities that the “exported” attributable deaths exceed the “imported” ones.

View Table (HTML Version)
View larger image (TIF File)

Figure 3 shows the log ratio between the posterior median of Ci and the posterior median of Bi. 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 - Map.Figure 3 – Log ratio between the posterior median of “exported” attributable deaths (Ci, attributable deaths among nonresidents due to exposure in city i) and the posterior median of the “imported” attributable deaths (Bi, attributable deaths among residents of city i due to exposure outside of city i), by municipality.

View larger image (TIF File)


Our analysis indicates that a 20% reduction in annual average PM10 concentrations > 40 μg/m3 would have substantially reduced the short-term impact of PM10 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 PM10 concentrations > 20 μg/m3 was in the province of Milan. For PM10 > 40 μg/m3, 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 PM10 exposures, we relaxed the strong hypothesis of a static population by considering commuting. This enabled us to identify areas where local PM10 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 PM10 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 PM10 effects, PM10 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 between-municipality 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 PM10 over the region is strongly supported by the fact that for the most municipalities the probability of a non-null 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 PM10 than the general population; if so, applying mortality rates and effect estimates derived for the general population to commuters could overestimate the commuting-related impact of exposure. We think that, although some degree of overestimation is possible, it may be offset by underestimation of PM10 exposures among commuters, who are probably exposed to higher levels of air pollution than noncommuters because they travel along highly polluted “corridors,” and because PM10 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 Ci and Bi 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 short-term effect of air pollution (homogeneity of exposure within municipality, log-linear 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 PM10 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 PM10 on health (for example, meteorological conditions, PM10 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.


Using data on population daily commuting, we estimated for each municipality the short-term impact of PM10 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 PM10 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 PM10 exposure on mortality, and also that different scenarios of PM10 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.


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 PM2.5 exposure simulator to explain birthweight. Environmetrics 22:553–571.

Besag J, Kooperberg C. 1995. On conditional and intrinsic auto-regressions. 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:​​:PDF [accessed 1 October 2014].

Flachs EM, Sørensen J, Bønløkke J, Brønnum-Hansen 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: [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/1476-069X-8-21.

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. Public-health impact of outdoor and traffic-related 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:​/ETCACM_TP_2012_2_PM_exposure_urban_comm​uting.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 PM10 and Ozone in 13 Italian Cities. Copenhagen:World Health Organization Regional Office for Europe. Available:​f_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 PM10 in 19 European cities. J Epidemiol Community Health 58:831–836.

Mesa-Frias 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/1476-069X-8-7.

Regione Lombardia. 2004. Pendolari in Lombardia. Censimento 2011 [in Italian]. Available: http://www.statistica.regione.lombardia.​it/DocumentiNews/Pendolari_Lomb_cens_200​1.pdf [accessed 1 October 2014].

Riley RD, Higgins JPT, Deeks JJ. 2011. Interpretation of random effects meta-analyses. 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 PM10 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 (El-Shaarawi 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:​f_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:​f_file/0004/193108/REVIHAAP-Final-techni​cal-report-final-version.pdf [accessed 1 October 2014].

WP-Backgrounds Lite by InoPlugs Web Design and Juwelier Schönmann 1010 Wien