A National Multicity Analysis of the Causal Effect of Local Pollution, [Formula: see text], and [Formula: see text] on Mortality.

BACKGROUND: Studies have long associated [Formula: see text] with daily mortality, but few applied causal-modeling methods, or at low exposures. Short-term exposure to [Formula: see text], a marker of local traffic, has also been associated with mortality but is less studied. We previously found a causal effect between local air pollution and mortality in Boston.
OBJECTIVES: We aimed to estimate the causal effects of local pollution, [Formula: see text], and [Formula: see text] on mortality in 135 U.S. cities.
METHODS: We used three methods which, under different assumptions, provide causal marginal estimates of effect: a marginal structural model, an instrumental variable analysis, and a negative exposure control. The instrumental approach used planetary boundary layer, wind speed, and air pressure as instruments for concentrations of local pollutants; the marginal structural model separated the effects of [Formula: see text] from the effects of [Formula: see text], and the negative exposure control provided protection against unmeasured confounders.
RESULTS: In 7.3 million deaths, the instrumental approach estimated that mortality increased 1.5% [95% confidence interval (CI): 1.1%, 2.0%] per [Formula: see text] increase in local pollution indexed as [Formula: see text]. The negative control exposure was not associated with mortality. Restricting our analysis to days with [Formula: see text] below [Formula: see text], we found a 1.70% (95% CI 1.11%, 2.29%) increase. With marginal structural models, we found positive significant increases in deaths with both [Formula: see text] and [Formula: see text]. On days with [Formula: see text] below [Formula: see text], we found a 0.83% (95% CI 0.39%, 1.27%) increase. Including negative exposure controls changed estimates minimally.
CONCLUSIONS: Causal-modeling techniques, each subject to different assumptions, demonstrated causal effects of locally generated pollutants on daily deaths with effects at concentrations below the current EPA daily [Formula: see text] standard. https://doi.org/10.1289/EHP2732.

Introduction

