This article appeared in a journal published by Elsevier. The attached

copy is furnished to the author for internal non-commercial research

and education use, including for instruction at the authors institution

and sharing with colleagues.

Other uses, including reproduction and distribution, or selling or

licensing copies, or posting to personal, institutional or third party

websites are prohibited.

In most cases authors are permitted to post their version of the

article (e.g. in Word or Tex form) to their personal website or

institutional repository. Authors requiring further information

regarding Elsevier’s archiving and manuscript policies are

encouraged to visit:

http://www.elsevier.com/copyright
Author's personal copy

Cumulative effects and threshold levels in air pollution mortality: Data analysis

of nine large US cities using the NMMAPS dataset

Mario Stylianou a,1, Mark J. Nicolich b,*

a 9517 Lindale Drive, Bethesda, MD 20817, USA

b Epidemiology and Health Surveillance, ExxonMobil Biomedical Sciences, Inc., 1545 Route 22 East, P.O. Box 971, Annandale, NJ 08801-0971, USA

Exploratory analyses suggest that a linear, no threshold, model for pollution studies is not consistent with the observed data.

article info

Article history:

Received 21 June 2008

Received in revised form

2 April 2009

Accepted 9 April 2009

Keywords:

Air pollution

Cumulative effects

Thresholds

RR

abstract

We examined the existence of thresholds, cumulative effects and the homogeneity of five air pollutants

on the relative risk of three mortality outcomes using data from nine major US cities using data from

NMMAPS. Overall, PM10 (usually 200-day accumulation) and ozone (3-day accumulation) were the two

important predictors of outcome but their effect was not uniform across the nine cities. Many models

exhibited thresholds (25–45 mm g/m3 for PM10, and 10–45 ppb for O3). Our preliminary exploratory

analyses suggest that the use of a linear, no threshold, model for pollution studies is not consistent with

the observed data. The heterogeneity in the risk estimates across the nine cities suggests combining the

local risk estimates to obtain a national risk estimate may not be justifiable and the estimate is likely to

be confounded.

2009 Elsevier Ltd. All rights reserved.

1. Introduction

Many authors have developed statistical models to investigate

the association between air pollutant levels, such as ozone (O3),

particulate matter of less than 10 mm or 2.5 mm diameter (PM10 or

PM2.5), nitrogen dioxide (NO2), sulphur dioxide (SO2), and carbon

monoxide (CO) with morbidity or mortality outcomes. Much of the

statistical research has focused on time series related methods using

lag times of some short length to evaluate the association between

a single pollutant and the outcome of interest, adjusting for weather,

age, day of the week and time trends (Schwartz, 2000; Bell et al.,

2004, 2006; Roberts and Martin, 2006; among many others).

Although most of the literature deals with the short-term effects

of air pollutants on mortality, it seems reasonable to assume that

pollutants may increase the risk for adverse events over longer

duration of exposure because a narrow moving window may not

adequately capture larger effect variations of pollutants under

study. Several researchers have considered longer-term effects;

Schwartz (2000) studied the effect of PM2.5 and weather variables

using moving averages of up to 60 days. Janes et al. (2007)

considered mean monthly exposures over the previous 12 months

with an approximate 30-day window. More recently, Schwartz

et al. (2008) explored the effect of PM2.5 exposure to mortality for

up to the previous five years. There is no consensus in the literature

regarding the lag or cumulative time needed to evaluate the

pollutants effect.

Another consideration in air pollution models is the combined

effect of several pollutants. Bell et al. (2004) estimated same day

non-accidental mortality using cumulative ozone exposure during

the previous week considering PM10 and the usual confounders of

weather, day of the week, and calendar time. They report that ozone

is associated with daily non-accidental mortality (and the subsets of

cardiovascular and respiratory mortality), and these results are

independent of adjustment for PM10, weather, seasonality and longterm trends. However, PM10 was considered only with a 1-day lag

which is different from most published models and may explain why

PM10 had no effect in the model with ozone. Bell et al. (2005)

replicated their results in a subsequent meta-analysis of ozone

mortality studies, again with short lag times for exposures to ozone

and PM10 or PM2.5. It is possible that a PM effect on the death count

may be more pronounced when it is cumulated over longer times.

Many of the models assessing the effect of several pollutants on

death rates produce relatively poor fitting or biased estimates of

exposure–response association (Zeka and Schwartz, 2004). One

explanation is that these pollutants are usually correlated in the

* Corresponding author. Current address: 24 Lakeview Rd, Lambertville,

NJ 08530, USA. Tel.: þ1 609 397 4089.

E-mail addresses: stylianm@gmail.com (M. Stylianou), mark.nicolich@gmail.com

(M.J. Nicolich).

1 Tel.: þ1 301 435 0441; fax: þ1 301 480 1862.

Contents lists available at ScienceDirect

Environmental Pollution

journal homepage: www.elsevier.com/locate/envpol

0269-7491/$ – see front matter 2009 Elsevier Ltd. All rights reserved.

doi:10.1016/j.envpol.2009.04.011

Environmental Pollution 157 (2009) 2216–2223
Author's personal copy

short-term and the correlation leads to unstable estimates of the

pollution effects. One approach to the problem of correlated shortterm data is to use long-term average or cumulative pollution data

so the effect of correlation among long-term cumulative pollutants

in the modeling process would be minimized. This is especially true

when these correlations are relatively small, such as with PM10 and

ozone. In general, the correlation with longer-term cumulated

series is not a large problem compared to the problems of nonlinearity or model over-fitting.

