Air pollution, especially airborne particles, has been consistently reported to be associated with daily deaths in reports from all over the world. ^{1–8} More recently, systematic multicity analyses have confirmed these findings. ^{9–12} Nevertheless, some have questioned the public health significance of these associations, arguing that if these deaths are occurring only in those who would have died in a few days anyway, the public health significance of exposure is small. Were that the case, the increase in deaths during and immediately after exposure would be counterbalanced by a deficit in daily deaths a few days later, when those deaths would have otherwise occurred. If such a pattern were true, the positive correlation seen between daily deaths and exposure shortly before the death would be counterbalanced by a negative correlation between exposure and daily deaths at some longer lag. An example of such a hypothetical pattern, called mortality displacement or harvesting effect, is seen in Figure 1. Were such a phenomenon to exist, it should be detected readily in studies of acute episodes, but those patterns have not been observed in air pollution episodes. ^{13}

It is useful to examine the reason for such a phenomenon. Assume there is a pool of people at high risk of dying at any given time. An air pollution episode, by increasing the risk in that pool, would increase the death rate out of the pool and result in a smaller pool size. The finite size of the risk pool creates the possibility of a negative association with pollution at some lags. This rebound (*ie*, drop in the number of deaths, after an initial increase) presupposes that air pollution does not affect recruitment into the pool. Yet numerous epidemiologic studies have shown particulate air pollution to be associated with exacerbation of illness, including increased hospitalizations, ^{14} decreased heart rate variability, ^{15} etc, thus suggesting that increased recruitment is possible. Recently, Zelikoff et al ^{16} have shown that particle exposure exacerbates pneumonia in animals. Hence, air pollution may intensify some illnesses, increasing the size of the risk pool. Further, this may occur with a different lag than that between exposure and death out of the risk pool. Hence, the direction of the effect of an air pollution episode on the size of the risk pool, and the effect of the risk pool on the death rate over time, may be positive or negative.

Recently, three papers have examined this issue indirectly, by estimating the association between air pollution and daily deaths in Philadelphia, ^{17} Boston, ^{18} and Chicago ^{19} after filtering out such rebounds. None of the studies found any evidence that the effect size for air pollution was reduced as a result of the mortality displacement, and indeed all three studies reported that the effect size approximately doubled. Schwartz ^{18} interpreted this as suggesting that, far from depleting the pool of critically ill people, air pollution increased the size of the pool over longer time scales by increasing the intensity of illness in general. None of these studies provided any direct estimate of what the time course of the rise and fall of mortality after exposure might be (*eg*, Figure 1).

One additional analysis has recently been published. ^{20} These authors assumed a model in which air pollution could only deplete the pool of susceptible individuals at high risk of dying and could not increase recruitment into that pool. This is equivalent to assuming that the correlation between air pollution and daily deaths must become negative after a lag of several days. That assumption is a testable hypothesis.

Another recent paper ^{21} applied a different approach that explicitly tests this hypothesis. Zanobetti *et al*^{21} estimated the association of air pollution at multiple lags simultaneously, providing a direct estimate of Figure 1. Because air pollution is generally correlated, putting a large number of lags of a pollutant into a model produces high levels of multicolinearity and unstable results. To counter this problem, these authors used a nonparametric smoothed distributed lag, looking out to 40 days after exposure, to estimate the effect of air pollution on daily deaths in Milan between 1980 and 1989. This constrained the estimated effects of air pollution to vary smoothly with the number of days of lag between exposure and death. This required special software that is not generally available. However, in a sensitivity analysis, they showed that essentially identical results could be obtained using a cubic polynomial distributed lag model, which can be implemented in any Poisson regression package. In both cases, the coefficients of air pollution at each lag are constrained to fit a smooth shape, in which the latter case is a polynomial. If the polynomial is flexible enough to fit the true pattern of the data reasonably well, little bias will be introduced.