Hundreds of studies have reported associations between short-term exposure to air pollution and daily deaths (Baccini et al. 2006; Bell et al. 2004; Braga et al. 2001; Chen et al. 2014; Jhun et al. 2014; Katsouyanni et al. 1997; Katsouyanni et al. 2009; Levy et al. 2012; Peng et al. 2005; Peng et al. 2013; Samet et al. 2000; Schwartz 1991; Tao et al. 2012; Zanobetti et al. 2002; Zanobetti and Schwartz 2008, 2009). The most common findings are that associations with particulate air pollution and ozone exist and that these two exposures do not confound each other. Many toxicology and controlled human-exposure studies showing associations of these pollutants with changes in intermediary outcomes (blood pressure, inflammation, autonomic function, endothelial function, thrombosis, etc.) support those findings. (Bartoli et al. 2009b; Calderón-Garcidueñas et al. 2008a; Fakhri et al. 2009; Langrish et al. 2009; Lundbäck et al. 2009; Matsumoto et al. 2010; O'Toole et al. 2010; Peretz et al. 2008).

More recently, studies have reported associations with nitrogen dioxide () with daily deaths, often while controlling for other pollutants (Mills et al. 2015). is a byproduct of combustion and local traffic, particularly Diesel traffic, is a major source. This information has raised questions as to whether the findings represent health effects of itself, or if it acts as a surrogate for some other pollutant from traffic. To date, there has been less toxicological investigation as to how might influence the processes that rapidly generate respiratory and cardiovascular deaths. However, clearly deserves more attention than it has received.

Few time-series studies of acute effects of air pollution applied modern causal-modeling techniques. Causal modeling seeks to analyze observational data in a way that simulates conducting a randomized experiment. Randomization makes exposure independent of all potential confounders, and causal methods seek to replicate that situation, rather than to control for the confounders in the outcome regressions, as conventional analysis does. Under specified assumptions, including ones that are untestable in the data and rely on external knowledge, causal methods yield causal estimates of the effects of exposure. Often, they provide marginal estimates of the effects of exposure, that is, ones that are not conditional on the distribution of covariates and are therefore more generalizable.

Marginal structural models are the best known causal models in epidemiology, estimating the marginal effects of exposure by using inverse probability weights of time-varying exposures to render the exposure independent of the measured covariates. If the exposure is independent of covariates, its effect on the outcome cannot be confounded by them and resulting estimates do not depend on the distributions of confounders. If all important covariates are measured, these models provide causal estimates of the marginal effects of exposure.

Recently, we used an instrumental variable analysis to estimate the causal effects of locally generated air pollution in Boston (Schwartz et al. 2016). The analysis used an instrument for the part of the daily fluctuations in air pollution caused by changes in the mixing height and wind speed, which modify the build-up of locally generated pollution but do not have other plausible connections to daily changes in mortality except through air pollution. If that assumption is true, then the instrument represents variations in local pollutants that are randomized with respect to confounders, measured or unmeasured, and therefore provides a causal estimate of the effect of local air pollution concentrations. However, a lower mixing height increases the concentration of all locally emitted pollutants. Thus, these models do not provide much guidance on the relative importance of those pollutants. They do, however, control for unmeasured confounders and are thus complementary to the marginal structural model, which can estimate independent effects of and , but rely on the assumption of no unmeasured confounders to provide unbiased estimates of effect.

Negative exposure controls are also used in causal modeling. A negative exposure control identifies a negative exposure variable, which is likely to be correlated with unmeasured potential confounders but could not be a cause of the outcome of interest. For example, a negative control can be the exposure of interest after the outcome has occurred. In this case, negative exposure controls serve as instruments for the unmeasured confounders. If such confounders exist, control for the negative exposure would be expected to reduce or eliminate the estimated effect of the exposure of interest.

Here we expand the instrumental variable approach we used in Boston to 135 cities across the United States to gain a more robust understanding of the causal effects of local pollution, supplemented with a negative exposure control analysis and extended by adding marginal structural models to estimate separate effects of and . We also implement the negative exposure control in the marginal structural model to provide greater assurance of the assumption of no important unmeasured confounders.

Because data were available in fewer cities than , we fit the instrumental variable analysis in the larger number of cities and calibrated the instrument to . Then, we fit marginal structural models for both pollutants.

Data and Methods

Data

Mortality and atmospheric data.

We obtained data for 135 cities between 1999 and 2010 on daily deaths, , and weather variables, and with wide geographic spread in the United States. Of these, 105 had measurements. Figure 1 shows a map of the city locations. Daily deaths from any natural cause (ICD-10: V01-Y98, ICD-9: 1–799) of persons who resided in the city where they died were obtained from the National Center for Health Statistics (NCHS), from 1999 through 2006 after first obtaining permission from each state health department. After 2006, the National Center for Health Statistics stopped providing date of death, so we obtained data from public health departments in individual states.

Figure 1.

A map of the United States with the cities included in the analysis. The cities with a triangle are cities with both and .

We obtained and concentration data from the U.S. EPA Air Quality System Technology Transfer Network (Air Quality System Technology Transfer Network 2012).

Many cities have more than one monitor for or , and values at individual monitors can be missing on some days. Hence, daily means can change on days when a particular monitor was present or absent simply because its values tended to be higher or lower than average. This circumstance would not represent true changes in exposure. Therefore, we used a standard algorithm that adjusts for differences in means and standard deviations among monitors (Schwartz 2000). Daily mean temperature, wind speed, and sea-level pressure data in every city were obtained from the National Oceanic and Atmospheric Administration (NOAA) (National Centers for Environmental Information (NCEI) 2012). Height of the planetary boundary layer data was obtained from the NOAA reanalysis dataset (NOAA 2010).

Statistical Analysis

The analysis was conducted on a city-specific level and then aggregated across cities using a random effects meta-analysis. All analyses were conducted with R (version 3.3.1; R Core Team).

Causal Modeling

Causal modeling contrasts the results of two potential outcomes: what would have been observed had the entire population been exposed to exposure a, vs. observations made had they all been exposed to a'. At most, one potential outcome is observed, and various methods provide legitimate surrogates for the unobserved potential outcome under certain assumptions, some of which are untestable (Hernán et al. 2008). Because those assumptions differ between methods, using several different methods provides additional assurance for the robustness of the results. In this paper, we apply three approaches: an instrumental variable, marginal structural models, and negative exposure controls.

I. Instrumental variable with planetary boundary layer, wind speed, and sea level pressure as instruments.

Let be the potential outcome (daily deaths) in the population of a city exposed to on day t, and let be the potential outcome under an alternative exposure a′. We assume the potential outcome depends on predictors as follows:

Where represents the potential outcome at time t under exposure a, and are the intercept and the slope of exposure, respectively, and represents all the other predictors of outcome. E denotes expected value and log is the natural logarithm.

Suppose there is a variable Z that is a source of variation in exposure, and Z is associated with Y only through A. Z is called an instrumental variable. We can then express as follows: where, represents the other sources of variation in exposure, and in particular, all the exposure variations that are associated with other measured or unmeasured predictors of outcome which are included in . Then let Z1 and Z2 be the values of Z such that: and and therefore,

As a result, if we use Z as an instrument for A, we obtain a causal estimate for , which is the log rate ratio. This estimate is true even if there are unmeasured confounders, provided they are not correlated with Z. This scenario is displayed in the Directed Acyclic Graph in Figure 2A.

Figure 2.

Directed Acyclic Graphs (DAGs) showing the causal paths in A) the Instrumental Variable Analysis, B) the Marginal Structural Models, and C) the Negative Exposure Control models. Solid lines with arrowheads indicate directed causal paths. Conditioning is indicated by a box around the variable. In A), the association between A and Y induced by C is blocked by conditioning on C; the association between A and Y induced by U is open. However, the Instrument I is independent of U, and allows estimation of a causal path to Y. In B), A is independent of C after inverse probability of exposure weighting, indicated by the dotted line. In C), conditioning on descendent blocks the association between and through under the assumption that is the sole descendent of , and partially blocks it otherwise.