The work of Xia and Tong (2006) suggests a method to account

for the cumulative effects of longer-term exposures to air pollution

by using a weighted cumulative sum for each pollutant without the

restriction of a linear dose–response function. The weights are

derived by a restricted minimization method where the weights are

non-increasing (or decreasing) the further we go into the past.

Their analysis of data from Hong Kong, El Paso, Chicago, and

Pittsburgh suggests that cardiovascular disease (CVD) mortality is

affected by a cumulative effect of PM10 for about 200 days, with

roughly uniform weights, and that thresholds exist for all pollutants, noticeably for PM10 and ozone. Their study suggests that

different pollutants may require different cumulative periods.

Bell et al. (2006) have investigated the exposure–response

relationship for ozone and risk of total mortality and explored the

existence of thresholds using linear and non-linear approaches.

They concluded that if an ozone threshold exists it must be at low

levels (15–19 ppb average daily levels). Although they used a twostage approach to model each community separately and then

combine information across communities, they only reported the

combined national average risk. However, Janes et al. (2007) have

reported that a national risk average may be confounded. Others

have also investigated the existence of thresholds in PM and ozone

with mixed results (Daniels et al., 2000; Koop and Tole, 2006;

Leitenstorfer and Tutz, 2007; Schwartz et al., 2008).

The objective of this exploratory analysis is to investigate the

simultaneous cumulative effects of O3 and PM on non-accidental

mortality and the existence of thresholds for these pollutants. We

explore the assumption of monotonicity between the risk function

and pollutant levels. We think that a monotonic (but not necessarily linear) assumption for a pollutant effect is reasonable in the

absence of plausible mechanisms for beneficial effects and their

likelihood that these effects at low concentrations are possibly due

to chance. In addition, we explore the cumulative effect of other

pollutants, such as NO2, SO2, and CO on death outcomes. Finally

we question the usefulness of combining the city-specific risk

estimates to form a national average. This investigation is based

on daily time-series data from nine major US cities; Baltimore,

Chicago, Dallas/Fort Worth, Los Angeles, Miami, New York,

Philadelphia, Pittsburgh, and Seattle, which have a combined

population of more than 35 million people.

2. Methods

2.1. Data collection

The data are extracted from the updated National Morbidity, Mortality, and Air

Pollution Study (NMMAPS) for 1987–2000.

The response variables considered for these analyses are counts of:

a) total non-injury related mortality

b) cardiovascular disease (CVD) mortality

c) respiratory disease mortality

Because we are doing an exploratory data analysis, we did not use all the cities in

the NMMAPS dataset. As has been done by many researchers analyzing NMMAPS

data (Ren et al., 2008; Samet et al., 2000; Schwartz, 1999; Dockery et al., 1993), we

selected a subset of the dataset that included major cities with large populations,

geographically dispersed throughout the US, and having generally complete data

records. The nine cities were selected before the statistical analysis began and no

cities were dropped after selection. The analyses of Bell et al. (2004), based on 95 US

cities, suggest that the effect of ozone on mortality for the youngest and oldest age

groups were similar, and therefore, we think it is reasonable to assume that age is

not a major factor in estimating the RRs, so we have combined the three NMMAPS

age groups in order to simplify our presentation. A major reason for choosing

a subset of the available data rather than the entire set was to alleviate the

computational requirements imposed by some analyses described in the statistical

section. The non-pollutant explanatory variables considered include temperature,

dewpoint, day of the week (day), and calendar time trend (time)

2.2. Statistical approach

Our methods are motivated by the newly developed approach of Xia and Tong

(2006) that considers weighted cumulative sums of exposure. These authors have

estimated the effect of various air pollutants over time for Chicago, Pittsburgh, El Paso,

and Hong Kong among other cities. Their results show long-term effect (about 200

days) of PM10 on CVD in Chicago and Pittsburgh, and on respiratory diseases in Hong

Kong. Their results suggest that the ozone effect is short term. Given that these results

were consistent among the different cities they studied, we have implemented those

assumptions into our analyses. However, to reduce the heavy computational burden

of deriving exact daily weights for the cumulative effects, we use their results and

linear interpolation to approximate the weights. This process resulted in uniform

weights for the last 200 days for PM10 and linearly decreasing weights for O3, NO2,

SO2, and CO for the last three days (current plus previous two days). Also, linearly

decreasing weights were used for temperature and dewpoint, for the last 30 days.

Uniform weights suggest that all days considered in the cumulative estimate have

equal weights. In contrast, the linearly decreasing weights assume that the further

into the past we look the smaller the assigned weights are. Beyond these cutoff points,

the weight is 0. Therefore, the number of consecutive days with a non-zero weight

defines the lag time for each variable. We calculated the cumulative weighed total

exposure (or a weighted mean) for each explanatory variable. We also investigated the

robustness of the above weights using different cumulative periods.

For each cause of death and for each city, we first fit a model that considers

calendar variables (day, time) and meteorological variables (temperature, dewpoint). Only those variables that are significantly associated with the outcome

(at p 0.05) remained in the model to form the ‘‘basic model’’. This basic model may

be different between cities and among different causes of death.

We added each pollutant individually to assess its importance in predicting the

number of events. Each pollutant associated with the response in the basic model

would be further considered for a multivariate model.

We modeled the effects of pollutants on each of the daily causes of death using

the generalized additive model (GAM) as implemented in the R library mgcv.

Specifically, we fit a Poisson model with a log-link. For example, a ‘‘full’’ model in