We have adopted that approach for a systematic examination of the lag between air pollution and daily deaths in the Air Pollution and Health: A European Approach (APHEA-2) study. ^{22,23} This analysis focuses on particulate air pollution in a multicity hierarchic model.

#### Subjects and Methods

##### Health Data

The APHEA-2 study is a comprehensive, multicenter study that examines the association between air pollution and daily deaths in 30 cities across Europe and associated regions (*eg*, Tel Aviv). Data collection included daily counts of all-cause mortality, excluding deaths from external causes (*International Classification of Diseases*, 9th revision, code >800). The years of study were 1990 through 1997, although mortality data in most cities were available only through 1995 or 1996. In some cases, air pollution data were available only for part of the period.

Because of resource and time constraints, it was decided *a priori* to limit the analysis of mortality displacement to ten cities. To maximize the power of the study, we chose the largest cities in the study, with the stipulation that only one city could be chosen in each country. The ten cities selected were Athens, Budapest, Lodz, London, Madrid, Paris, Prague, Rome, Stockholm, and Tel Aviv. Together, they comprise a population of about 28 million people, which is two-thirds of the population in the full study, and they represent northern Europe, central Europe, and the Mediterranean region. An earlier paper ^{23} examined the association of particulate air pollution in all available cities and addressed the issue of heterogeneity in response. That analysis did not examine the “harvesting” issue addressed in this paper.

Daily measurements of particulate air pollution were provided by each city participating in the APHEA-2 project. Particulate matter was measured as PM_{10} (particulate air matter less than 10 μM in aerodynamic diameter) in four cities, as PM_{13} (particulate air matter with aerodynamic diameter less than 13 μM) in Paris, and PM_{15} (particulate air matter with aerodynamic diameter less than 15 μM) in Rome. The Paris data were assumed to be equivalent to PM_{10} in this study. Rome data were converted to PM_{10} using a site-specific conversion factor based on colocated measurements. ^{24} In Athens, data were routinely collected only on black smoke. Because traffic is the dominant source of particles in Athens, there were some days of colocated PM_{10} and black smoke monitoring that allowed the establishment of a site-specific selective conversion. Also in Lodz only data for black smoke were available, whereas in Budapest the original data were measured as total suspended particulate. In these three cities, data were converted to PM_{10} as a function of both black smoke (total suspended particulate for Budapest) and season, again on the basis of regression modeling with limited PM_{10} data.

We conducted a weighted metaregression with a dummy variable equal to 1 for cities where the other particle measures were converted to PM_{10} on the basis of site-specific calibration. We found a somewhat higher coefficient in the converted cities (1.98% per 10 μg/m^{3} increase in PM_{10} compared with 1.48% in the cities that measured PM_{10}), but the confidence interval for the incremental 0.5% effect was ±1.93%. These results indicate that the coefficients could in fact be 0. Further, three of the five cities where the conversion occurred were in southern Europe, where a previous hierarchic model of all 29 cities in APHEA-2 showed larger coefficients. We conclude that there is little reason to believe the effect estimates differ between the cities where the air pollutant measurement has been converted and the other cities. Hence, results were reported as the effect of PM_{10}. Further details have been previously reported. ^{23}

##### Covariate Control

Generalized additive regression models ^{25} were fitted in each of the ten cities, controlling for seasonal patterns, long-term time trends for weather, influenza epidemics, holidays, and day of the week. The models were built following the APHEA-2 methodology. ^{23} Because of the substantial variability in seasonal patterns and weather between, for example, Stockholm and Tel Aviv, separate models were chosen in each city. All models controlled for temperature and humidity on the same day using nonparametric smooth function. ^{27} In addition, we examined whether nonparametric functions of weather variables on the previous day or up to 3 previous days or the average of a few days improved model fit (defined as lowering the Akaike information criterion ^{28} for the model). We similarly chose the number of degrees of freedom for each weather variable to minimize the Akaike information criterion. This approach has been used and discussed previously. ^{29,30}

