Time-Series Analyses of Air Pollution and Mortality in the United States: A Subsampling Approach

Background: Hierarchical Bayesian methods have been used in previous papers to estimate national mean effects of air pollutants on daily deaths in time-series analyses. Objectives: We obtained maximum likelihood estimates of the common national effects of the criteria pollutants on mortality based on time-series data from ≤ 108 metropolitan areas in the United States. Methods: We used a subsampling bootstrap procedure to obtain the maximum likelihood estimates and confidence bounds for common national effects of the criteria pollutants, as measured by the percentage increase in daily mortality associated with a unit increase in daily 24-hr mean pollutant concentration on the previous day, while controlling for weather and temporal trends. We considered five pollutants [PM10, ozone (O3), carbon monoxide (CO), nitrogen dioxide (NO2), and sulfur dioxide (SO2)] in single- and multipollutant analyses. Flexible ambient concentration–response models for the pollutant effects were considered as well. We performed limited sensitivity analyses with different degrees of freedom for time trends. Results: In single-pollutant models, we observed significant associations of daily deaths with all pollutants. The O3 coefficient was highly sensitive to the degree of smoothing of time trends. Among the gases, SO2 and NO2 were most strongly associated with mortality. The flexible ambient concentration–response curve for O3 showed evidence of nonlinearity and a threshold at about 30 ppb. Conclusions: Differences between the results of our analyses and those reported from using the Bayesian approach suggest that estimates of the quantitative impact of pollutants depend on the choice of statistical approach, although results are not directly comparable because they are based on different data. In addition, the estimate of the O3-mortality coefficient depends on the amount of smoothing of time trends.

