# Using Meta-Smoothing to Estimate Dose-Response Trends across Multiple Studies, with Application to Air Pollution and Daily Death

Air pollution has been associated with daily mortality in numerous studies over the last decade. Although considerable attention has focused on issues of potential confounding in these associations, little has been done to address the question of what the shape of the dose-response relation looks like. The question of whether a threshold exists for these relations is of particular concern, with regard to both this application and many other epidemiologic questions. Nonparametric smoothing is widely used to control for the potentially nonlinear relations between covariates and daily deaths but has been little used to model the air pollution associations. Because sampling variability, among other factors, can introduce considerable noise into the estimates of linear dose-response curves, quantitative summaries have been widely used to come up with best linear fits. The same ability of meta-analytic techniques to average out noise applies to nonparametric smooth estimates in individual cities. We have developed a method of applying these techniques to combining nonparametric smooths. Using simulation studies, we show that this method can detect threshold and other nonlinear relations in epidemiologic studies, and we then apply it to analyze the association between PM_{10} and daily deaths in ten U.S. cities. We find that the association appears linear down to the lowest levels observed in the study. This method is generally applicable in settings where data from multiple studies can be combined.

From the Environmental Epidemiology Program, Department of Environmental Health, Harvard School of Public Health, Boston, MA.

Supported in part by Environmental Protection Agency Grant R827353 and National Institute of Environmental Health Sciences Grant ES07410.

Submitted November 29, 1999; final version accepted April 14, 2000.

Address correspondence to: Joel Schwartz, Environmental Epidemiology Program, Department of Environmental Health, Harvard School of Public Health, 665 Huntington Avenue, Boston, MA 02115.

Numerous studies have reported associations between particulate air pollution and daily deaths. ^{1–3} The magnitude of the regression coefficients, when combined with the mean exposure levels reported by the U.S. Environmental Protection Agency, suggests that ambient particles are associated with about 60,000 early deaths per year. ^{1} This figure is large relative to the effects of many diseases and suggests that particle control is an important public health problem. The above estimate of effect, however, depends on the assumption of a linear no-threshold relation between daily levels of airborne particles and daily deaths. The public health significance also depends on the assumption that the deaths are being advanced by more than a few days. Two recent papers have addressed the question of the extent of mortality displacement ^{4,5} due to airborne particles and demonstrated that the deaths are displaced by nontrivial amounts. Relatively few papers have addressed the first question, however.

Assessments of the overall literature have generally involved quantitative summaries of the air pollution slopes, and several meta-analyses have been reported. ^{6} Meta-analysis is a quantitative tool to combine slopes or effect size estimates. In either case, what is ultimately being combined is the regression coefficient of a linear predictor, which may be air pollution or a transform of air pollution (for example, its logarithm). The advantage of meta-analysis is that it provides a quantitative summary of the literature. By averaging over many estimates, it cancels out noise and provides a more reliable estimate of the effect size. By combining results from many studies, it is less subject than single studies to accidental confounding, and it allows the use of meta-regression to examine factors that influence the exposure-response slope. This approach assumes that the dose-response relation is linear. If this is not the case, biased estimates can result. Meta-analyses often have to deal with issues relating to comparability of studies, publication bias, etc. If meta-analysis is merely the second stage of a hierarchical model, in which a series of studies is planned to be undertaken using an identical methodology for each study and combined in a second step, these objections vanish. The problem of potential nonlinearity remains, however.

If the dose-response relation is not linear, there may be important public health implications. For example, if there is a threshold for the effects of an air pollutant, measures that reduce exposure that is already below the threshold will have little or no public health benefit. If the threshold level is high relative to typical ambient concentrations, control measures should focus primarily on eliminating high-exposure days. Hence, there is considerable interest in determining the shape of the concentration-response relation.

Nonparametric smoothing represents an attractive approach for modeling the dependence of one variable on another in single studies. Instead of assuming a particular form for the relation, one can let the data determine the shape. This approach is already the standard method for control of weather and seasonality in time-series studies of daily deaths. ^{7} By extending that approach to the pollutant variable as well, we can obtain an estimate of the concentration-response curve in a particular study.