Seasonal patterns are controlled because there are unmeasured predictors of death, such as diet, which vary seasonally and have long-term trends over time. Because air pollution also shows seasonal variations and long-term trends, this creates a potential for confounding. Shorter-term fluctuations in diet are unlikely to be correlated with air pollution. Hence, the goal of our smooth function of time is to remove seasonal and long-term fluctuations.

Various smoothing parameters exist for producing residuals with no seasonality. To choose among them, we examined the partial autocorrelation function of the residuals. This is because, although each death is an independent event, seasonal patterns in the mortality data produce correlations between the number of deaths on one day and on the previous day. Eliminating short-term serial correlation is therefore a measure of how successful our seasonal control has been. On the other hand, the use of excessive degrees of freedom for seasonal control induces negative serial correlation in the residuals of the mortality series, ^{31} which can distort the association with air pollution. Therefore, we chose a smoothing parameter for time to reduce the residuals to white noise. Sometimes it was necessary to introduce autoregressive terms to accomplish this. ^{32} This approach has been used in a number of recent studies. ^{6,12,30}

##### Distributed Lag Model

The goal of our analysis was to estimate the dependence of daily deaths (on day *t*) on PM_{10} on that day and up to the previous 40 days. If the pollution-related deaths are only being advanced by a few days to a few weeks, we will see this effect as a negative association between air pollution and deaths several days to several weeks subsequently. The net effect of air pollution, net of any such short-term rebound up to 40 days, is the sum of the effect estimates for all 41 days. In addition, plotting individual effect size estimates *vs* lag number gives us a direct estimate of what Figure 1 really looks like. This is an example of a distributed lag model, which has been described previously. ^{33,34}

For Poisson regression, the unconstrained distributed lag model may be written as:

where *Z*_{t} = pollution variable delayed over time, for *j* = 0 . . . *q* days.

Because this model produces unstable estimates for large *q*, it is common to constrain the coefficients to vary smoothly with lag number. ^{33} A polynomial distributed lag constrains the β_{j} to follow a polynomial pattern in the lag number, that is:

where *j* is the number of lag of delay and *k* is the degree of the polynomial. Further details, including how to estimate the η_{k} in a Poisson model, have been published previously. ^{34} Too much constraint risks bias, producing a distorted shape, whereas too little constraint produces estimates that are too noisy to be informative. Although a cubic polynomial was sufficient to match the results of the smoothed distributed lag in Milan, ^{21} we have chosen a fourth-degree polynomial in this study, to ensure enough degrees of freedom to fit the pattern of response over time. Such a polynomial has enough degrees of freedom to model a curve such as that shown in Figure 1, or any other plausible shape. Therefore, we estimated in each city the five coefficients η_{0} . . . η_{4} for the fourth-degree polynomial that defines the shape of the distributed lag. As a sensitivity analysis, we used a cubic polynomial and an unconstrained distributed lag model. The unconstrained distributed lag model is too noisy to provide any information about the shape of the effect size *vs* lag, but it does give an unbiased estimate of the overall effect. A separate distributed lag model was fit for each of the ten cities.

##### Second-Stage Modeling

The hierarchic model has two stages. In the first stage, the &eegr_{ik} values are estimated in each city *i*, as described in Eqs 1 and 2.

Equation 1 Image Tools |
Equation 2 Image Tools |

In the second stage, we combined the city-specific coefficients η_{ik}, using the multivariate maximum likelihood method. ^{35}

We assume that:

where &eegr_{i} is the vector of η_{k} in city *i*, S_{i} is the estimated variance-covariance matrix in city *i*, and *D* is the random variance-covariance matrix component, reflecting heterogeneity in response among the cities.

After combining the coefficients &_{ik} by city, the combined coefficients by lag (&bgr_{j}) for the distributed lag model were obtained from Eq 2.