a multivariate setting with a basic model plus all five pollutants could be:

EðlogðYtÞÞ 1⁄4 m þ X6

k 1⁄4 1

bkDk;t þ g1

 X29

s 1⁄4 0

q1;sTts

!

þ g2

 X29

s 1⁄4 0

q2;sDPts

!

þ g3ðTimetÞ þ g4

 X199

s 1⁄4 0

q4;sPMts

!

þ g5

 X2

s 1⁄4 0

q5;sOts

!

þ g6

 X2

s 1⁄4 0

q6;sNO2ts

!

þ g7

 X2

s 1⁄4 0

q7;sSO2ts

!

þ g8

 X2

s 1⁄4 0

q8;sCOts

!

ð1Þ

where m is the intercept term, the subscript t represents study time, the summations

are over the cumulative period of interest (s), Yt is the daily number of deaths at time

t, Dk,t is the indicator function for day of the week, bk is the mortality adjustment for

the kth day of the week, Tt is the average daily temperature, DPt is the average daily

dew point, and Timet is calendar time. PMt and Ot are the PM10 and ozone levels at

time t, while, NO2t, SO2t and COt are the corresponding pollutants at time t. g1(x), .,

g8(x) are unknown smooth functions, qk;s is the weight for the kth term at s days

prior to the current day, PL

s 1⁄4 0 qk;s 1⁄4 1 and 0 qk;L qk;L1 / qk;0. The number

of knots used for temperature, dew point, and time are 6, 6, and 7/year respectively,

which are the same as used by many NMMAPS investigators (Daniels et al., 2000;

Samet et al., 2000).

We use p-values not as a strict decision value, but in the sense of Weinberg

(2001) and Poole (2001) to provide a useful scale to evaluate the relative usefulness

of terms in the models of different cities.

Model (1) is flexible, so one can incorporate both constrained (monotone) and

unconstrained (no monotonic assumption) smooth terms. The unconstrained

smooth terms are estimated using the gam function as implemented in the R library

mgcv with penalized regression splines (Eilers and Marx, 1996), and a manual

backward elimination procedure was used to remove factors from the full model

that do not contribute to the prediction of the outcome based on p-values, AIC

change, and other residual analyses.

For the constrained model, we expected to use a similar procedure as in

Leitenstorfer and Tutz (2007) employing B-splines and monotone regression.

However, the approach involves inverting an N  N matrix where N is the number of

daily measurements (w5000), but computational restrictions limited the number of

measurements we could use. One possible solution was to take a subset of the dataset

M. Stylianou, M.J. Nicolich / Environmental Pollution 157 (2009) 2216–2223 2217
Author's personal copy

(such as the most recent or oldest) although doing so might limit the variability of

each pollutant, resulting in a less parsimonious model. Also, a time trend might be lost

when using a shorter period. We chose to take a random sample of the entire dataset

to cover the whole time-span of NMMAPS, thereby incorporating all possible time

trends. Because we needed to have consecutive days to obtain the cumulative effects

of each factor, we preprocessed the dataset before a random sample was selected.

To preprocess the dataset, we first imputed any missing pollutant value with the

average of the two adjacent observations. Then, for each day, we derived the

weighted cumulative levels for each factor listed in equation (1) except Time. This

method incorporates the lag time for each variable, since the lag is determined by

the number of consecutive days with non-zero weights. These daily cumulative

values then became our ‘‘new’’ predictors. Approximately 15–20% of the dataset was

used in the selected random samples that spanned the 14-year-period.

The model fit was assessed both quantitatively and graphically using residual

and partial residual plots. Autocorrelation and partial autocorrelation plots were

used to assess the existence of any autocorrelation remaining in the residuals.

Estimates of the RR were derived using the final model. The RR at exposure x is

the risk from being exposed at level x relative to the risk at the base level. When the

effect-function is linear any base can be used because the effect change is the same

at all points along the line. Risk can then be reported as (Bell et al., 2004) (p. 2375) do

‘‘a 0.15% percent increase in daily mortality for a 10-ppb increase in ozone .’’ If the

form of the effect-function is non-linear, the risk must be reported relative to

a selected base. Some authors choose the first quartile of exposure as the base level

(Miller et al., 2007). We have chosen to report risk relative to the lowest observed

exposure as do Leitenstorfer and Tutz (2007):

RRðtÞ 1⁄4 EðeventjPðtÞÞ

EðeventjPðminÞÞ 1⁄4 expðgPðtÞÞexpðgðPðminÞÞÞ (2)

where P(t) is the pollutant level at time t, P(min) is the minimum pollutant level

observed, and g is a smooth function.

A visual inspection of the RR curves and effect-functions by pollutant levels may

indicate the existence of thresholds. A threshold is suggested when a ‘‘hockey stick’’

shape is observed. One could also use the 95 percent confidence band for RR to

determine the existence of a threshold.

Lastly, in order to compare and contrast the individual city models with a model

that is the same across cities, we have used a unified model that incorporated the

same predictors for all cities. These predictors were time and meteorological variables and pollutants that are associated with death events across many cities.

3. Results

We present the model form and estimated RR for the basic

model by city. Then we present results for the unified model when

the same variables are used in all city-based models. Lastly, we

present the results by death outcome.

3.1. Basic models

The basic models and which pollutants remain in the final

models differ by city (Table 1). Day, time, and the cumulative effects

of weather variables are the most frequent predictors of death in

most cities. If the cumulative effect of PM10 was not associated over