The meta-analytic approach allows us to combine evidence from multiple cities, producing results that are more stable and robust than those from single-city studies and less subject to stochastic variability. The nonparametric regression approach allows us to estimate the shape of the concentration-response relation in a single city. Necessarily, that information will be noisier in any given exposure range than a simple linear slope. We here propose a methodology, which we call meta-smoothing, for combining smooth curves across cities in a hierarchical model to obtain the benefits of both approaches. We demonstrate its applicability with a simulation study, and we then apply it to estimate the concentration-response relation between particulate air pollution and daily deaths in ten U.S. cities.

## Subjects and Methods

### Data

We selected ten U.S. cities with roughly daily monitoring of particulate matter less than 10 μm in aerodynamic diameter (PM_{10}) to provide adequate numbers for the combined analysis. Most other U.S. cities do not have daily PM_{10} monitoring. The selected 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 computed from National Center for Health Statistics mortality tapes for the years 1986 through 1993. All-cause mortality was used for this analysis. Daily weather data were obtained from the nearest airport weather station, and daily concentrations of PM_{10} were obtained from the U.S. Environmental Protection Agency’s Aerometric Information Retrieval System (AIRS) monitoring network. Further details on these cities and data have been published previously. ^{8}

The computation of PM_{10} exposure in each county was complicated by differing monitoring schedules. Some monitors operate on a daily basis, with the others operating every second, third, or sixth day. If the monitors were simply averaged, the daily mean would jump on days when new monitors were included merely because their annual averages differ from those of the monitoring stations that operate on a daily basis.

The variance of PM_{10} measurements also can differ between monitoring locations. Because the monitors contributing to the average change from day to day, this will result in variations 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, we used the following algorithm. The annual mean was computed for each monitor for each year and subtracted from the daily values of that monitor. We 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. We multiplied this by the standard deviation of all of the (centered) monitor readings for the entire year and added back in the average of all of the monitor’s annual averages. This approach has been discussed previously. ^{4}

The mean of the PM_{10} level on the concurrent and previous days was used as the exposure index. This approach has been shown to give results that are similar to those obtained by fitting a full distributed lag model including all lags out to 5 days. ^{4}

### Analytical Methods

### Simulation Study

To demonstrate that our meta-smoothing approach produces unbiased estimates of the dose-response curves, we conducted a simulation study. We used three different true dose-response curves: a linear relation (to assure that the approach recovers the true relation in the basic regression scenario); a piecewise linear dose-response curve, with a slope of 0 for observations below 20 μg/m ^{3} (which represent the hypothesized threshold relation); and a logarithmic curve, with a decreasing slope at higher concentrations (which represent the type of results that have been reported in many individual city studies with high air pollution levels). In addition, we examined the effect of measurement error in exposure on the threshold model.

Each of these three hypothetical curves was used to generate count data as follows. First we generated random PM_{10} exposure data that was log normally distributed with a geometric mean of 27 and a geometric standard deviation of 2. From that, we generated our daily death counts as a random Poisson variable with mean of 22*exp[*f* (*x*)], where *f* (*x*) was one of the three dose-response curves described above. This calculation was done for ten cities to mimic our study design and for 1,096 days (3 years) of data in each city. For each city, we fit a generalized additive model with a smoothed function of exposure, using a span of 0.7. The predicted values (and their pointwise standard errors) were computed for each 2 μg/m ^{3} increment of PM_{10}, and a meta-analysis was done separately for each increment. This estimated the value of *f* (*x*) at each exposure level by the inverse variance weighted average of the ten predicted values.

To test whether the meta-smooth successfully replicated the true dose-response curve, we regressed the predicted values in each interval against the midpoint of the interval to see whether we recovered the original results. For the linear-no-threshold curve, we expected (on average) to recover the true slope; for the piecewise linear curve, the true slopes above and below the knot point; and in the logarithmic curve, the coefficients of the logarithm of exposure. These steps were repeated 500 times and averaged to test that assumption. These results are called the meta-smooth results.

