The Distributed Lag between Air Pollution and Daily Deaths

Abstract

Many studies have reported associations between air pollution and daily deaths. Those studies have not consistently specified the lag between exposure and response, although most have found associations that persisted for more than 1 day. A systematic approach to specifying the lag association would allow better comparison across sites and give insight into the nature of the relation. To examine this question, I fit unconstrained and constrained distributed lag relations to the association between daily deaths of persons 65 years of age and older with PM₁₀ in 10 U.S. cities (New Haven, Birmingham, Pittsburgh, Canton, Detroit, Chicago, Minneapolis, Colorado Springs, Spokane, and Seattle) that had daily monitoring for PM₁₀. After control for temperature, humidity, barometric pressure, day of the week, and seasonal patterns, I found evidence in each city that the effect of a single day’s exposure to PM₁₀ was manifested across several days. Averaging over the 10 cities, the overall effect of an increase in exposure of 10 μg/m³ on a single day was a 1.4% increase in deaths (95% confidence intervals (CI) = 1.15–1.68) using a quadratic distributed lag model, and a 1.3% increase (95% CI = 1.04–1.56) using an unconstrained distributed lag model. In contrast, constraining the model to assume the effect all occurs in one day resulted in an estimate of only 0.65% (95% CI = 0.49–0.81), indicating that this constraint leads to a substantial underestimate of effect. Combining the estimated effect at each day’s lag across the 10 cities showed that the effect was spread over several days and did not reach zero until 5 days after the exposure. Given the distribution of sensitivities likely in the general population, this result is biologically plausible. I also found a protective effect of barometric pressure in all 10 locations.

Numerous studies in the last decade have reported associations between day-to-day fluctuations in air pollution and day-to-day fluctuations in daily deaths in cities. ^1–6 These associations have been seen with pollution levels on the same day, or within 5 previous days. ⁷ The studies have been criticized because the same lag relation between air pollution and deaths was not used in all studies. Instead, the best single day relation has usually been fit in each study, leading to inconsistencies in the lag relation specified. A more consistent approach would facilitate combining evidence across studies.

Most of the studies that considered associations with a multi-day moving average ^1,5,6 found that a 2- or 3-day moving average of air pollution has fit better than any single day’s pollution or longer moving average. In some cases longer (up to 5 days) moving averages have been fit. This finding suggests that the effect of an increase in pollution concentration on a single day is distributed across several subsequent days. Toxicological data also suggest that effects of exposure may be seen over several subsequent days. For example, Clarke and co-workers reported changes in tidal volume of rats immediately after exposure to concentrated ambient air particles, and increases in inflammatory markers approximately 36 hours after exposure. ⁸

In general, we might suppose the shape of the distribution of effect of air pollution over time looks something like Figure 1. If the time scale in this figure is short (hours, for example), then a single 24-hour average pollution concentration will be an acceptable proxy for the true relation shown in the figure. If the time scale extends over several days, a 2-day moving average will probably be a reasonable proxy, but a single day’s pollution will be a poorer surrogate. In both those examples, we are approximating the true shape of the relation by a simple step function. If the relation extends over a longer time scale, then the extent to which a 1- or 2-day average concentration will represent the true relation depends on the serial correlation in the air pollution measure. Air pollution concentration measurements taken within a few days of each other tend to be correlated. That correlation varies from city to city, however. Hence some of the variation in the effect sizes and best fitting lags of the reported associations among cities may be due to fitting a simplistic model for the distribution of the exposure-response over time.

FIGURE 1: A hypothesized curve showing the impact of an environmental toxin over time. The effect rises, and then falls, possibly with a long tail. The goal of this analysis is to determine what the actual shape of the curve representing the time course of deaths after exposure to PM₁₀ is.

This issue is increasing in importance. Previous studies of the association between daily pollution and daily deaths have been done in cities with daily monitoring. Recent studies have begun to examine large numbers of cities in which airborne particles are monitored for only 1 of 6 days. ⁹ Two- or 3-day averages cannot be used in those cities, and it is important to quantify the impact of this limitation.

A systematic approach to investigating the distribution of effect over time offers the possibility of explaining some of that variation noted above, of testing the potential bias associated with the use of a single day’s pollution, of providing estimates in individual locations that can more appropriately be combined, and of indicating what the nature of the lag structure between air pollution and daily deaths is. This study demonstrates the methodology for such an approach, and applies it to a study of PM₁₀ pollution and daily deaths in 10 cities across the United States.