a period of 200 days, we then examined shorter durations of 90, 30,

or 7 days. For example, the 200-day cumulative effect of PM10 was

not an important predictor of CVD deaths for NY and Seattle, but the

90-day cumulative effect of PM10 was. Among the five pollutants

studied, ozone and PM10 were the most frequent important

predictors of death. However these results are not consistent across

different cities. For example, neither ozone nor PM10 was associated

with any type of deaths in Miami. In addition, PM10 was not associated with any type of death in Baltimore, Dallas/Fort Worth, and

Philadelphia. Neither is ozone a consistent predictor, although it

significantly predicts at least one type of death for most cities. CO

was neither a consistent nor a major predictor of death in the cities

studied. All the final models were free of residual autocorrelation or

any residual temporal effects.

Of the three types of deaths assessed (total, CVD, and respiratory), respiratory deaths have the weakest pollutant associations

and total deaths outcome has the strongest associations.

Using the models of Table 1, we estimated for each city the RR and

threshold levels of predictors by type of death (Table 2) under the

unconstrained and monotone assumptions. We also have evaluated

the random samples used in the analyses under the monotonic

assumptions and found them to be representative of the full data.

Table 2 indicates that 3-day cumulative ozone is the most likely

pollutant to have a threshold with an average value of 32 ppb using

the non-restricted models for the six cities that appear to have

a threshold with a total death outcome. Using a monotonic

assumption, we obtain an average threshold level at 21 ppb. Using

CVD death as an outcome, the average threshold with a nonrestricted and monotonic models is at 28 and 31 ppb. With the 200-

day cumulative PM10, a threshold is suggested only for Chicago

using the unconstrained model 29 and 30 mg/m3 for total and CVD

deaths. Under monotonic model the average threshold is 31 and

30 mg/m3 for total and CVD deaths. Figures of risk functions for the

non-restricted and monotonic models are given in the online

supplement. The monotonic risk estimates are always 1.

3.2. Unified model analyses

The basic model analyses suggest that PM10 and ozone are the

most likely pollutants predictive of different accountings of

mortality. We evaluated the consistency of the pollutant effect

across the different cities using a model that adjusts for the same

variables, that is, the basic model plus PM10 and ozone;

EðlogðYtÞÞ 1⁄4 m þ X

6

k 1⁄4 1

bkDk;t þ g1

 X

29

s 1⁄4 0

q1;sTts

!

þ g2

 X

29

s 1⁄4 0

q2;sDPts

!

þ g3ðTimetÞ

þ g4

 X

199

s 1⁄4 0

q4;sPMts

!

þ g5

 X

2

s 1⁄4 0

q5;sOts

!

ð3Þ

Table 1

Final pollution models for each of the nine US cities. Each city model may be

different.

Event City Terms in the Final Model

dow dew Time Temp PM10 O3 NO2 SO2 CO

Total Deaths Baltimore XX X

Chicago X X X X XX X

Dallas/FW XXX X X

Los Angeles XXX X X X

Miami X X

New York X X X X X XX X

Philadelphia X XX X X

Pittsburgh X XX X X X

Seattle X X X X(90) X X

CVD death Baltimore X X

Chicago XXX X X X

Dallas/FW XXX X X

Los Angeles XXX X X

Miami X XX

New York X X X X X(90) X X X

Philadelphia X XX X X

Pittsburgh X XX X

Seattle X X X X(90) X X

Resp death Baltimore XXX X

Chicago XX X

Dallas/FW X X

Los Angeles X XX X X

Miami XXX

New York XX X

Philadelphia XX X

Pittsburgh X X(30)

Seattle XXX

The (number) following PM10 indicates a cumulative length used other than 200

days. Final model for each city may be different.

2218 M. Stylianou, M.J. Nicolich / Environmental Pollution 157 (2009) 2216–2223
Author's personal copy

Fig. 1 graphically displays the statistical significance levels of

PM10 and ozone in the above model. Neither PM10 nor ozone is

a predictor of total deaths in Miami, CVD deaths in Pittsburgh,

Baltimore, and Miami, or respiratory deaths in Pittsburgh, Miami,

and Dallas/Fort Worth.

However, both PM10 and ozone are associated with the total

deaths in Chicago, NY, Pittsburgh, and LA, CVD deaths in Chicago,

and respiratory deaths in NY and LA.

The highest percent variation explained by the above twopollutant model (for R2

adj>30%) with the outcome of total death

were NY, LA, Chicago, and Dallas/Forth Worth with 64%, 57%, 39%,

and 39%, respectively. (Table 3).

The three highest R2

adj with CVD outcome were LA, NY, and

Chicago with 54%, 53%, and 37%, respectively.

The best two R2

adj in the model with respiratory deaths were LA

and NY with 49% and 40%, respectively.

Adding either or both pollutants to the basic model results in

only slight (if any) model improvement. Overall, the model with

total mortality performed better than CVD or respiratory deaths,

although most cities had low R2

adj. Caution is necessary when using

R2 for quantifying predictive ability because in a large dataset, such

as NMMAPS, a factor that changes risk slightly may be statistically

significant but may not be useful in predicting the response of an

individual day.

3.3. RR cause of death

Fig. 2 (and Figs. S3 and S4 in the online supplement) display the

results of the unified analyses. The emerging pattern for PM10 and

ozone is that for a given outcome the shape and direction of the RR

functions are not consistent across cities. Consistent with our

previous analyses, RRs increase for some cities, decrease for others,

and are flat for the rest.

3.3.1. All deaths