For comparison, we used the same simulated data but fit the true regression scenario in each of the three situations. That is, for the data generated using a linear dose-response relation, we fit a linear term for exposure in the Poisson regression; for the data generated with a threshold relation, we fit a piecewise linear term with a breakpoint of 20 μg/m ^{3}; and for the data generated with a logarithmic dose-response relation, we regressed the daily counts against the logarithm of exposure. These results are referred to as the model-based results.

Finally, for the threshold analysis, after we simulated the count data from the true exposure, we added two types of measurement error to the exposure: additive error with a standard deviation of 2.7 μg/m ^{3}, and multiplicative error with a standard deviation of 20% of exposure. In addition, we simulated the results using spans of 0.5 and 0.7 for the exposure variable, instead of the 0.7 we used in our baseline analysis.

### Ten-City Analysis

We then applied the meta-smoothing analysis to the ten-city mortality data. For each city, a generalized additive Poisson regression was fit, ^{9,10} modeling the logarithm of the expected value of daily deaths as a sum of smooth functions of the predictor variables. This approach allows the regressions to include nonparametric smooth functions that model the potential nonlinear dependence of daily admissions on weather and season. It assumes that:MATH where *Y* is the daily count of deaths, *E* (*Y*) is the expected value of that count, *Xi values are the covariates, and**S*_{i} values are the smooth (that is, continuously differentiable) functions. For the *S*_{i} values we used loess, ^{11} a moving regression smoother. This approach is now standard in air pollution time series. ^{12}

For each covariate, it is necessary to choose a smoothing parameter that determines how smooth the function of that covariate should be. Three classes 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_{10}. 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. We chose a separate smoothing parameter in each city to remove seasonality and to reduce the residuals of the regression to “white noise”^{13} (that is, to remove serial correlation). We used this approach because each death is an independent event, and autocorrelation in residuals indicates there are omitted, time-dependent covariates, the variation of which 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. Because seasonal patterns induce serial correlation, this approach removes seasonal patterns for the data and has been described previously. ^{4,13} Sometimes it was necessary to incorporate autoregressive terms ^{14} to eliminate serial correlation from the residuals.

The other covariates were temperature, dew point temperature, and barometric pressure on the same day; 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 Criterion. ^{15}

To ensure comparability across the ten cities, we used the same span in each city for our smooth function of PM_{10}. We chose a span of 0.7, which corresponds to about 2.5 degrees of freedom, because we believe that should be sufficient to capture the degree of nonlinearity in the PM_{10}-mortality relation. More degrees of freedom would risk introducing unnecessary noise.

As in the simulation study, we computed the predicted values of the log relative risk of daily death in each city for 2-μg/m ^{3} increments between 5.5 and 69.5 μg/m ^{3}, along with their standard errors. The lower increment was chosen so that all of the cities had data within that increment, and we stopped at 69.5 μg/m ^{3} because above that point each increment did not have data from all ten cities. These predicted values at each increment were then combined across cities using inverse variance weighting. This approach is essentially a form of a hierarchical model. ^{16} The explanatory factor for the effect size at each level is the pollution concentration at that level, and the dose-response curve is not constrained to be linear. Rather, in this case, we simply display it, because the key question is what it looks like.

## Results

### Simulation Study

Table 1 shows the true coefficients used in the simulation study and the estimated coefficients obtained by fitting the regressions to the meta-smooth curves. The empirical 95% confidence limits are also shown. There was no evidence that our approach was biased for the linear and threshold models, but there appeared to be a slight downward bias for the logarithmic model with a span of 0.7.

For each scenario (true linear, true threshold, and true logarithmic), we compared the distribution of results from the 500 simulations using the meta-smooth with the distribution of results we obtained fitting the true model for each of the 500 simulations. We defined efficiency as the ratio of the variance of the model-based results to the variance of the meta-smooth. For a true linear model, our meta-smoothing approach was 41% efficient; for the threshold case, it was 58% efficient in estimating the slope above the threshold; and for the logarithmic case, it was 90% efficient. Hence, using locations with large numbers of days of monitoring data available will be important in applying this approach, because the greater power will overcome the loss of efficiency.