Creation of instrument.

The air above a city contains both locally emitted and transported pollutants. The lower atmosphere has substantial vertical mixing that dilutes local emissions into the air above. This mixing mostly disappears above a certain elevation, called the height of the planetary boundary layer (PBLH) (Finlayson-Pitts and Pitts 1986). Therefore, the impact of local emissions on pollutant concentrations varies inversely with the PBLH (i.e., for the same local emissions, concentrations are higher when the boundary layer is low and vice versa) (Seinfeld and Pandis 1998). The mean PBLH varies substantially from day to day, and PBLH is mainly determined by the thermal convection from the earth’s surface (rising during the day and sinking at night); by the balance between turbulent kinetic energy production (from wind flow) and its dissipation; by the vertical wind shear; by the density stratification of the atmosphere above the boundary layer; and by the water content of the air from evaporation and blown in moisture. These processes are unlikely to be associated with daily deaths except through air pollution changes, and therefore, PBLH is an attractive option for an instrumental variable.

Besides vertical transport, locally emitted air pollutants are also transported horizontally, where the influence of local sources increases with decreasing wind speed and vice versa. Higher wind speed also produces more turbulent mixing of the pollutants into the air above the surface (but below the PBLH), thereby reducing concentrations. Outside of extreme events, wind speed is an unlikely predictor of health other than through air pollution, making it another potential instrumental variable.

Finally, high atmospheric pressure generally brings weather conditions, such as lower vertical temperature gradients, that inhibit vertical and horizontal mixing of pollutants. Again, there is no obvious direct connection between atmospheric pressure and daily deaths. As such, PBLH, wind speed, and atmospheric pressure represent attractive options as instruments for local pollution. Although correlated, each may capture some variation of air pollution that is missed by the others, so constructing an instrument by combining the three should improve power and avoid the problems of a weak instrument.

However, PBLH, atmospheric pressure, and wind speed may vary seasonally and with temperature, and have independent associations with mortality. Within strata of month-by-year and day of the week, and after control for temperature, further association with predictors of health is unlikely. Hence, we looked only at pollution variation within month-by-year strata and with control for temperature and day of the week, and then calibrated that variation with our instruments—that is, we assume that short-term predictors of mortality such as smoking, psychological stress, etc., are uncorrelated with PBL height on a day-to-day level, within strata of month-by-year, and after control for day of the week and nonlinear control for temperature. Specifically, first we fit the following model in each city: where s(x) denotes a penalized adaptive spline for the variable x. The residuals (e) from that model are independent of time trend, season, and temperature.

There may be days when no monitored values are present. One reason we chose indicator variables for each month of each year as our control for time trends and season, rather than natural splines of date, is that this approach does not suffer the problem that missingness creates for splines. The knots of natural splines are located by placing them equally across the nonmissing days. If half the days in January are missing, and we use 4 df per year in our spline, instead of placing a knot at the end of March, it will be placed in mid-April, possibly undercontrolling for season. With indicator variables for each month of each year, as long as some days are present in a month, we can control for it, and this method does not influence the number of degrees of freedom used to control for other months.

To produce a single pollution-calibrated instrumental variable, we combined information from PBLH, atmospheric pressure, and wind speed on the day of death (lag 0) and the day before death (lag 1). Because their effect on particle concentrations need not be linear and may interact, we used a support vector regression (SVM) (Cortes and Vapnik 1995) with a radial kernel to estimate the variation of the residuals that was explained by these variables. We used the SVM function in the R package e1071. Consistent with most previous literature, we used the mean of on the day of death and the day preceding death in deriving the instrument. By calibrating to the residuals of . This method allows the coefficient of the instrument to be interpreted on the same scale as .