DATA AND METHODS

To analyze effectively the distributed lag between PM₁₀ and daily deaths, we need daily PM₁₀ measurements. Most U.S. cities only measured PM₁₀ for 1 of 6 days, but a number of locations had monitors on a daily schedule. I selected 10 U.S. cities with roughly daily PM₁₀ monitoring to provide a reasonable number of locations for a combined analysis. The cities were New Haven, Pittsburgh, Birmingham, Detroit, Canton, Chicago, Minneapolis-St. Paul, Colorado Springs, Spokane, and Seattle. Daily deaths in the metropolitan county containing each city were extracted from National Center for Health Statistics mortality tapes for the years 1986 through 1993. Deaths caused by external causes (ICD-9 800–999) were excluded. Because a previous study suggested the seasonal pattern of mortality might differ by age, ¹⁰ this study is limited to deaths of persons 65 years of age and older. Minneapolis and St. Paul were combined and treated as one city. Daily weather data were obtained from the nearest airport weather station, and daily concentrations of PM₁₀ were obtained from the U.S. Environmental Protection Agency’s (EPA) AIRS monitoring network.

The assignment of PM₁₀ exposure raised a number of issues. Many of the locations have more than one monitoring location, but typically only one monitor operates on a daily basis, with the others operating every 3rd or 6th day. If the monitors were simply averaged, the daily mean would jump on days when new monitors were included merely because their annual average differs from the monitoring station that operates on a daily basis.

The variance of PM₁₀ measurements also can differ from monitoring location to location. Day-to-day changes in which monitors are in the daily average would also result in changes in the day-to-day variation in the exposure measure that do not represent true changes in exposure, but only changes in the sampling of monitors. To remove these influences, I used the following algorithm. The annual mean was computed for each monitor for each year, and subtracted from the daily values of that monitor. I then standardized these daily deviances from each monitor’s annual average by dividing by the standard deviation for that monitor. The daily standardized deviations for each monitor on each day were averaged, producing a daily averaged standardized deviation. I multiplied this by the standard deviation of all of the centered monitor readings for the entire year, and added back in the annual average of all the monitors.

Analytical Approach

For each city, a generalized additive Poisson regression was fit ^11,12 modeling the logarithm of the expected value of daily deaths as a sum of smooth functions of the predictor variables. The generalized additive model allows regressions to include nonparametric smooth functions to model the potential nonlinear dependence of daily admissions on weather and season. It assumes that: log(E(Y)) = β₀ + S₁(X₁)+..+S_p(X_p)

where Y is the daily count of deaths, E (Y) is the expected value of that count, the X_i are the covariates and the S_i are the smooth (ie, continuously differentiable) functions. For the S_i, I used loess, ¹³ a moving regression smoother. This approach is now standard in air pollution time series. ¹⁴ For each covariate, it is necessary to choose a smoothing parameter that determines how smooth the function of that covariate should be. Three sets of predictor variables were used: a smooth function of time to capture seasonal and other long-term trends in the data, weather, and day of the week variables to capture shorter term potential confounding, and PM₁₀. The choice of smoothing parameter for each set of variables is described below.

The purpose of the smooth function of time is to remove the basic long-term pattern from the data. Seasonal patterns can vary greatly between, for example, Birmingham and Spokane, and a separate smoothing parameter was chosen in each city to reduce the residuals of the regression to “white noise”¹⁵ (ie, remove serial correlation). This approach was used because each death is an independent event, and autocorrelation in residuals indicates there are omitted, time-dependent covariates whose variation may confound air pollution. If the autocorrelation is removed, remaining variation in omitted covariates has no systematic temporal pattern, and hence confounding is less likely. Sometimes it was necessary to incorporate autoregressive terms ¹⁶ to eliminate serial correlation from the residuals.

The other covariates were temperature, relative humidity, and barometric pressure on the same day, the previous day’s temperature, and day of the week. To allow for city-specific differences, the smoothing parameters for these covariates were also optimized separately in each location. The criterion used was to choose the parameter for each variable that minimized Akaike’s Information Criteria. ¹⁷