To see how the results compare with more traditional models, we fit the same model in each city using as our exposure index the mean PM_{10} concentration on the day of death and the previous day. ^{11,34,36,37} Note that this model is a highly constrained variant of our distributed lag model, with the constraints forcing β_{1} = β_{0}, and β_{2} = β_{3} = . . . = β_{40} = 0. All analyses were done using the S-plus software (Mathsoft Inc, Seattle, WA).

#### Results

Table 1 shows the ten cities, their populations, the study period in each location, and the mean and standard deviation of the number of daily deaths and environmental variables. Further details of the baseline models for each city have been published previously. ^{23}

Table 2 shows, for each city, the estimated regression coefficients of PM_{10} (per 10 μg/m^{3} and its 95% confidence interval) for the traditional model (mean of the current and previous day), and the overall effect from the fourth-degree polynomial, the cubic, and unrestricted distributed lag models. The overall effect is the sum of the β_{j} per 10 μg/m^{3}. It also shows the combined effect estimates across all of the ten cities, based on a random-effect model to combine results across cities.

Apart from Rome, the estimated effect of PM_{10} increased, and in many cities was more than doubled, when the lagged effects were considered, rather than reduced. These results are seen in all of the distributed lag models that we applied, including the unconstrained model.

The reason for this increase is clear from Figure 2, which shows the estimated effect at each lag, and its confidence interval from the fourth-degree polynomial. It shows that the effect of PM_{10} does decrease to close to 0 with a lag of 10 days, but remains positive, and rises again to a second smaller peak, before dying out to 0 by lag 40.

Figure 3 shows the combined effect for the cubic polynomial. The PM_{10} effect decreases with a minimum at 14 days of lag and then rises again. Although they differ in some detail, both figures show the same general pattern. The initial effect declines to 0 with a lag of 1–2 weeks and then shows a second peak.

To test whether the effect at longer lags made an important contribution to the overall effect, we computed the overall effect (and its standard error) for the first 10 days and for days 11–40 before the death. The effect estimate (×1000) was 0.922 ± 0.184 for the first 10 days of exposure, and 0.688 ± 0.261 for the deaths associated with PM_{10} 11–40 days before. Hence, although the exposure in the first week (and indeed the first 2 days) before the event had a stronger impact, the exposure in the preceding month substantially increased the estimate of the overall effect.

#### Discussion

Previous studies have addressed the mortality displacement issue in single-city analysis. Although these studies were both methodologically innovative and produced valuable information on the issue, the heterogeneity of response to air pollution that has been reported in single-city results ^{23} suggests that a multicity approach, in various locations and using a predefined sampling framework, would be quite valuable in furthering discussion of this issue. Such a study would be necessary to obtain reliable estimates of effect size by lag. Our study is the first report to obtain such stable estimates of effect size by lag in multiple locations.

Qualitatively, our study confirms the basic finding of the previous four studies that did not force harvesting to occur: we do not find that most of the effect of air pollution is short-term harvesting. These results have now been shown in five studies using three different methodologies and in 13 of 14 cities, suggesting that the finding is robust. These findings are also consistent with the results of the episode studies. ^{13} Quantitatively, our study also confirms the previous results by showing that the effect size estimate for airborne particles more than doubles when longer-term effects are taken into consideration.

Our study adds several things to the previous literature. One is the weight of ten cities, which were not selected haphazardly or according to having positive results. This gives considerable assurance that the results are not due to a chance selection of the study locations or selection bias. Second, our study provides insight into the shape of the longer-term response to particulate air pollution. In particular, it suggests that the adverse response to pollution persists up to a month or longer. Moreover, the smoothed distributed lag model of Zanobetti *et al*^{21} produced a very similar curve of effect over time in Milan. There was a prolonged response out to a month in that study as well, with the same dip after 1–2 weeks.