To illustrate the ability of this method to detect the existence of a threshold, Figure 1 shows the estimated dose-response relation from the final simulation of the true threshold dose-response relation. It clearly illustrates the ability to detect such a threshold, in a study of ten locations with 3 years of data each.

Table 1 also shows the results of the two simulations with a true threshold relation and either additive or multiplicative measurement error. Importantly, there was no evidence of upward bias in the relation below the threshold. Hence, measurement error will not falsely lead us to find an association at concentrations below a threshold using our meta-smoothing technique. Further, our meta-smooth was not more sensitive to measurement error, in terms of the amount of downward bias in the coefficient above the threshold, than the model-based approach and did not require the knowledge of the true threshold point to be estimated.

Finally, our sensitivity analysis showed that the threshold results were not very sensitive to the choice of span for the exposure variable. A variance-bias tradeoff was seen. If we used a span of 0.90, the efficiency improved to 63%, but the estimated effect size (0.00143) was slightly biased downward. Choosing a span of 0.5 had the opposite effect. The efficiency was reduced to only 32%, but there was no evidence of bias. The logarithmic model showed more sensitivity. With a span of 0.7, there was some downward bias, with an estimate of 0.0367 *vs* a true 0.0400. This bias was reduced (0.379) for a span of 0.50.

Table 2 shows characteristics of our ten study locations, including the mean daily deaths, and the means of the environmental parameters. In total, there were 914,368 deaths in the study locations during the study period.

Figure 2 shows the estimated dose-response curve between PM_{10} and daily deaths in our ten study locations. There is no evidence of a threshold in the relation, down to the lowest levels of airborne particles observed in our study. The curve indicates that an increase from 5 μg/m ^{3} (its baseline) to 30 μg/m ^{3} is associated with about a 2% increase in daily deaths. This is consistent with the results from regressions that have assumed a linear relation.

## Discussion

We have demonstrated that our meta-smoothing approach is an unbiased and moderately efficient method for estimating the dose-response relation between exposure and outcome, when applied to a hierarchical analysis of multiple locations or studies. Understanding the shape of the dose-response curve is important in environmental and occupational health studies, in which exposure can rarely be reduced to zero, and risk assessments are usually performed to estimate the size of the problem. It is of increasing interest in other areas of epidemiology. Some studies have reported that dietary intakes of vitamins or other micronutrients have beneficial effects, and this finding has led to a substantial increase in the number of people who take nutritional supplements, sometimes at high doses. Understanding of the shape of the dose-response curve in these cases will be important as well.

Other methods are available to examine the shape of dose-response curves, notably the use of quintile plots. As the boundaries of the quintiles are unlikely to be identical across studies, however, it does not facilitate the combining information across studies to obtain more stable results.

One clear limitation of meta-smoothing is that this combination across studies can occur only for the range of exposure that is common to the studies. In the example we have given, this overlap included the range of interest, which went from intermediate to low exposures, and allowed us to look for a threshold. It did not allow us to examine the shape of the high-dose part of the dose-response curve. The moderate efficiency is also a limitation of the analysis. The simulation involved ten studies and 1,096 observations each. The 95% confidence intervals in that simulation were reasonably tight. Nevertheless, in situations in which the individual studies have less power, whether because of fewer observations or for other reasons, the confidence bands will increase. In the specific example, there were roughly twice as many days in each study as in the simulation. In that case, ten cities was a reasonable choice for meta-smoothing. For cases with less power in the individual studies, such as those relating environmental tobacco smoke to lung cancer, more studies would probably be needed. In many cases, however (including that one), more studies are available, and it is not impossible to imagine the reanalysis of 20 studies to address this issue.

In recent years, multilocation analyses of the relation between air pollution and daily counts of deaths or hospital admissions have been published ^{3,13} or are underway [APHEA (Air Pollution Health Effects: A European Approach) II ^{17} and NMMAPS (National Mortality and Morbidity Air Pollution Study) ^{18}]. Hence, there is considerable potential for the application of this methodology to air pollution studies.