PM₁₀ was treated as having a linear association in this analysis to facilitate the combination of coefficients across cities and examination of lag structure. Robust regression was used to reduce sensitivity to outliers in the dependent variable. These regressions were done using the generalized additive model function in Splus, and M-estimation was the robust regression method. To reduce sensitivity to outliers in the pollution variable, analysis was restricted to days when PM₁₀ levels were below 150 μg/m,³ the currently enforced ambient standard. This approach also ensures that the results are unambiguously relevant to questions of revision of those standards.

Distributed Lag Models

Distributed lag models have been used for decades in the social sciences ¹⁸ and Pope and Schwartz ¹⁹ recently described the use of this approach in epidemiology.

The motivation for the distributed lag model is the realization that air pollution can affect not merely deaths occurring on the same day, but on several subsequent days. The converse is therefore also true: Deaths today will depend on the “same day” effect of today’s pollution levels, the “1-day lag” effects of yesterday’s PM₁₀ concentrations, etc. Therefore, suppressing covariates and assuming Gaussian data for the moment, the unconstrained distributed lag model assumes: where X_t-q is the PM₁₀ concentration q days before the deaths. The overall effect of a unit increase in air pollution on a single day is its impact on that day plus its impact on subsequent days. That is, it is the sum of β₀ + ..+β_q. ¹² To see this more easily, note that Equation 1 can be recast as: Y_t = α + β*(w₀X_t +..+w_qX_t−q) + ε_t

where the wi are weights that sum to one, and β* is β₀ +..+ β_q. That is, β* is also interpretable as the marginal effect of a unit increase in a weighted average pollution variable. Since a unit increase in pollution on a single day increases the weighted average on all q subsequent days, the effect of that single day’s increase will be β*wi on each of the q subsequent days, or β* overall.

Since there is substantial correlation between air pollution concentrations on days close together, the above regression will have a high degree of collinearity. This collinearity will result in unstable estimates of the individual β_q’s and hence poor estimates of the shape of Figure 1. Nevertheless, the sum of the individual β_q’s will be an unbiased estimate of the overall effect of a unit increase in pollution, albeit an inefficient one.

To gain more efficiency and more insight into the shape of the distributed effect of air pollution over time, it is useful to constrain the β_q’s. If this is done flexibly, substantial gains in reducing the noise of the unconstrained distributed lag model can be obtained, with minimal bias. This approach has been widely applied in the social sciences, using Gaussian data. The most common approach is to constrain the shape of the variation of the β_q’s with lag number (ie, the shape of Figure 1) to fit some polynomial function. That is, the polynomial distributed lag model (PDL(q, d) with q lags and degree d) is the Model 1 above, subject to the restriction:

This approach originated with Almon. ²⁰ Here I extend that model to the Generalized Additive Model case. Assume that: log(E(Y)) = covariates + α₀ + S₁(X₁)+..+S_p(X_p)+β₀Z₀+..+β_qZ_q

where Z₀ is the exposure on the concurrent day, Z₁ on the previous day, etc. If we impose the constraints in Equation 2, and suppress the covariates, we can write this as Log(E(Y)) = covariates + η₀Z₀ + (η₀ + η₁ + η₂ + .. +η_d) Z₁ + (η₀ + 2η₁ + 4η₂ + .. + 2^d η_d) Z₂ + (η₀ +q η₁ +q²η₂ + .. +q^dη_d) Z_q

We can rewrite this by collecting terms in each of the η’s. This yields: Log(E(Y)) = covariates + η₀(Z₀ +Z₁ + .. +Z_q) + .. + η_d(Z₁+ 2^dZ₂+ .. + q^dZ_q)

Hence, if we define d +1 new variables W_d to be weighted sums of the exposure variable Z and its lags, with W_d=Z₁+2^dZ₂+..+q^dZ_q and W₀=Z₀ + Z₁+..+ Z_q

we can estimate the model Log(E(Y)) = covariates = η₀W₀ + .. + η_dW_d

and the coefficients of the W ’s will be the parameters of the polynomial distributed lag.

Note that the use of a single day’s exposure is also a constrained lag model. In that case, we are fitting Equation 1, with the constraint that β₁=β₂=··=β_q=0. If we are not quite sure the effects of pollution are limited to a single day, these constraints are much more restrictive than those of Equation 2 are, and much more likely to introduce bias into the estimated overall effect.