The risk functions of PM10 on total deaths (Fig. 2) are mostly flat

near one, except Chicago where the RR for PM10 increases and

flattens around 35 mg/m3

. Some risk estimates are not increasing

with the PM10 level, such as NY or LA. Some even inexplicably

decrease, as in Pittsburgh or Baltimore. This behavior could justify

the monotonic (non-decreasing) assumption for the effect or risk

function. The association and shape of the adjusted PM10 RR function

(adjusted for ozone, day, time, temperature, and dp) are inconsistent

across the nine cities; the only risk function significantly larger than

one is for Chicago at 29 mg/m3

. Over half the functions appear linear

but are either flat or declining. Chicago clearly appears to be nonlinear. The observed range of PM10 average cumulative levels and

variability differ considerably between cities. Likewise, the shape of

the ozone risk function is also inconsistent; thresholds appear to

occur in the range of 20–40 ppb for all cities but Miami and LA.

The estimated ozone risk function is not always a non-decreasing

function (LA, Philadelphia, Pittsburgh, and Seattle).

3.3.2. CVD deaths

The PM10 risk estimates for CVD (Fig. S3) are non-linear for

Chicago, and LA where the risk function increases and then flattens

at about 30 mg/m3

. The risk functions for Dallas and Philadelphia

tend to increase linearly, and the rest are mostly flat near one. The

ozone risk-functions suggest thresholds at about 25–30 ppb. Dallas,

Philadelphia, and Pittsburgh seem to have a linear ozone riskfunction, but in Seattle and Miami, it is curvilinear mostly below

one, suggesting a beneficial effect, again a justification for a monotonic assumption.

3.3.3. Respiratory deaths

PM10 risk-functions for respiratory death (Fig. S4) are also

inconsistent. Some cities show a constant negative risk (Dallas, LA,

Miami, Seattle), some a risk near one (Baltimore, NY, Philadelphia,

Pittsburgh), and Chicago is suggestive of a hockey-stick shape with

an upward trend near 40 mg/m3

.

For ozone, most cities tend to have a risk that seems to be linear

except LA, NY, and Seattle that have a non-linear shape. A negative

risk seems to occur in LA, Miami, and Seattle for most of the

observed range.

Table 2

Estimated threshold levels by cause of death, pollutant and city based on the city-specific models of Table 1.

City Total Deaths CVD Deaths Respiratory Deaths

PM10 O3 SO2 NO2 PM10 O3 SO2 NO2 PM10 O3 SO2 NO2

Baltimore GAM RR – ND – – – – – – – 20 – –

m.RR – 39 – – – – – – – ND – –

Chicago GAM RR 29 25 4 22 30 38 – – 28 – – –

m.RR 26 11 12 1 30 50 – – ND – – –

Dallas/Fort Worth GAM RR– 30– – – – 5 – – – – –

m.RR – ND – – – – ND – – – – –

Los Angeles GAM RR ND – – 22 – 38 – – – – – 29

m.RR 43 – – 19 – 36 – – – – – 19

Miami GAM RR– – – 10– – – – – – – –

m.RR – – – ND – – – – – – – –

New York GAM RR ND 21 8 37 ND 20 8 30 – 39 – –

m.RR 24 2 7 10 ND 2 5 27 – 2 – –

Philadelphia GAM RR – 36 7 – – 16 7 – – 17 – –

m.RR – 35 1 – – 36 16 – – ND – –

Pittsburgh GAM RR ND 44 ND – – – – 22 ND – – –

m.RR ND ND ND – – – – 29 ND – – –

Seattle GAM RR ND 38 ND – ND ND ND – – – – –

m.RR ND 17 ND – ND ND ND – – – – –

Thresholds obtained based on the models listed in Table 1. Cells with ‘‘–’’ indicate that the pollutant was not associated with the outcome, where ‘‘ND’’ indicates that no

threshold is detected. The GAM RR row lists the concentration level at which the RR is significantly larger than one using the whole dataset. m.RR lists the concentration levels

at which RR departs one under the non-decreasing monotonic assumption (based on a random sample of the available data).

M. Stylianou, M.J. Nicolich / Environmental Pollution 157 (2009) 2216–2223 2219
Author's personal copy

4. Discussion

We studied the cumulative effect of PM10, O3, NO2, SO2, and CO

on CVD, respiratory, and total non-accidental deaths. Overall, we

found that the two most frequently death-associated pollutant

metrics are cumulative 200-day PM10 and 3-day cumulative ozone.

We believe that there is a biological plausibility related to the more

extended cumulative exposure of PM versus the gaseous pollutants

beyond the statistical related reasons. For PM, there is mechanistic

related information that suggests possible systemic toxic effects

including possible cardiovascular effects that include biological

events that may proceed over longer periods, i.e., not only acute

events such as changes in ‘‘heart rate variability’’ or increased

infarction but also chronic events such as contributing to atherosclerosis (Sun et al., 2005; Pope et al., 2004). In contrast, with the

exception of chronic inflammation, the mechanistic data for ozone

are dominated by acute respiratory related events (Costa, 2001).

The literature suggesting cardiac mechanisms for ozone is limited

(Gong et al., 1998). Further, in the premier chronic air pollution

epidemiology study, there was no statistically significant positive

association for ozone whereas there was such a signal for PM

(Pope et al., 2002). Hence, there is a biologic basis for examining

longer duration effects for PM versus ozone.

Our analyses suggest the effect of these pollutants across the

nine cities studied is inconsistent. Associations are more likely for