Using the instrument as exposure, we fit a quasi-Poisson regression in each city predicting all-cause mortality. We stratified by each month of each year using indicator variables, and estimated the rate ratio for the instrument. The control for month-by-year in the model is to correct for the substantial overdispersion of counts of daily deaths due to seasonal variation and time trends in mortality rates.

We tested whether the instrument was independent of the covariates by calculating the correlation coefficient with temperature and by comparing boxplots of the instrument by year of study and month of study.

Effects below the National Ambient Air Quality Standard (NAAQS).

To see if the association remained when our analysis was restricted to days below the NAAQS, we limited our sample to days that when was at least below the current standard (that is, below ). This limitation ensures that even with measurement error, the association cannot be influenced by exposures above the NAAQS. The instrumental variable was constructed on the restricted data set, and the analysis of its association with daily mortality was repeated.

II. Marginal structural models.

Marginal structural models assume that there are no unmeasured confounders. With this assumption, inverse probability of exposure weighting (IPW) can provide causal estimates of the marginal effect of an exposure because the weighting renders the exposure independent of the measured confounders. This scenario is shown in the Directed Acyclic Graph in Figure 2B, where the dashed line between exposure and confounder indicates that the weighting has removed the association between the two. Those weights are derived from propensity scores.

The propensity score is the probability, given the confounders, that an individual would have received the exposure that they got. For a dichotomous exposure, it can be estimated using a logistic regression predicting exposure category as a function of the confounders. In the case of a continuous exposure, we estimate the probability density of receiving that exposure, given the covariates, by fitting a linear regression predicting the exposure as a function of the covariates. The predicted value of exposure from this regression is the expected exposure given the confounders, and the measured value is the actual exposure. The difference is the residual, and the probability density of the residuals gives us the probability density of receiving the actual exposure, given the confounders (Imai and van Dyke 2004). This is generally normalized by the marginal density of the exposure.

Propensity score models may be sensitive to correctly specifying the relation between exposure and covariates. Hence, we used a flexible model that estimates city-specific propensity score models predicting exposure or exposure each day as a function of indicator variables for each month of each year, day-of-the-week, and penalized splines for temperature, previous day’s temperature, and, for each pollutant, the other pollutant. To ensure adequate capture of the association with continuous covariates, we used adaptive splines. For , in each city we fit the model: where PM2501 is the mean of on day lag 0 and lag1, and NO201 is the two-day mean for . From these models, we extracted the residuals and their standard deviation (SD), and computed the probability density of each observation (and, hence, each residual). The inverse probability weight is the inverse of this density, stabilized the marginal probability of in that city.

The computation of the IPW for was identical, but with an adaptive spline for in the propensity score model.

Generalized propensity scores are subject to possible outlier values. To avoid this, these scores were truncated at the 2.fifth and 97.fifth percentile of their distribution in each city (Cole and Hernán 2008).

Positivity.

Positivity is a major concern in such models. That is, unbiased estimates rely on the assumption that any unit of observation (in our case the day), regardless of its covariate values, has some positive probability of obtaining any exposure. To address this aspect, we computed, for all observations, the probability (given the covariates) that they could have had an exposure at or below the 2.5th percentile in that city, or at or above the 97.5th percentile of the exposure distribution, and then excluded all days whose probability was . We fit Poisson regressions for each city in these trimmed samples and combined across cities with a random-effects meta-analysis. To examine if positivity was an issue before the exclusion, we computed the average probability that an exposure day that was at or above 97.5th percentile in that city could have been (given the covariates) below the 2.5th percentile in that city, and vice versa, for each pollutant in each city.

We implemented a negative exposure control in the marginal structural model by including the mean pollution two and three days after the deaths in the propensity score model for each pollutant and examined what change in the coefficient of exposure occurred.

Effects below the NAAQS.

To assess effects below the NAAQS, we restricted the marginal structural model to days when was below . The propensity score model was refit, the weights were recomputed, and the truncation of weights and positivity subsetting were all repeated on this restricted data, and the marginal structural models refit.

III. Negative exposure control.

The concept of a negative control exposure is straightforward. Suppose there is an unmeasured confounder. Because we are conducting a time-series analysis, the confounder must be correlated with exposure and outcome over time (but within month). For example, perhaps diet fluctuates in correlation with air pollution. It seems reasonable to assume that any such variable that is correlated with pollution the day before a death occurs will also be correlated with pollution the day after the death, unless it tracks pollution extraordinarily tightly. Failure to control for it will bias the association of exposure with outcome. However, it will also induce an association between exposure after the death occurs and the death, which cannot be causal. This situation is illustrated in Figure 2C, Hence, failure to find a positive association between exposure after the event and mortality is evidence that no such confounder exists, and control for exposure after the event should capture much of the confounding effect of the unmeasured covariate. Failure to find much of a change in the association of exposure with outcome after control for future exposure suggests that little confounding is occurring.