I have chosen a maximum lag of 5 days before the deaths for the air pollution variable, because previous studies in the literature have shown that lags of more than a few days had little correlation with daily deaths, and because the goal of this analysis is to estimate the short-term effects of air pollution. Because the explanatory power of air pollution on daily deaths is modest, parsimony in the degree of the polynomial is necessary. Therefore I have chosen a second-degree polynomial in this analysis. To test the sensitivity of the conclusions regarding the overall effect of air pollution on daily deaths to the degree of the polynomial, I have also fit the unconstrained distributed lag model in each city. To see if the traditional approaches of using a 1- or 2-day moving average resulted in a downward bias in the estimated effects of PM₁₀, I have also fit those models in each city.

To combine results across cities, I used inverse variance weighted averages. For the distributed lag model, the effect at each day’s lag (and its variance) was estimated from the parameters of the polynomial (ie, η₀, η₁, and η₂) and their covariances, as was the overall effect and its variance.

RESULTS

Table 1 shows the populations, mean daily deaths of people 65 years of age and over, and means of the environmental variables in the 10 study locations. Table 2 shows the correlation between PM₁₀ and the weather variables. The correlation between PM₁₀ and barometric pressure was quite small and mixed in sign. The correlation between PM₁₀ and relative humidity was generally negative and moderately low. The correlation with temperature varied considerably across the locations, ranging from −0.34 in Colorado Springs to 0.45 in Pittsburgh.

Table 1: Characteristics of the Study Cities (1986–1994)

Table 2: Correlation between PM₁₀ and Other Environmental Variables in the 10 Cities

Weather Results

The most notable result using weather variables is the consistent association between barometric pressure and daily deaths. In every city, higher barometric pressure was associated with fewer deaths. In six of the cities, Akaike’s Information Criteria (AIC), which trades off improvement in model fit vs the number of degrees of freedom that produced that improvement, was reduced by including this term in the model. A lower AIC is generally taken to indicate a better fitting model. In three smaller locations (Canton, Birmingham, and New Haven) AIC was lower without including barometric pressure, and in Detroit, AIC was the same for both models. For comparison, temperature, a term traditionally included in models relating air pollution to daily death, also improved (lowered) AIC in 6 out of the 10 locations, and relative humidity only improved model fit in 4 of the 10 locations. The association with barometric pressure was not linear—the protective effect tended to flatten out at high pressure. This effect is illustrated in Figure 2, which shows the results for Minneapolis-St. Paul. To illustrate better the magnitude of the barometric pressure effect, I also fit linear terms for barometric pressure and performed a meta-analysis of the linear coefficients. Days with high barometric pressure had lower deaths. A 0.25-inches of mercury increase in barometric pressure (which is about the mean interquartile range in the 10 cities) was associated with a 1.58% decrease in daily deaths (95% CI = 1.29–1.86 decrease).

FIGURE 2: The nonparametric smooth plot of the Log Relative Risk of Death vs Barometric Pressure (mmHg) in Minneapolis-St. Paul.

PM₁₀ Results

Table 3 shows the estimated effect of a 10 μg/m³ increase in PM₁₀ in each city using the concurrent day’s pollution, the 2-day moving average, the polynomial distributed lag model, and the unconstrained distributed lag model. The effect size estimates using the unconstrained distributed lag were similar to those estimated using the polynomial distributed lag model, suggesting the constraint introduced little bias. Both of the distributed lag models had substantially greater overall impacts than models using only a single day’s exposure, and moderately larger effects than the 2-day average models.

Table 3: Estimated Percent Increase in Daily Deaths (Standard Errors) for a 10-μg/m³ Increase in PM₁₀ under Different Distributed Lag Models

The distributed lag model explained some of the heterogeneity among cities in their effect estimates. For example, using the 2-day average model, the variance of the estimated coefficients of PM₁₀ across the 10 cities was 3.41 × 10⁻⁷. The average within-city variance of the coefficients was 2.39 × 10⁻⁷, suggesting a heterogeneity not attributable to sampling variability of 1.02 × 10⁻⁷. The across city variance of the estimated overall effect from the polynomial distributed lag model was 3.13 × 10⁻⁷. This value suggests that failure to account for the distributed lag properly accounted for about a quarter of the unexplained variation in effect size among the cities.

Figure 3 shows the combined estimate of the distributed lag between air pollution and daily deaths. It clearly remains positive for several days, before falling toward zero, which explains why the single day’s pollutant model does not do a good job as a proxy for the overall effect of air pollution.