CVD or total death in older large industrialized cities such as

Chicago and NY. A possible reason as to why PM10 and ozone were

not as strongly associated with respiratory deaths could be the low

daily events. While high levels of ozone and PM10 are associated

with increased non-accidental mortality, our analyses suggest that

these effects are not uniform across the US, that dose–response

relationships are not necessarily linear, and that some relationship

shapes are suggestive of no-effect thresholds.

Higher average pollutant concentration levels do not necessarily

predict that the pollutant is associated with increased mortality. For

example, Miami has one of the highest average concentrations of

ozone but this pollutant is not associated with any cause of death.

Philadelphia had the second highest PM10 levels but this pollutant

is not associated with mortality. This finding contrasts with the

current belief that higher pollution levels are associated with

higher death counts. This appears to be the case for some cities, but

not others. We have no scientific or biologic explanation as to why

this might be, but it does indicate that these results are not

homogeneous across the different cities. Similar observations were

made by Janes et al. (2007), that is, decreasing levels of PM2.5 were

not associated with decreasing rates in mortality.

The cumulative effects of the pollutants, especially PM10, were

not associated with death events for some cities. When PM10 was

a significant predictor, the 200-day cumulative period was, generally, a better predictor than shorter durations.

We believe that a monotonic (but not necessarily linear)

assumption for a pollutant effect is reasonable in the absence of

plausible mechanisms for beneficial effects, and the likelihood that

effects at low concentrations are likely due to chance. RR estimates

using the monotonic assumption are based on relatively small

samples of the data, usually 15–20 percent depending on the

number of factors used for each analysis. Therefore, the constrained

analyses were less powerful than the unconstrained analyses.

However, repeating some of the analyses with different random

samples produced similar results.

The pollution averaging times used in these models differ from

the averaging times for the current standards. For example, current

PM standards are based on a 24-h or annual averages, while the

averaging times used here are primarily 200-day weighted means

with occasional 30 or 90-day weighted means (Table 4). Care

should be taken to not compare the thresholds and responses

presented in this paper to the pollutant averaging times used in the

NAAQS. In general, because short-term extremes flatten out over

longer periods, the average levels tend to be lower as the cumulative period increase. For example, one would expect the 8-h

average ozone level to be higher than that of daily average level

The RR is often assumed or reported to be linear (Bell et al.,

2004; Schwartz et al., 2001), implying that a specific increase of

a given pollutant has a fixed risk. In contrast, our cumulative

pollutant community-specific analyses indicate that the pollutant

RR is often non-linear and the rate of risk increase is not constant.

Bell et al. (2006) have explored the possibility of non-linearity or

thresholds, but they only reported the combined national average,

leaving the reader wondering about the shape of the communityspecific risk estimates. Even with the combined risk estimate they

did find the possibility of thresholds at low levels of ozone.

However, they dismissed the threshold models because the

improvement in the model AIC was minimal (<1%). However, we

have demonstrated (Table 3) that the inclusion of statistically

significant PM10 or ozone in the basic model does not provide

a considerable model improvement either (<0.4%). Our analyses

0

0.05

.1

.2

.3

.4

.5

O3 p−value

Total Deaths

Balt

Chic D/FW LA

Miam

NY Phil Pitt Seat

0

0.05

.1

.2

.3

.4

.5

O3 p−value

CVD Deaths

Balt

Chic

D/FW

LA

Miam

NY

Phil

Pitt

Seat

0 0.05 .1 .2 .3 .4 .5

0

0.05

.1

.2

.3

.4

.5

PM10 p−value

O3 p−value

Resp Deaths

Balt

Chic

D/FW

LA

Miam

NY Phil

Pitt

Seat

P−values of Ozone and PM10

(Values <.001 or >.5 are truncated to .001 and .5 respectively)

Fig. 1. Significance levels of PM10 and ozone in the unified model, by city.

2220 M. Stylianou, M.J. Nicolich / Environmental Pollution 157 (2009) 2216–2223
Author's personal copy

suggest that thresholds exist for some of the cities. We found that

PM10 risk-function for Chicago tends to increase and then flatten at

about a 200-day mean of 30 mg/m3 in contrast to ozone risk-functions for most cities where the risk is flat at low levels and then

starts to increase at about a 3-day mean of 20–30 ppb. These

conclusions remain the same, even when we run the same basic

model plus both PM10 and ozone. In another sensitivity analysis, we

have compared the model R2

adj and AIC with an ozone term as

either linear or with a threshold at 30 ppb (Table S1). For most

cities, the threshold model performs slightly better than the model

using the linear ozone term.

Many authors attempted to combine city-specific RRs into

a single national estimate (Bell et al., 2004, 2006; Dominici et al.,

2002). This approach can increase the statistical power and give an

overall RR perspective. However, several underlying assumptions

must be true for the combined results to be valid. For example, all

the cities must be comparable in terms of other confounders, such

as location or socioeconomic variables. Some of these confounders

may be partly addressed by the weather (for location) and calendar

time adjustments. Since RR estimates are conditional on other

predictors in the model, it is necessary to include the same variables in order to have the same risk interpretation. Our analyses

suggest that even with the same variables in the model, the shapes

of the risk-functions from such models are inconsistent across the

different cities. This indicates that a rigorous examination of the

results from each city is necessary before including specific cities

into combined estimates of the RR. Similar conclusions regarding

associations between 12-month PM2.5 exposure and mortality were

recently reached by Janes et al. (2007), that is, the association

between trends in PM2.5 and mortality at the national scale is more