We applied the negative exposure control to both the instrumental variable model and to the marginal structural models for and . Specifically, we fit a model with both the instrument on the day of and day before the death, and with the instrument two and three days after the death. We allowed a one-day gap between the two exposures because has a high correlation between two consecutive days. Similarly, we fit marginal structural models with future (and ) two and three days after the death included in the propensity score model for (and ) before the death.

Results

There was a total of 318,109 d with no missing data for in the study, and on those days, there were 7,277,274 deaths. The mean value of during the study was with an interquartile range from 7.5 to (Table 1). The city-specific mean ranged from in Fresno, California, to in Tucson, Arizona, and in Albuquerque, New Mexico [see Table S1 for city-specific means (±SD) for all-cause mortality, , and , and Table S2 for city-specific means for weather variables]. For , 105 cities were included in the analysis, with a mean value of , and an interquartile range from 9.5 to . The city-specific means ranged from in Elizabeth, New Jersey, to in Bath, New York.

Table 1

Descriptive statistics for air pollution, meteorological variables, and daily mortality in 135 U.S. cities, 1999–2010.

Variable	Mean	25% Percentile	75% Percentile
Temperature (°F)	57	45	72
PM_2.5 (μg/m^3)	12.8	7.5	16.1
Wind Speed (m/s)	6.7	4.4	8.6
Sea level Pressure (mmHg)	1,017	1,013	1,021
NO_2 (ppb)*	15.3	9.5	19.8
Daily Deaths (# per day)	22.8	8.0	27.0
Height of PBL (m)	854	572	1,041

Note: All values are daily averages.

was available only in 105 cities.

Instrumental Variable Analysis

Temperature, previous day’s temperature, day of the week, time trend, and season explained an average of 36.7% of the variation in , and this was removed before fitting the instrumental variable. The instrumental variable explained, on average 18.4% of the remaining variation of . We found a significant association with the instrument, with an estimated causal increase in daily mortality of 1.54% per (95% CI: 1.12%, 1.97%) (Table 2). There was some heterogeneity across cities in the estimate (). The correlation of the instrument with temperature was 0.05, and with previous day’s temperature was 0.003, and there was no pattern of change in the instrument by year, or month of year (Figure 3A and 3B).

Table 2

Estimated percentage change in daily mortality with an increase in the mean value of the instrumental variable (), (), or (), respectively, on the day of and day before death (pooled city-specific estimates derived by random effects meta-analysis).

Regression Model	% Change	95% CI
Instrumental Variable (10 μg/m^3)a	1.54%	1.12%, 1.97%
Instrumental Variable (10 μg/m^3) With Negative Controla,b	1.54%	1.12%, 1.97%
Marginal Structural Modelsc
PM_2.5 (10 μg/m^3)	0.75%	0.35%, 1.15%
PM_2.5 with Negative Control (10 μg/m^3)b	0.79%	0.36%, 1.23%
PM_2.5<25 μg/m^3d	0.83%	0.39%, 1.27%
NO_2 (10 ppb)	2.59%	1.78%, 3.40%
NO_2 with Negative Control (10 ppb)b	2.62%	1.81%, 3.43%
Conventional Time Seriese
PM_2.5 (10 μg/m^3)	0.60%	0.34%, 0.85%
NO_2 (10 ppb)	0.38%	0.08%, 0.69%
PM_2.5<25 μg/m^3	0.62%	0.32%, 0.93%

Instrumental Variable models: quasi-Poisson regression models stratified on month-by-year.

Negative Controls: Models with negative controls are adjusted for mean IV, , or , respectively, on the second and third day after death, in addition to the exposure on the day of and day before death.

Marginal Structural Models: Fit with city-specific inverse probability weights based on month, day-of-the-week, temperature, previous day’s temperature, and, for each pollutant, the other pollutant.

: Percentage change in daily mortality with a increase in on the day of and day before death, restricted to days with below the .

Conventional Time Series: Models of or with penalized splines for temperature (same day and previous day) and indicator variables for the month-of-year and day-of-week.

Figure 3.

A boxplot of the value of the instrumental variable for by month of the year (A) and by year (B) in the study, to examine whether exposure is balanced by month and year. The central line in a boxplot is the median, the top and bottom of the box are the upper and lower quartiles, and the circles show the extreme values.

Effects below National Ambient Air Quality Standards

We repeated the instrumental variable analysis, restricting the sample to only days when concentrations were below . This restriction left us with 6,705,626 deaths on 296,096 days. The instrument for was associated with daily deaths, with an estimated causal effect size of 1.70% per (95% CI: 1.11%, 2.29%).

Negative Exposure Control

Using the negative control exposure, there was no change in the estimated effect or confidence interval for the instrument, and the negative control exposure was not associated with mortality ( change in mortality per increase in negative control exposure, 95% CI: , 0.3%).

Marginal Structural Model

Positivity analysis.

For , before positivity exclusions, the probability that a high exposure () could have been low () was not trivial. Across cities, the median probability was 0.48%, and the first and third quartiles were 0.17% and 0.91%, respectively. However, the probability that a low exposure could have been high was lower. Across cities, the median probability was 0.05%, with first and third quartiles of 0.008% and 0.15%.

For , violation of positivity seemed to be a greater concern. Across cities, the mean probability that a high exposure could have been low was only 0.0014% at the first quartile of cities, 0.014% at the median, and 0.05% at the third quartile. The probability that a low exposure could have been high was even lower: 0.0006% at the first quartile, 0.0035% at the median, and 0.018% at the upper quartile. These are low enough to indicate that positivity was a serious concern when looking at and daily deaths. This noncomparability of high and low exposure days even after control for covariates is concerning for standard outcome regressions as well as propensity score methods. Our analyses were done on the trimmed data that excluded days where positivity might be problematic.

Regression results.

For , we found an effect size estimate of 0.75% (95% CI: 0.35%, 1.15%) increase in daily deaths per (Table 2). When we incorporated negative exposure controls in the propensity score model and repeated the analyses, the effect size became 0.79% (95% CI: 0.36%, 1.23%) (Table 2).

For , the effect size estimate was a 2.59% (95% CI: 1.78%, 3.40%) increase in daily deaths per increase (Table 2). After incorporating the negative exposure control into the propensity score model, the effect size estimate became 2.62% (95% CI: 1.81%, 3.43%) (Table 2).

Effects below NAAQS.

We further restricted our analysis to days when was below , which is below the current U.S. EPA 24-h standard. We continued to find a significant association with in this analysis with an effect size of 0.83% per (95% CI: 0.39%, 1.27%) (Table 2). U.S. EPA does not have a 24-h standard.

Comparison with conventional analysis.

We fit standard time-series analyses with and in a model with penalized splines for temperature (lags 0 and 1) and indicator variables for every month of every year and day-of-week. We found similar effects for (0.60% CI: 0.34%, 0.85%), but considerably lower, but still significant estimates for (0.38% CI: 0.08%, 0.69%) (Table 2). Restricted to concentrations below , we still see a highly significant association between and daily deaths in that model, with a slight increase in the effect size estimate (0.62%, CI: 0.32%, 0.93%) (Table 2).

Discussion

We used an instrumental variable analysis, a negative exposure control, and marginal structural models to estimate the acute causal effects of and exposure on daily deaths. Using the instrumental variable, we found a significant association of local air pollution with daily deaths in the 135 cities. If the instrumental assumption is valid, this association provides a causal estimate of the effect of locally generated pollutants on daily mortality. This finding persisted after incorporating a negative exposure control, which should control for any unmeasured confounders that are correlated with local air pollution several days in the future as well as on the day of study. Restricting to days when was less than resulted in a somewhat larger estimated effect size.

Our instrument only explained a bit over 18% of the variation in after removing effects of season, time trend, and temperature. By design, we only want the instrumental variable to capture the changes in pollution associated with changes in boundary layer height, wind speed, and pressure because we believe those variations will be independent of confounders. However, a limitation of an instrumental variable is that if it does not capture enough of the variation in exposure, it will be too weak to capture an association. In our analysis, we have 135 cities and more than 7 million deaths, giving us substantial power. Hence, we view the moderate association of the instrument with as a strength (less likely to capture variations that may be correlated with other predictors of mortality), rather than a weakness.

With marginal structural models, we found independent associations of both and with daily deaths, which, under different assumptions, are also causal estimates. We again found that these effects remained at low concentrations well below the current EPA 24-h standard of . These effects persisted after incorporating negative exposure controls.

A key attribute of instrumental variable analysis is that it provides protection against unmeasured confounders. The large number of cities (135) and deaths (7.3 million) further supports this conclusion, because the relationship between omitted confounders and the instrumental variable, if any, is likely to vary considerably from city to city. Yet the statistic from the meta-analysis indicated little actual heterogeneity in the effect size estimate across cities. We cannot think of a causal factor for number of daily deaths that would be associated with the mixing height of the atmosphere except for season and temperature. For this reason, we removed seasonal variation and temperature associations from before calibrating the instrumental variables to the remaining variation of . We then confirmed that our instrument was not associated with temperature or seasonal patterns.

The negative exposure control analysis provides further assurance that the assumption of no unmeasured confounders is met. If a common cause of air pollution and acute mortality events fluctuated on time scales shorter than a month in correlation with air pollution, we would expect at least some correlation with air pollution two days in the future as well. Controlling for the negative exposure control would therefore be expected to reduce the effect estimate for air pollution, and, via a backdoor path through the omitted confounder, be associated with daily mortality. Neither of these happened, suggesting that there is no such confounder. Together, these two approaches controlling for measured and unmeasured confounding provide considerable evidence that local air pollution is causally associated with daily mortality.

The marginal structural models provide evidence of causal effects of both and with the assumption of no unmeasured confounders. As we have noted previously, time series have few potential confounders because most predictors of mortality do not vary day to day, or if they do, are unlikely to be correlated with air pollution (Schwartz et al. 2015). The major potential confounders (temperature and the other pollutant) were included in the propensity score. Again, the negative exposure control protects against unmeasured confounders. If unmeasured variables are correlated with exposure 2–3 d in the future as well as exposure on the day of the events, then the negative exposure control should partially control for them and change the estimated effect of or . However, we did not find changes in the estimated effect size estimates after including negative exposure in the propensity score.

Many studies have demonstrated associations between daily fluctuations in and daily fluctuations in mortality, including many large multicity studies (Bell et al. 2013; Katsouyanni et al. 2009; Katsouyanni 1997; Peng et al. 2005; Samet et al. 2000; Schwartz and Marcus 1990; Schwartz and Dockery 1992; Zanobetti et al. 2002; Zanobetti and Schwartz 2009; Zanobetti et al. 2014). A recent meta-analysis of time-series studies of reported an overall effect size of 1.04% (95% CI: 0.52%, 1.56%) (Atkinson et al. 2014). Our estimated effect size, controlling for and with positivity exclusions, was 0.75% (95% CI: 0.35%, 1.15%), which is similar to other studies that did not use causal modeling techniques. Given the previous literature on and daily deaths, and the extensive toxicological evidence, we believe our findings further supports that the link between air pollution and mortality is causal.

A recent systematic review and meta-analysis examined the association of short-term exposure with mortality and hospital admissions (Mills et al. 2015). They identified 204 time-series studies that had examined the association of with either daily deaths or hospital admissions. There were 101 studies that examined the association of with daily mortality and found a significant association with all-cause mortality: a increase was associated with a 0.71% increase in daily deaths (95% CI: 0.43%, 1.00%). Converting that estimate to an effect size per increase yields a point estimate of 1.33%, which is smaller than the estimated effect size of 2.59% increase we found for the same increase in . However, it is larger than the estimate of 0.38% we found with a conventional time-series analysis.

One key difference between the causal estimate and the conventional time-series estimate is that the propensity score for included a nonlinear spline for (and vice versa). The conventional time series has them both as linear predictors of daily deaths. If the association between them is not linear, this could result in residual confounding, which may have reduced the effect size estimate in the conventional analysis. Additionally, although many papers have shown essentially linear associations between and daily deaths, this is less clear for . As noted above, the positivity exclusions for were primarily days with high values, which had extremely low probabilities of being low. This exclusion resulted in an analysis more focused on low exposure days, where the slope may be different.

Our instrumental variable analysis found a 1.5% increase in daily deaths per increase. In contrast, the marginal structural model estimate for was 0.75%, for the same increment. Several factors may explain this difference. First, although the instrument was calibrated with , it likely captured variations in other locally generated exposures such as or other traffic pollutants. This finding is consistent with our finding independent effects of in the marginal structural models. Second, the particle variation that is being captured by the instrument is primarily that of elemental and organic carbon particles from local fuel combustion. These particles may be more toxic in comparison with the average particle, which includes those from long-range transport and which was the exposure of the marginal structural model.

Supporting these hypotheses, most previous time series-studies reported smaller coefficients in the United States than the instrumental variable analyses reported. For example, Zanobetti and Schwartz (2009) found an average effect of 0.98% (95% CI: 0.75, 1.22%) for a increment in a study of 112 cities, similar to the 0.75% increase in our marginal structural analysis, but well below the 1.54% (95% CI: 1.12%, 1.97%) from the instrumental variable analysis.

Causal modeling is an approach to estimating marginal effects of changes in exposures that are causal, conditional on some assumptions. The validity of some of the assumptions is not testable and require outside information to justify. In this study, we tested assumptions for our instrumental variable with a negative exposure control. Both approaches can control for measured and unmeasured confounders. We similarly fit marginal structural models with negative exposure controls to produce more robust causal estimates. Nevertheless, causality of an association is a conclusion of humans, not an output of a statistical model. Supporting evidence from other types of studies is needed to form such a conclusion. There is a substantial body of toxicological and controlled human exposure studies that have revealed plausible biological pathways for these associations, especially for at relevant doses.

In human exposure studies, 50 healthy subjects exposed to air from a busy street () vs. filtered air () for 5 h had a 25% reduction in nitroglycerin-induced vasodilation, increased sympathetic tone, and decreased parasympathetic tone at the higher concentration (Hemmingsen et al. 2015). Both exposures were below the current EPA 24-h NAAQS . An intervention trial in Beijing had subjects walk the streets for 2 h twice, once wearing a particle-filtering mask and once without. Blood pressure was measured continuously during the two 2-h walks, and it was lower when subjects wore the filter (Langrish et al. 2009). A randomized trial of air filtration in the elderly has shown improvements in microvascular function following a 48-h exposure to filtered air vs. unfiltered air (Bräuner et al. 2008). A recent randomized trial of air filtration in college students found sham filtration was associated with higher blood pressure, insulin resistance, serum lipids, fasting glucose, cortisol, epinephrine, and norepinephrine (Li et al. 2017). Finally, particle exposure was associated with increased sICAM-1 and sVCAM-1, markers of endothelial activation, in a controlled human exposure chamber study (Brook et al. 2009).

Animal and in vitro studies further support these findings. Dogs exposed to particles vs. filtered air experienced greater ischemia after temporary occlusion of the coronary artery (Bartoli et al. 2009a; Wellenius et al. 2003). Endothelin-1, a potent acute inducer of vascular contraction and blood pressure increase, increased directly in response to urban PM (Bouthillier et al. 1998; Calderón-Garcidueñas et al. 2007; Calderón-Garcidueñas et al. 2008b; Chauhan et al. 2005; Tamagawa et al. 2008; Thomson et al. 2004, 2007). Black soot particles from diesel exhaust increased oxidative stress in endothelial tissue and induced the production of heme oxygenase-1, a rapid-response defense against oxidative stress (Furuyama et al. 2006; Hirano et al. 2003). Finally, when animals breathing Boston air were placed in a chamber with filtered air, the concentrations of reactive oxygen species in the heart and lung fell by a third within days of removing the exposure (Evelson and González-Flecha 2000).

Toxicological and controlled exposure studies supporting short-term exposures to producing increases in daily deaths are much less common. Most such studies have examined asthma or asthma-related end points that are not relevant for mortality. Studies of exacerbation of COPD have been mixed with some reporting no association (Gong et al. 2005; Morrow et al. 1992), and others reporting small increases, but only at very high exposures not seen in ambient air (Linn et al. 1985; Vagaggini et al. 1996). In contrast, there is substantial toxicological support for increasing susceptibility and response to respiratory infections (Parker et al. 1989). However, these studies were mostly at concentrations of or higher.

With respect to cardiovascular outcomes, controlled human exposure studies have not found increases in blood pressure or changes in cardiac output following exposure (Folinsbee et al. 1978; Gong et al. 2005; Huang et al. 2012). Healthy adults performing intermittent exercise following exposure to (Scaife et al. 2012) were evaluated for electrocardiogram changes. No significant changes were seen, except for a marginally significant decrease in QT interval. Inflammatory markers in the blood were not significantly affected.

The lack of supporting toxicological evidence for the strong associations reported here and elsewhere between and daily deaths suggests that the association may represent the effect of other components of traffic exhaust for which is acting as a surrogate. Diesel particles, ultrafine particles, and polycyclic aromatic hydrocarbons are plausible hypotheses for what those other components might be. However, more toxicology is needed to address this question.

Conclusions

We applied a variety of causal modeling techniques to estimate associations between locally generated air pollution and daily deaths in up to 135 cities across the United States, encompassing over 7 million deaths. Each of these approaches provides a causal estimate subject to assumptions, and these required assumptions differed by the technique. The consistency of results across the different sets of assumptions supports a judgment that the association with is causal. This assertion is further supported by extensive existing toxicological literature. In contrast, although we find clear evidence for an association with that is independent of , the lack of toxicological support makes a causal conclusion more difficult, and this association may represent representing other, more toxic components of traffic exhaust. The associations continue when restricted solely to days well below the current EPA standard for , indicating that standard is not protective of public health.

Click here for additional data file.

Click here for additional data file.