Time-series analy ses of air pollution and various health end points, including daily mortality, using flexible generalized linear models (GLMs) or generalized additive models (GAMs), have become commonplace in the last decade. The National Morbidity, Mortality and Air Pollution Study (NMMAPS) was an ambitious effort undertaken by scientists at Johns Hopkins University (Baltimore, MD, USA) and Harvard University (Cambridge, MA, USA) to investigate the association between particulate pollution (i.e., particulate matter with an aerodynamic diameter ≤ 10 µm; PM 10 ) and morbidity and mortality in the 90 largest metropolitan areas in the United States over the period 1987-1994 using time-series methods (Dominici et al. 2003;Samet et al., 2000aSamet et al., , 2000bSamet et al., , 2000c. For the mortality analy ses, the investigators used a common approach for analy ses of the timeseries data on daily PM 10 levels and deaths for each city. Specifically, a national mean estimate for the association between PM 10 and mortality was obtained by combining the individual city-specific estimates using a hierarchical Bayes procedure, which assumed that the city-specific effects were normally distributed. Since then the analy ses have been extended to include more cities, more years of data, and a second pollu tant, ozone (O 3 ) (Bell et al. 2004;Smith et al. 2009). In this article, we propose an approach to analyzing multi city time-series data that is complementary to the hierarchical Bayes approach.
Arguably, for standard setting, one might want to estimate not a national mean effect, but a common pollu tant effect across the country. The most direct way to estimate a common pollu tant effect on mortality across multiple cities would be to analyze all cities simultaneously using a GLM or GAM model postulating a common pollu tant effect estimate across cities, but with control of confounders specific to each city. A statistical test for the hypothesis that there is a common pollu tant effect across cities could then be based on standard likelihood-based procedures or on the Akaike or Bayes information criterion (AIC or BIC). Alternatively, a computational approach to investigating the distribution of a common estimator of a pollu tant effect is to use the bootstrap or the jackknife procedure on the full complement of cities (Efron and Tibshirani 1993). However, with a large number of cities this conceptually simple approach presents formidable computational problems. One way around the computational problems is to use either the "delete-k jackknife" followed by the bootstrap procedure (Efron and Tibshirani 1993) or the subsampling procedure (Politis et al. 1999). In this study, we used the subsampling procedure to analyze the time-series data on all the criteria pollu tants, with the exception of lead, and mortality in 108 metropolitan areas in the United States over the 14-year period 1987-2000.

Methods
We downloaded the mortality and air pollution time-series data from the NMMAPS web site maintained by the Johns Hopkins investigators (http://www.ihapss.jhsph.edu/data/ data.htm). Using the same database allowed direct comparisons to be made to the previous work of the Johns Hopkins team and others using these data. Daily data on the number of deaths are availa ble for 108 metropolitan areas over the 14-year period 1987-2000. Daily concentrations of the criteria pollu tants-PM 10 , O 3 , carbon monoxide (CO), nitrogen dioxide (NO 2 ), and sulfur dioxide (SO 2 )-are also availa ble; although for each of the pollutants, information is availa ble only for a subset of the days. A limitation, shared with all epidemiologic analy ses of air pollution data, is that the data were collected from stationary monitors to determine compliance with air quality regulation. Thus, the data collected are not necessarily ideal for epidemiologic studies. We recognize also that ambient concentrations are imperfect surrogates for personal exposure.
For the usual bootstrap approach, for each bootstrap cycle, 108 cities would be chosen with replacement from the original set of Background: Hierarchical Bayesian methods have been used in previous papers to estimate national mean effects of air pollu tants on daily deaths in time-series analy ses. oBjectives: We obtained maximum likelihood estimates of the common national effects of the criteria pollu tants on mortality based on time-series data from ≤ 108 metropolitan areas in the United States. Methods: We used a subsampling bootstrap procedure to obtain the maximum likelihood estimates and confidence bounds for common national effects of the criteria pollu tants, as measured by the percentage increase in daily mortality associated with a unit increase in daily 24-hr mean pollu tant concentration on the previous day, while controlling for weather and temporal trends. We considered five pollu tants [PM 10 , ozone (O 3 ), carbon monoxide (CO), nitrogen dioxide (NO 2 ), and sulfur dioxide (SO 2 )] in single-and multi pollu tant analy ses. Flexible ambient concentrationresponse models for the pollu tant effects were considered as well. We performed limited sensitivity analy ses with different degrees of freedom for time trends. results: In single-pollu tant models, we observed significant associations of daily deaths with all pollu tants. The O 3 coefficient was highly sensitive to the degree of smoothing of time trends. Among the gases, SO 2 and NO 2 were most strongly associated with mortality. The flexible ambient concentration-response curve for O 3 showed evidence of non linearity and a threshold at about 30 ppb. conclusions: Differences between the results of our analy ses and those reported from using the Bayesian approach suggest that estimates of the quantitative impact of pollu tants depend on the choice of statistical approach, although results are not directly comparable because they are based on different data. In addition, the estimate of the O 3 -mortality coefficient depends on the amount of smoothing of time trends. volume 121 | number 1 | January 2013 • Environmental Health Perspectives 108 cities, and the maximum likelihood estimate (MLE) of the pollu tant effect computed using Poisson regression methods. This procedure would be repeated many times to obtain an estimate of the distribution of the common pollu tant effect. However, with the large number of cities considered in these analy ses, it is computationally infeasible to derive an MLE using this approach. Therefore, an alternative procedure, such as the delete-k jackknife (Efron and Tibshirani 1993) or the closely related subsampling procedure (Politis et al. 1999) must be implemented. For the delete-k jackknife followed by the bootstrap procedure, first the set of all possible subsets of size d = 108 -k is constructed from the 108 cities. Each bootstrap cycle then randomly selects a city from this set, and the common pollu tant effect is estimated for the specific chosen set of d cities.
We used a closely related method, the subsampling procedure described by Politis et al. (1999). In this procedure, for each bootstrap cycle, we randomly chose d cities without replacement out of the 108 cities with availa ble data, and estimated the common pollu tant effect for each sample of d cities. Politis et al. (1999) recommend that d be much smaller than 108, and we chose d = 4. The choice of d is arbitrary; however, the confidence intervals (CIs) for the parameter estimates have to be adjusted for this choice as discussed below. The distribution of the estimator of the common pollu tant effect was based on 5,000 bootstrap cycles.
Let X 1 , X 2 , …, X n be a sequence of independent observations (i.e., a sequence of realiza tions of independent random variables). In our case, each observation in the sequence represents the data on daily deaths, pollu tant concentrations, and weather variables over the period 1987-2000 in each of the n cities considered in the analy ses. The total number of cities, n, depends on the pollu tant or combination of pollu tants considered. For example, n = 102 for analy ses that involve PM 10 alone, and n = 56 for analy ses that involve all pollutants. Let Θ be the parameter representing the common effect of a pollu tant, and let Θ n be the MLE for Θ based on n observations and Θ d the MLE for Θ based on d<<n observations (i.e., based on a subset of d cities), then, by general maximum likelihood theory, Θ n and Θ d converge to Θ with rates √ n and √ d, respectively.
Under standard likelihood-based procedures, inferences would be based on the distribution of Θ n . However, computation of Θ n and its distribution is infeasible with large n. Therefore, we based our inferences on Θ d , where in our case, d = 4, and based our CIs on the subsampling distribution as described by Politis et al. (1999). The properties of the subsampling procedure hold under rather weak conditions. The sequence of observations is required to satisfy the α-mixing condition (Politis et al. 1999, p. 315). This condition is seen to be trivially satisfied by a sequence of independent observations. An anonymous referee suggested that sampling units other than entire cities might be more appropriate. We agree that other sampling schemes need to be explored.
CIs were computed after adjusting for the size of the subsample as follows. If F d * represents the empirical distribution function of √ d (Θ d -Θ), then for any significance level α, with probability close to 1 -α. It then follows that In our case, we were unable to estimate Θ n because of computational issues. We therefore approximated Θ n by (Σ Θ d )/N, where N = 5,000, the total number of subsamples drawn. That is, we approximated Θ n as the mean of the Θ d .
Thus, for our analy ses, for each bootstrap cycle, we drew a random sample of four cities without replacement from among the 108 cities with availa ble data. We then fit an over-dispersed Poisson model to the randomly chosen 4 cities to obtain the MLE of the common pollu tant effects on mortality in the 4 cities, but with confounders, such as temperature and relative humidity, being separately controlled in each of the 4 cities. Because a number of previous analy ses (e.g., Bell et al. 2004;Dominici et al. 2003) have considered the effect of the pollu tants with a 1-day lag, we have done the same in these analy ses. Likewise, the number of degrees of freedom (df) for time trends and weather are also consistent with those used in previous analy ses (Bell et al. 2004;Dominici et al. 2003;Moolgavkar 2003). Specifically, for each bootstrap cycle we modeled the number of deaths from all causes (with accidents and suicides removed) in a city on a specific day as a function of the 24-hr average pollu tant concentration on the previous day, temporal trends (50 or 100 df natural spline), day of the week (categorical variable), mean temperature on the previous day (6 df natural spline), and mean dew-point temperature on the previous day (6 df natural spline). Note that this model controlled for confounders, day of week effects, and time trends in a city-specific fashion.
To investigate the shape of the ambient concentration-response relationship, we used the same models with the pollu tant effects represented by natural splines with 6 df. Multiple pollu tants that have been concurrently measured can be easily added as covariates in these analy ses. Because of missing data, the size of the data set from which samples of cities were drawn for analy ses depended on the number of pollu tants considered.
Clearly, the choice of lags for specific pollu tants and for the weather variables-temperature, and relative humidity-should be based on biological considerations whenever possible. Unfortunately, there is little information to guide these choices. In previous publications, a 1-day lag has often been used. For example, in the revised NMMAPS analy ses, Dominici et al. (2003) used the same lag for each of the pollu tants in their multi pollu tant analy ses. Their results indicate that a 1-day lag yields close to the maximum impact on daily mortality. It would be possible also to consider other lag structures, such as distributed lag models. However, the purpose of this study was not to be a comprehensive re analysis of the NMMAPS data, but to provide an approach to the analysis of national data that complements the Bayesian approach. Similarly for the weather covariates, we used a 1-day lag because there is little biological information to suggest that any specific lag structure is better than any other.
The model was fit to the data using the R software package (http://www.r-project.org). The means of the 5,000 maximum likelihood estimates of the common pollu tant effects were approximately unbiased and consistent estimators of the common national effects of the pollu tants.
We conducted simulations to investigate the coverage properties of the CIs constructed as described above. Specifically, we generated 100 observations based on a Poisson variate with an intercept and a slope. From these 100 observations, we drew 5,000 subsamples of size 4 without replacement and computed 90% and 95% CIs as described above in addition to computing the usual likelihood-based CIs. We repeated this entire procedure 1,000 times to investigate the coverage properties of the CIs. We found that the mean of the 5,000 subsample estimates was an excellent approximation to the MLE (data not shown). While the usual likelihood-based CIs covered the true values of the parameters with the nominal coverage probabilities, the subsampling CIs were conservative, that is, their coverage probabilities were larger than the nominal coverage probabilities. The 95% CIs covered the true values of the parameters approximately 98% of the time and the 90% CIs covered the true values of the parameters approximately 95% of the time (data not shown). Therefore, because the estimated CIs were highly conservative (i.e., too wide), we present both 95% and 90% CIs in our tables in the present study and note that our tests of significance (i.e., whether the CIs contain 0) are also conservative (i.e., the actual level of significance is smaller than the nominal alpha). The coverage properties of both of the CIs constructed using the subsampling approach and the credible intervals constructed in the hierarchical Bayes approach need to be explored in realistic simulation scenarios.

Results
Bootstrap means and 90% and 95% CIs after small-sample corrections for common nationwide estimates of effects of pollu tants represent associations with incremental changes in the 24-hr average concentration on the previous day of 10 µg/m 3 for PM 10 ; 10 ppb for O 3 , SO 2 , and NO 2 ; and 1 ppm for CO. The unit measures for the individual pollu tants were chosen to facilitate comparisons with previous estimates rather than typical day-to-day variations in ambient concentrations, which vary among cities. With time trends smoothed using 100 df natural splines, all pollu tants were significantly (at the 0.05 level) associated with mortality in single-pollu tant models (Table 1). With 50 df natural splines for time trends, the estimated coefficient for O 3 was greatly attenuated and statistically insignificant.
The magnitudes of associations with incremental increases in CO, NO 2 , and SO 2 were greater than for associations with PM 10 and O 3 . When CO, NO 2 , and SO 2 were included in the same model with 100 df splines for time trends, all three associations were attenuated and the association of CO with mortality was no longer statistically significant (Table 2).
Two-pollu tant analy ses of PM 10 and each of the gases with 100 df splines for time trends indicated significant associations for PM 10 in all four models, although estimated effects were attenuated except in the model adjusted for O 3 (Table 3). Effect estimates for CO, NO 2 , and SO 2 were significant and consistent with estimates that were not adjusted for PM 10 , but the estimated effect of O 3 on mortality was attenuated and was not statistically significant based on the 95% CI or even the 90% CI (Table 3).
Bootstrap means for flexible ambient concentration-response relationships, using 6 df natural splines, between pollu tant concentrations and deaths on the following day estimated from single-pollu tant models suggest non linearity and threshold-like behavior for NO 2 , PM 10 , and O 3 (Figure 1). However, the CIs for NO 2 and PM 10 are wide and the concentration-response relationships are consistent with linearity.
Our approach assumes a common national pollu tant effect. However, if there is a common national pollu tant effect, then for a randomly chosen set of four cities, one would expect a common shared-effects model to have a lower AIC than a model postulating separate effects. Therefore, for each bootstrap cycle, we computed the AIC for the common pollu tant effect model and for the individual models of each sample of four cities selected for each cycle. For all pollu tants, the AIC was lower for the separate four-city fits than for the model postulating common pollu tant effects for more than half the cycles (data not shown).

Discussion
In this study, we used a subsampling bootstrap approach to estimating maximum likelihood estimates for common national effects for associations of individual pollu tants with mortality. Previous analy ses of NMMAPS data have used a hierarchical Bayes approach (e.g., Bell et al. 2004;Samet et al. 2000aSamet et al. , 2000bSamet et al. , 2000cSmith et al. 2009) to estimate a national mean effect.
The procedure described in this study is one approach to estimating common pollutant effects, if they exist; however, in a vast geographically and climatically diverse country such as the United States, it is not unreasonable to expect heterogeneity of pollu tant effects across the country. First, PM is a complex mixture whose composition varies by region and season. O 3 could also be considered a mixture because it is generally present with other oxidants. Second, for all pollutants, any effects on human health would be expected to be modified by weather and by the circumstances of exposure, which clearly vary by region and season. A comparison of AIC statistics for sample-specific (four city) models and for common estimate models suggests that the four-city models fit the data better than the common model in most cases, which is consistent with heterogeneous effects among cities. This finding is consistent with the results of the hierarchical Bayes multi city analy ses (Bell and Dominici 2008;Dominici et al. 2003;Smith et al. 2009). The heterogeneity of pollutant effects suggests that any single national estimate may not provide a reliable measure of the health benefits that would accrue from a reduction in pollu tant concentrations.
An advantage of the Bayesian approach is that it allows for heterogeneity of cityspecific coefficients, albeit with the simplifying assumptions that these are independent and identically distributed. That said, the hier archical Bayes approach and the subsampling approach described in this study are complementary approaches to estimating a single national pollu tant effect. Under both approaches, however, further investigation of heterogeneity requires that regional analy ses be performed.
The use of Bayesian methods in multi city analy ses was pioneered by the first investigators of NMMAPS (Dominici et al. 2000;Samet et al. 2000aSamet et al. , 2000bSamet et al. , 2000c. Their approach to the analy ses of multi city data used a two-stage procedure, with analyses of  Units are 10 ppb for NO 2 and SO 2 and 1 ppm for CO. Temperature and relative humidity on the previous day are controlled using 6 df natural splines; time trends are controlled using 100 df natural splines and day of week is controlled as a categorical variable. These joint-pollu tant analy ses are based on data from 58 cities. volume 121 | number 1 | January 2013 • Environmental Health Perspectives single cities at the first stage, followed by a hierarchical Bayes analysis of the first-stage results to arrive at a single estimate of a national mean for the pollu tant effect. Once the first-stage analy ses were completed, the estimated pollu tant effects were considered in isolation from the other covariates and combined in a second stage, resulting in a procedure that is operationally similar to a meta-analysis. This procedure assumes that the asymptotic distribution of the pollu tant effect estimate has been achieved, and that a simple approximation to a full Bayesian analysis yields valid results in the second stage. The first assumption regarding the asymptotic distribution of parameter estimates is widely made in statistical procedures; in fact, one study (Dominici et al. 2000) reports that the asymptotic properties of the MLE are well approximated in the first-stage analy ses. Our procedure also relied on asymptotic results. The assumptions underlying the secondstage procedure may be more problematic. A fully Bayesian analysis would put prior distributions on all the parameters of the model, not just the pollu tant coefficients, and allow the parameter estimates along with their covariance structure from the first stage to be carried forward to the second-stage analy ses. Although this procedure may be possible in theory, it would be much more computationally intensive and may not be practical. Although it is not clear that a fully Bayesian analysis would make a substantial difference to the results, these issues do not arise with the approach proposed in this study, and common national pollu tant-mortality coefficients and their CIs can be estimated in a single-stage analysis even when multiple pollu tants are modeled using splines. However, the simple simulation that we performed indicated that the CIs generated using our approach were conservative. The coverage properties of CIs and credible intervals generated, respectively, by the subsampling and the hierarchical Bayes procedures need to be investigated using realistic and comprehensive simulation scenarios. Smith et al. (2009) observed that the national mean from the hierarchical Bayes analysis approximates the mean (not the inverse variance weighted mean) of the individual city-specific estimates as the variance of assumed distribution of the city-specific coefficients goes to infinity. They cautioned against any use of a national statistic, but advocated the use of a population-weighted mean if a national statistic is computed. The mean of the distribution of the bootstrap samples we generated implicitly has a population weight built in because the common estimate for each bootstrap cycle is influenced by the sizes of the populations of the cities in that cycle.
Our CIs for the estimated common national effects of the pollu tants are wider than the credible intervals reported in earlier studies using the hierarchical Bayes approach. This could be a consequence of hierarchical Bayes credible intervals that may be too narrow (because the fully Bayesian procedure is replaced by an approximation that ignores estimates of weather and time-trend parameters in the second-stage analy ses), in addition to subsampling CIs that are too wide (conservative), as suggested by our simulations.
Previous multi city analy ses have focused on PM 10 or O 3 (Bell et al. 2004(Bell et al. , 2006Dominici et al. 2003;Katsouyanni and Samet 2009;Smith et al. 2009). Analyses that addressed regional heterogeneity (Bell and Dominici 2008;Dominici et al. 2003;Smith et al. 2009) reported considerable hetero geneity of estimated coefficients for both PM and O 3 . Dominici et al. (2003) reported the largest PM 10 coefficients were in the Northeast region among the seven regions they considered, followed by Southern California. The smallest PM 10 coefficients were reported in the Upper Midwest. The estimated coefficients in the Northeast were twice as high as those reported for the Upper Midwest. Smith et al. (2009) estimated the highest O 3 coefficients for the Northeast and the Industrial Midwest but did not estimate significant effects for Southern   (Moolgavkar 2003). Previous time-series analy ses of air pollution and mortality have focused on PM 10 or O 3 , with other pollu tants addressed as confounders, if at all (Dominici et al. 2003). Our findings suggest that the emphasis on PM and O 3 may deserve reconsideration because estimated associations were strongest for CO, NO 2 , and SO 2 . In addition, coeffi cients for all three gases remained highly significant in joint pollu tant analy ses with PM 10 , and coefficients for NO 2 , and SO 2 remained significant in a model that included all three gases.
PM 10. The NMMAPS analy ses estimated a national mean increase of 0.27% in mortality for a 10-µg/m 3 increase in PM 10 on the previous day using GAMs, and a 0.22% increase using a GLM. Our estimate of the mean with 50 df natural splines for time trends is 0.4% (95% CI: 0.30, 0.53). This estimate is not directly comparable to the estimates from the hierarchical Bayes analy ses because we used somewhat different city-specific models and had many more cities with more years of data, although it is of interest to note that our estimate is virtually identical to the estimate reported for the Northeast (0.41; 95% posterior interval: 0.04, 0.78) in the NMMAPS re analy ses (Dominici et al. 2003).
In two-pollu tant analy ses of PM 10 with CO, NO 2 , or SO 2 , the PM 10 coefficients were somewhat attenuated but continued to be significantly associated with mortality. The association with PM 10 was attenuated most in joint pollu tant analy ses with NO 2 . With the exception of the two-pollu tant model of PM 10 and O 3 , which produced a highly insignificant coefficient for O 3 when a 50 df smoother was used, results were robust to the degree of smoothing of time trends (data not shown). Therefore, we report only 100-df time-trend smoothers for joint pollu tant analy ses in this study.
O 3 . There are several fairly recent metaanaly ses (Bell et al. 2005;Ito et al. 2005;Levy et al. 2005) and multi city analy ses (Bell et al. 2004;Katsouyanni and Samet 2009;Smith et al. 2009) of the association between O 3 and daily deaths. As Smith et al. (2009) point out, there are two kinds of potential biases in meta-analy ses: publication bias and model selection bias (in which the investigators test many models and report the results of only those models that show positive and significant associations). Therefore, multi city analyses are preferable to meta-analy ses.
The first multi city analy ses were conducted by Dominici et al. (2003). These analy ses were based on 8 years of data in 90 cities in the United States. For a 10-ppb increase in daily O 3 , Dominici et al. (2003) reported a national mean increase of approximately 0.25% in deaths on the following day, and this estimate was statistically significant. In an update of these analy ses in 95 cities and with 14 years of data, Bell et al. (2004) reported a statistically significant national mean of approximately 0.2% (approximate 95% posterior interval: 0.1, 0.3%) for a 1-day lag O 3 effect on mortality. In a reanalysis of data from the NMMAPS (Katsouyanni and Samet 2009), the O 3 mortality estimates were highly sensitive to the degree and type of smoothing used for seasonality control, and inclusion of PM 10 markedly reduced the O 3 mortality estimates. Our estimated mean with 100 df smoothers for time trends is 0.40% (95% CI: 0.27, 0.56). With 50 df smoothers for time trends, the estimated mean for O 3 is 0.08% (95% CI: -0.16, 0.38). Thus, the estimated O 3 effect is sensitive to the degree of smoothing of temporal trends, in agreement with the results reported by Katsouyanni and Samet (2009). Because O 3 is a highly seasonal pollu tant, adequate control of seasonality is important, and using only 50 df for time trends may not be sufficient. To investigate this issue more thoroughly, further analy ses of the data stratified by season should be undertaken.
In joint pollu tant analy ses (100 df smoothers for time trends) with PM 10 , the O 3 coefficient is substantially attenuated and becomes insignificant, as also reported by Smith et al. (2009) and Katsouyanni and Samet (2009). One possible reason for the widening of CIs for O 3 effects is the greatly reduced data set used in joint pollu tant analy ses. If this were the sole reason, however, one would expect to observe the same phenomenon for the PM 10 effect, which was virtually unchanged and continued to remain highly significant.
In an update to the Bell et al. (2004) study, Bell et al. (2006 examined the ambient concentration-response curve for O 3 and mortality, where their O 3 metric was the average of the 24-hr mean concentrations on the current and previous days. Among various approaches evaluating ambient concentration-response relationships, Bell et al. (2006) used natural splines followed by a hierarchical Bayes procedure to obtain a national estimate. The results shown in Figure 3 of their article suggest the possibility of a threshold at about 15 ppb. Smith et al. (2009) fit a piece-wise linear ambient concentration-response relationship with breaks at 40 ppb and 60 ppb, and they showed that a hockey-stick-shaped curve with a break at 40 ppb is consistent with the data. Our analysis shows evidence of a threshold at a little over 30 ppb ( Figure 1B). However, the data are too sparse to draw any firm conclusions regarding the shape of the concentrationresponse curve at low concentrations. We emphasize here that our analy ses are based on 24-hr average concentrations of O 3 , whereas regulation is based on the maximum 8-hr averages. Pollutant effect estimates cannot easily be converted from one averaging time to another. Smith et al. (2009) showed that a simple scaling procedure is inadequate. It is, therefore, important to repeat these analy ses with other measures of exposure.
CO, NO 2 , SO 2 . The only previous multicity analy ses of these three pollu tants were conducted during re analy ses of NMMAPS (Dominici et al. 2003). For increases of 1 ppm for CO and 10 ppb for NO 2 and SO 2 , respectively, the NMMAPS investigators reported statistically significant increases in all-cause mortality at a 1-day lag of approximately 0.5%, 0.25%, and 0.6%. Our single-pollu tant estimates of the mean are higher for each of these pollu tants (e.g., 1.3%, 1.0%, and 1.5%, respectively, using 100 df splines for time trends). In a multi pollu tant model including these three gases, CO was no longer significantly associated with mortality (mean 0.64%; 95% CI: -0.20, 1.09), but NO 2 and SO 2 remained so (mean 0.62; 95% CI: 0.40, 0.98 and mean 0.82; 95% CI: 0.48, 1.15, respectively). In interpreting this result it should be kept in mind that, because of missing data, multi pollutant analy ses are based on a smaller set of cities than single-pollu tant analy ses. Concentrationresponse relationships for CO and SO 2 are consistent with linearity. Although there is the suggestion of a threshold at around 20 ppb for NO 2 , as for PM 10 , the confidence bounds are too wide to rule out linearity.

Conclusions
We used the subsampling bootstrap procedure to derive maximum likelihood estimators of national effects of exposures to criteria pollutants on deaths the following day. The AIC for the fitted models provided little evidence of a common effect estimate across the United States, however (data not shown). While the focus of much air pollution research for the past couple of decades has been on PM and O 3 , we find stronger associations between CO, NO 2 , and SO 2 and mortality. The ambient concentration-response relationship for O 3 shows evidence of non linearity. Regional and seasonal analy ses using the methods described in this study may offer further insight. Previous publications have reported that the results of time-series analy ses of air pollution data in individual cities can be highly sensitive to choice of statistical model (e.g., Clyde 2000; Koop and Tole 2004;Moolgavkar 2003). While our analy ses are based on the most recent data availa ble to us (which are not identical to the data used in previous analyses), our results suggest that different statistical approaches to multi city analy ses can yield disparate results.