likely to be confounded. On the other hand our conclusions about

linearity and thresholds are in contrast to those of Schwartz et al.

(2008) using long term PM2.5 exposure and mortality from the

Harvard Six Cities Study.

We note that some final models developed for each city did not

have high predictive power (R2

adj), which is also true for most

published results. As suggested by a reviewer, larger cities having

higher R2

adj may indicate that the daily count is less susceptible to

a random variation and therefore more amenable to explanation.

Each individual term is significantly associated with the response

variable but other factors not considered or accounted for might

contribute significantly to the prediction of these death events.

Our analyses are exploratory because of the computational

restrictions. As with Janes et al. (2007) conclusions, we believe that

further research in needed to explore the potential confounders

that may vary at different cities such as the possible meteorological

confounding. Additional sensitivity analyses could be performed

using the whole NMMAPS data, though one could not do all the

sensitivity checks possible. A limitation of our study is the fact that

we have imputed all missing measurements which we believe

resulted in an overestimation of relative risks such as that of ozone

levels for the winter time. Studying pollutants effect during the

summer months would give us a better understanding of their

effects. Based on our experience with the data we analyzed we feel

the results and conclusions will be similar, but with more insight

into the inter-city differences.

5. Conclusions

We conclude that there is not a uniform effect, or model, for

PM10 and O3 across the nine major US cities we studied; we believe

this precludes the use of a single national pollution model. Because

Table 3

Comparison of R2

adj * and AICy for PM10 and ozone when added individually to the basic model.

Event City ‘‘Basic’’ model þPM10 þO3 þPM10 þ O3

R2

adja AICb R2adj AIC diff R2adj AIC diff R2adj AIC diff

Total deaths Baltimore 15.1 28820 15.1 1 15.2 4 15.2 3

Chicago 37.9 37877 38.0 12 38.4 52 38.6 62

Dallas/FW 38.9 33922 38.9 1 39.0 9 39.0 8

Los Angeles 56.7 39157 56.8 3 56.9 23 57.0 29

Miami 12.0 32944 12.0 1 11.9 1 12.1 4

New York 63.1 41188 63.2 18 64.0 170 64.1 178

Philadelphia 28.3 32776 28.3 0 29.0 45 28.9 44

Pittsburgh 21.7 31992 21.8 4 22.0 16 22.1 20

Seattle 15.0 30457 15.0 0 15.2 8 15.2 7

CVD deaths Baltimore 9.5 23853 9.6 0 9.8 4 9.8 2

Chicago 36.5 33512 36.7 10 36.9 28 37.1 34

Dallas/FW 20.7 29664 20.7 1 20.7 0 20.7 1

Los Angeles 53.4 35305 53.5 2 53.5 6 53.6 5

Miami 9.2 28752 8.9 8 8.8 3 8.9 7

New York 52.1 36902 52.1 0 52.7 76 52.7 79

Philadelphia 21.5 28204 21.5 1 21.7 14 21.7 11

Pittsburgh 17.5 27934 17.5 0 17.6 3 17.6 3

Seattle 7.1 25661 7.1 1 7.2 5 7.2 4

Resp deaths Baltimore 4.0 14976 4.1 1 4.3 14 4.3 12

Chicago 22.4 24999 22.5 0 22.5 0 22.6 1

Dallas/FW 26.7 21093 26.7 1 26.7 2 26.7 3

Los Angeles 48.8 26895 48.9 3 48.9 5 49.0 7

Miami 7.0 19658 7.0 2 7.0 0 7.0 1

New York 39.2 26878 39.5 12 39.5 21 39.8 34

Philadelphia 12.1 18808 12.1 2 12.2 4 12.1 3

Pittsburgh 12.0 19182 12.0 1 12.0 1 12.0 2

Seattle 10.3 18357 10.4 0 10.5 10 10.6 10

‘‘Basic model’’ consists of the same terms for all cities. ‘‘þPM10’’ or ‘‘þO3’’ shows the effect of PM10 or O3 when added individually to the basic model. The last two columns

correspond to model (3).

a R2

adj is the proportion of explained variation, adjusted for the number of terms in the model; larger is better. b AIC is related to the reduction in the likelihood function, adjusted for the number of terms in the model; smaller is better. AIC diff is the AIC difference from the basic model.

M. Stylianou, M.J. Nicolich / Environmental Pollution 157 (2009) 2216–2223 2221
Author's personal copy

Fig. 2. Relative risk functions and 95% confidence intervals for total mortality associated with PM10 and ozone levels for individual cities using the unified model.

2222 M. Stylianou, M.J. Nicolich / Environmental Pollution 157 (2009) 2216–2223
Author's personal copy

many models exhibited thresholds (25–45 mg/m3 for PM10; and

10–45 ppb for O3) under a variety of assumptions including longterm cumulative effects, we feel that the use of a linear nonthreshold pollutant model is inconsistent with the observed data.

Acknowledgements

The authors thank Mr. Larry Gephart and the referees for their

constructive comments that significantly improved the presentation of this paper. Funding provided by the American Petroleum

Institute (API) and ExxonMobil Biomedical Sciences, Inc.

Appendix. Supplementary information

Supplementary data associated with this article can be found, in

the online version, at doi:10.1016/j.envpol.2009.04.011

References

Bell, M.L., McDermott, A., Zeger, S.L., et al., 2004. Ozone and short-term mortality in

95 US urban communities, 1987–2000. Journal of the American Medical

Association 292, 2371–2378.