The curves shown in Figures 2 and 3 reflect two processes. One is the pattern of risk over time that occurs in an individual after exposure. This is presumably positive definite, as pollution cannot be expected to improve health. The second is the effect of pollution on the sensitive pool, which can be to expand or shrink that pool. One possible explanation for the observed results is that the effects of air pollution persist for over a month (*ie*, longer-term average exposures have cumulative effects), but that this is partially countered by a drop in the size of the frail pool in the week or two after exposure. A second possibility is that the direct effects of air pollution trail off by a week or so, but that enhanced recruitment into the frail pool results in a long tail of excess deaths triggered by other factors. This is an important issue that remains to be investigated. If there is a prolonged increase in individual risks, it should be possible to identify intermediary biomarkers that remain elevated for some time.

The two-fold increase in risk associated with longer time scales is consistent with the report of higher risk estimates in cohort studies ^{38,39} than in previous time-series studies, given that the cohort studies incorporate effects of longer-term exposure. Together with those studies, it suggests that risk assessment based on the short-term associations likely underestimate the number of early deaths that are advanced by a significant amount, and that estimates based on the cohort studies, or studies such as this one, would more accurately assess the public health impact. Nevertheless, it is important to note that the exposure on the day of death and the immediately preceding day have the greatest impact. This finding suggests that there are important short-term influences at work, which is consistent with recent reports of changes in electrocardiogram patterns within hours of exposure to airborne particles. ^{15}

We note that there appears to be heterogeneity in the response to particles evident in Table 2. This heterogeneity in response has been noted in several studies recently. ^{11,37} Exploration of the cause of such heterogeneity is now a major priority. Demographic factors do not appear to be major predictors. ^{11,37} Chronic obstructive pulmonary disease has been noted as an effect modifier in one study. ^{40} The factors responsible for this heterogeneity in the APHEA-2 cities was the focus of an earlier paper ^{23} (which did not address harvesting), and the mean concentration of NO_{2} and the mean temperature appeared to explain most of the variability. Because this analysis is more limited, we have not attempted to repeat those analyses.

#### Acknowledgments

The APHEA-2 collaborative group consists of: K. Katsouyanni, G. Touloumi, E. Samoli, A. Gryparis, Y. Monopolis, E. Aga, and D. Panagiotakos (Greece, coordinating center); C. Spix, A. Zanobetti, and H. E. Wichmann (Germany); H. R. Anderson, R. Atkinson, and J. Ayres (U.K.); S. Medina, A. Le Tertre, P. Quenel, L. Pascale, and A. Boumghar (Paris); J. Sunyer, M. Saez, F. Ballester, S. Perez-Hoyos, J. M. Tenias, E. Alonso, K. Kambra, E. Aranguez, A. Gandarillas, I. Galan, J. M. Ordonez (Spain); M. A. Vigotti, G. Rossi, E. Cadum, G. Costa, L. Albano, D. Mirabelli, P. Natale, L. Bisanti, A. Bellini, M. Baccini, A. Biggeri, P. Michelozzi, V. Fano, A. Barca, and F. Forastiere (Italy); D. Zmirou and F. Balducci (Grenoble, France); J. Schouten and J. Vonk (The Netherlands); J. Pekkanen and P. Tittanen (Finland); L. Clancy and P. Goodman (Ireland); A. Goren and R. Braunstein (Israel); C. Schindler (Switzerland); B. Wojtyniak, D. Rabczenko, and K. Szafraniek (Poland); B. Kriz, M. Celko, and J. Danova (Prague); A. Paldy, J. Bobvos, A. Vamos, G. Nador, I. Vincze, P. Rudnai, and A. Pinter (Hungary); E. Niciu, V. Frunza, and V. Bunda, (Romania); M. Macarol-Hitti and P. Otorepec (Slovenia); Z. Dörtbudak and F. Erkan (Turkey); B. Forsberg and B. Segerstedt, (Sweden); F. Kotesovec and J. Skorkovski (Teplice, Czech Republic).