FIGURE 3: The distribution over time of the increase in Daily Deaths (shown as percent increase) associated with an increase of 10 μg/m³ in PM₁₀ on a single day. The curve comes from the polynomial distributed lag model, and shows that the effects are distributed over multiple days, falling to zero at about 4 to 5 days after exposure.

DISCUSSION

In every city in this study evidence was seen that the effect of an incremental increase in particulate air pollution on a given day was spread over several succeeding days. First, the unconstrained distributed lag models in each location always showed greater total effects than the concurrent day models. In addition, Figure 2, which combines data across cities, shows the effects of pollution were spread over multiple days, and do not reach zero until a lag of 5 days. In plots of the distributed lag in each city (not shown), the effect is spread over multiple days in each of the 10 locations.

These results are biologically plausible. Given a distribution of sensitivity to air pollution in the general population, which seems likely, and a distribution of severity of pre-existing illness, one would expect some variation in the time between exposure and response. This expectation seems even more likely for a summary measure like all-cause mortality, which mixes deaths from myocardial infarctions, which have been shown to be acutely triggered by immediate exposure to stress or certain activities, ²¹ with deaths from respiratory disease, in which it may take more time for an exacerbatory event to result in the cascade of biological responses leading to death.

In a recent study of air pollution and daily deaths in Milan, Italy, causes of death were examined separately. ²² Some causes were more strongly associated with particle exposure on the concurrent day, but other causes were more strongly associated with exposure 2 days before. The mixing of such deaths in all cause mortality would naturally result in a distributed lag between exposure and response.

Support is also provided by the London smog episode of 1952. ²³ Air pollution concentrations shot up on December 5, 1952, and there was an immediate, same day, increase in deaths. The curve of increase and decrease of daily deaths in general, however, lagged behind the increase and decrease of air pollution, with peak deaths coming a day or two after peak exposure. This finding suggests that there were substantial lagged effects.

There is also toxicological support for effects that persist longer than 1 day. In addition to the report of Clarke et al cited earlier, ⁸ Lay and co-workers have reported that particles instilled in the lung induced an inflammation that took up to 4 days postexposure to resolve. ²⁴ The implications of this result are also evident in Table 3. The estimated effect of a 10-μg/m³ increase in PM₁₀ on 1 day using the 1-day moving average model was only about 40% of the estimate using an unconstrained or constrained model that estimates effects at longer lags. Notably, the unconstrained model makes no assumption about the shape of the distributed lag, or even the existence of effects at lags greater than 1 day. It merely allows for the possibility that they exist, and estimates them based on the data. Hence, these results indicate that studies that rely on single-day exposure averages will on average substantially underestimate the effect of particle exposure. It is possible that this underestimation may not apply for other pollutants, but it would not be prudent to make that assumption, when the methodology to test it is straightforward.

A 2-day moving average did substantially better than the 1-day average in estimating the effects of PM₁₀, but still underestimated the overall impacts by about 40%. For other outcomes, such as cause-specific mortality or hospital admissions, the distribution of effect over time may have a different pattern, and separate evaluations, using distributed lag models, will be needed to assess the adequacy of simple single- or multiple-day averages. For multi-time series studies, the unconstrained lag model avoids any risk of bias. For single time series, the polynomial distributed lag model appears to risk little bias, and should be the method of choice. The advantage of these distributed lag models is that we do not need to leave the question of how the effects are distributed over time to hazard. By fitting a model that allows but does not require the effect of pollution to be distributed over several days, we can make that question part of our investigation. By using the simple transformation shown in the Data and Methods section, this approach can be implemented in any Poisson regression package. Hence, polynomial distributed lag models should become standard practice in air pollution epidemiology, unless there is clear biological reason for assuming that the response is limited to 1 day. These models can equally well be applied to other acute triggers. For example, aeroallergen exposure and acute asthmatic response, or triggers of acute myocardial infarction may represent areas in which these models can be usefully applied.

The finding that barometric pressure is consistently associated with lower mortality is also of considerable interest. Previously, we have shown that barometric pressure was associated with oxygen saturation. ²⁵ These changes might plausibly influence mortality risk; however, they were noted in a high altitude location (Provo-Orem, Utah). The findings in this study of an association even in cities closer to sea level, such as Minneapolis (Figure 2), suggest further attention should be paid to this variable.