Our results indicate that there is no threshold for the effect of PM_{10} on daily deaths. This finding is not unexpected. Classical toxicology assumes that all toxicants have thresholds. This presumption has led to a widespread assumption that a threshold also exists for air pollution exposure in the general population. There are, however, considerable differences between general population studies and toxicology studies. Toxicology studies try to measure the results of a given exposure holding every other risk factor as constant as possible. Hence, they are often done on animals that are clones and of the exact same age. Epidemiology studies, in contrast, aim to take advantage of the diversity of the population with respect to other risk factors, to capture the effects of the specified toxicants on the population as a whole, which may include persons with widely varying risks.

Suppose, for example, that each person had a threshold for serious effects from particulate air pollution exposure, but that those thresholds were different. These differences arise from differences in the existence and intensity of current acute illnesses (for example, pneumonia), differences in the existence and intensity of chronic illnesses (for example, dysrhythmias), and differences in genetic factors predisposing to sensitivity. Given a distribution of thresholds for mortality after exposure to particulate air pollution, the number of deaths seen in a population at a given level of exposure will be the sum of all individuals whose thresholds are at or below that exposure. That is, the exposure-mortality dose response curve will be a cumulative probability curve. Because the distribution of thresholds in the general population is the sum of distributions due to multiple genetic factors, multiple degrees of pre-existing disease and multiple degrees of coexposure that cumulative distribution will tend, by the central limit theorem, to approach the cumulative normal distribution. That is, the dose-response relation will tend to look like a logistic curve. We are clearly dealing in the low-probability limit (daily death rates in U.S. cities are typically 2–3 × 10^{−}^{5}) and low-exposure limit (estimated effects of particles are on the order of 1–2% increases in deaths for an interquartile range change in exposure). Hence, if we use a Taylor series expansion of the logistic curve, we are at the low end where the linear term dominates. That is, we expect to find that the low-exposure dose-response relation is linear.

Considerable empirical evidence exists to support this theoretical argument. For example, only a minority of smokers develops bronchitis even after a lifetime of smoking. Some smokers, however, develop bronchitis after only a few years of exposure. Cigarette smoke is rich in particles, which are thought to play a role in the development of cigarette-induced lung disease. It seems plausible to expect a similar heterogeneity in response to ambient particulate air pollution.

Certainly, there is evidence that pre-existing illness modulates the response to airborne particles. In the London smog episode of 1952, the relative increase in deaths from chronic lung disease was much greater than the relative increase in deaths from most other causes. ^{19} Respiratory illness has been shown to predispose to cardiovascular mortality in response to airborne particles. ^{20} Similarly, animal studies have shown that exposure to combustion particles exacerbates *Staphylococcus aureus* pneumonia in rats ^{21} and influenza infections in mice, ^{22} again suggesting effect modification by infection status. Godleski *et al*^{23} have reported that exposure to concentrated ambient particles produced a 37% mortality rate in bronchitic rats but little response in healthy rats. Clearly, there is a wide difference in thresholds depending on the presence of (and presumably the extent of) lung inflammation due to chronic or acute illnesses. Other factors not yet identified may also be predisposing to the response.

There is also support for these results from previous epidemiologic studies. In a study of air pollution and daily deaths in six U.S. cities, the association persisted when restricted to days when concentrations of PM_{2.5} were below 30 μg/m ^{3}, with no diminution (actually, a modest increase) in the estimated effect size. ^{24} Other studies have also occasionally fit models restricting linear dose-response relations to low concentrations. That approach, however, rapidly runs into size limitations. And although it provides evidence of the persistence of an association at low levels, it provides little information about the shape of the dose-response curve.

In conclusion, we have demonstrated an effective methodology for making use of multicity studies to examine the shape of the dose-response relation between air pollution and daily counts of deaths or morbid events, applied it to a study of ten U.S. cities, and found no evidence of a threshold in the association between PM_{10} and daily deaths. Instead, the association appeared quite linear. This latter result is consistent with theoretical observations and with the experimentally observed wide variations in response to particles by level of pre-existing disease.

A copy of Splus Scripts to fit the meta-smooth is available on request from the authors.

## References

_{10}induced effects upon host mortality. Health Effects Institute Annual Meeting, 1997. p 23.

**Keywords:**

air pollution; smoothing; meta-analysis; PM_{10}; particulate matter; models