Bell, M.L., Dominici, F., Samet, J.M., 2005. A meta-analysis of time-series studies of

ozone and mortality with comparison to the national morbidity, mortality, and

air pollution study. Epidemiology 16 (4), 436–445.

Bell, M.L., Peng, R.D., Dominici, F., 2006. The exposure–response curve for ozone and

risk of mortality and the adequacy of current ozone regulations. Environmental

Health Perspectives 114 (4), 532–536.

Costa, D., 2001. Air pollution. In: Klaassen, Curtis D (Ed.), Casarett and Doull’s

Toxicology, the Basic Science of Poisons. McGraw-Hill Publishing, New York.

Daniels, M., Dominici, F., Zeger, S.L., et al., 2000. Estimating PM10-mortality doseresponse curves and threshold levels: an analysis of daily time-series for the 20

largest US cities (with discussion). American Journal of Epidemiology 152 (5),

397–412.

Dockery, D.W., Pope, C.A., Xu, X., 1993. An association between air pollution and

mortality in six US cities. New England Journal of Medicine 329, 1753–1759.

Dominici, F., McDermott, A., Zeger, S.L., et al., 2002. On the use of generalized

additive models in time-series studies of air pollution and health. American

Journal of Epidemiology 156 (3), 193.

Eilers, P.H.C., Marx, B.D., 1996. Flexible smoothing with B-splines and penalties.

Statistical Science 11, 89–121.

Gong, H., Worn, R., Sarma, R.J., et al., 1998. Cardiovascular effects of ozone exposure

in human volunteers. American Journal of Respiratory and Critical Care Medicine 158, 538–546.

Janes, H., Dominici, F., Zeger, S.L., 2007. Trends in air pollution mortality, an

approach to the assessment of unmeasured confounding. Epidemiology 18 (4),

416–423.

Koop, G., Tole, L., 2006. An investigation of thresholds in air pollution-mortality

effects. Environmental Modelling & Software 21, 1662–1673.

Leitenstorfer, F., Tutz, G., 2007. Generalized monotonic regression based on

B-splines with an application to air pollution data. Biostatistics 8 (3), 654–673.

Miller, K.A., Siscovick, D.S., Sheppard, L., et al., 2007. Long-term exposure to air

pollution and incidence of cardiovascular events in women. New England

Journal of Medicine 356, 447–458.

Poole, C., 2001. Low P-values or narrow confidence intervals: which are more

durable? Epidemiology 12 (3), 291–294.

Pope III, C.A., Burnett, R.T., Thun, M.J., et al., 2002. Lung cancer, cardiopulmonary

mortality, and long-term exposure to fine particulate air pollution. Journal of

the American Medical Association 287, 1132–1141.

Pope III, C.A., Burnett, R.T., Thurston, G.D., et al., 2004. Cardiovascular mortality and

long-term exposure to particulate matter air pollution: epidemiological

evidence of general pathopsyiological pathways of disease. Circulation 109,

717–777.

Ren, C., Williams, G.M., Mengersen, K., 2008. Does temperature modify short-term

effects of ozone on total mortality in 60 large eastern US communities? – an

assessment using the NMMAPS data. Environment International 34, 451–458.

Roberts, S., Martin, M.A., 2006. Applying a moving total mortality count to the cities

in the NMMAPS database to estimate the mortality effects of particulate matter

air pollution. Occupational and Environmental Medicine 63, 193–197.

Samet, J.M., Dominici, F., Curriero, C., et al., 2000. Fine particulate air pollution and

mortality in 20 US cities, 1987–1994. New England Journal of Medicine 343,

1742–1749.

Schwartz, J., 1999. Air pollution and hospital admissions for heart disease in eight

US countries. Epidemiology 10, 17–22.

Schwartz, J., 2000. Mortality displacement and long-term exposure effects related

to air pollution and mortality. Part I. Health Effects Institute Research Report 94,

51–59.

Schwartz, J., Ballester, F., Saez, M., et al., 2001. The concentration–response relation

between air pollution and daily deaths. Environmental Health Perspectives 109

(10), 1001–1006.

Schwartz, J., Coull, B., Laden, F., et al., 2008. Environmental Health Perspectives 116

(1), 64–69.

Sun, Q., Aixia, W., Ximei, J., et al., 2005. Long-term air pollution exposure and

acceleration of atherosclerosis and vascular inflammation in an animal model.

Journal of the American Medical Association 21 (294), 3003–3009.

Weinberg, C.R., 2001. It’s time to rehabilitate the P-value. Epidemiology 12 (3),

288–290.

Xia, Y., Tong, H., 2006. Cumulative effects of air pollution on public health. Statistics

in Medicine 25 (20), 3548–3559.

Zeka, A., Schwartz, J., 2004. Estimating the independent effects of multiple

pollutants in the presence of measurement error: an application of

a measurement-error-resistant technique. Environmental Health Perspectives

112 (17), 1686–1690.

Table 4

Averaging times currently used in the NAAQS and in the current paper.

Averaging Times

NAAQSa Current Paper

PM10 24 h mean 200-day weighted mean

Annual mean (or 30 or 90-day weighted mean)

O3 1 h mean 3-day weighted mean

8 h mean

NO2 Annual mean 3-day weighted mean

SO2 3 h mean 3-day weighted mean

24 h mean

Annual mean

CO 1 h mean 3-day weighted mean

8 h mean

a All means are arithmetic means.

M. Stylianou, M.J. Nicolich / Environmental Pollution 157 (2009) 2216–2223 2223