Open Access

Research Article

Data and Methods Top

Data. The data used to illustrate the methods of this article are concurrent daily time series of mortality, weather, and particulate air pollution from Cook County for 1987-1994. These are the same Cook County data used by Samet et al. (2000).

The mortality time series are nonaccidental daily deaths of individuals ≥ 65 years old. The mortality time-series data are available from the National Center for Health Statistics (Hyattsville, MD). The weather time series are the daily 24-hr mean temperature and mean dew point temperature. The temperature and dew point temperature are 24-hr averages of hourly observations taken at O'Hare International Airport (Chicago, IL). The weather time-series data are available from the EarthInfo (2003) database.

The air pollution time series is based on the ambient 24-hr concentration of particulate matter < 10 µm in diameter (PM₁₀). The PM₁₀ data are available from the U.S. Environmental Protection Agency (2004) Aerometric Information Retrieval System (AIRS; now Air Quality System Database). The PM₁₀ time series used here is the average of the PM₁₀ concentration of the current day and PM₁₀ concentration of the previous day. If either the current or previous day's PM₁₀ concentration was missing, the average was set equal to the nonmissing value.

The largest 1% of PM₁₀ concentrations was removed from the study to avoid the potential for these points to have an undue influence on the estimated concentration-response functions. This leaves 2,838 days of data for Cook County, with a maximum PM₁₀ concentration of 89 µg/m³.

Methods. In many community time-series studies on the effect of PM₁₀ on mortality, an additive Poisson log-linear model is fit to the time series of observed mortality. Under this model, the daily mortality counts are modeled as independent Poisson random variables with a time varying mean log(µ_t), where

log(µ_t) = µ + confounders_t + ßPM_t . [1]

Here, the subscript t refers to the day of the study; µ_t is the mean number of deaths on day t; µ is the intercept term; confounders_t represents other time-varying variables related to daily mortality. Typical confounders include temperature, humidity, longer-term mortality trends, and seasonality; PM_t is the time series of PM₁₀ concentrations; and ß is the effect of PM₁₀ on mortality. It gives the increase in log(µ_t) per unit increase in PM₁₀.

I allow the effect of PM₁₀ on the logarithm of daily mortality to be a piecewise linear function with one or two change points. At the location of the change points, the slope of the piecewise linear concentration-response function is allowed to change. The difference between this study and other studies is that the piecewise linear functions are constrained to ensure the resulting concentration-response function is biologically plausible.

Piecewise linear functions are commonly used to model the PM₁₀ mortality concentration-response function. These are useful for investigating the possibility of thresholds and other nonlinearities in the concentration-response function. This could yield important information about the effect of particulate air pollution on mortality and be helpful to regulators. For example, a piecewise linear function with one change point allows for the following forms of the concentration-response function:

• Zero slope before the change point and positive slope after the change point. This would indicate that there is a threshold below which PM₁₀ does not affect mortality.

• Positive slope before the change point and zero slope after the change point. This would indicate that there is a saturation point above which there is no incremental effect of PM₁₀ on mortality.

• Positive slope before the change point and smaller positive slope after the change point. This would indicate that there is a point above which the incremental effect of PM₁₀ on mortality is reduced.

• Positive slope before the change point and larger positive slope after the change point. This would indicate that there is a point above which the incremental effect of PM₁₀ on mortality is increased.

A piecewise linear function with two change points allows even more flexibility. For example, it can approximate the logistic curve that has been suggested as a plausible shape for the particulate air pollution concentration-response function (Schwartz and Zanobetti 2000).

To allow for more general forms of the PM₁₀ concentration-response function, model 1 may be extended:

log(µ_t) = µ + confounders_t + g(PM_t). [2]

Where g(PM_t) is some function of PM₁₀, two forms of g(PM_t) will be used in this article:

• Piecewise linear function with one change point:

ß₁ and ß₂ are constrained to be non-negative. These constraints ensure that the piecewise linear function is nondecreasing. θ is the change point and gives the PM₁₀ concentration at which the slope of the piecewise linear concentration-response function is allowed to change. This form of g(PM_t) is referred to as C1. This model fit without constraints on ß₁ and ß₂ is referred to as UC1.

• Piecewise linear function with two change points:

ß₁, ß₂, and ß₃ are constrained to be non-negative. These constraints ensure that the piecewise linear function is nondecreasing. θ₁ and θ₂ are the change points and give the PM₁₀ concentrations at which the slope of the piecewise linear concentration-response function is allowed to change. This form of g(PM_t) is referred to as C2. This model fit without constraints on ß₁, ß₂, and ß₃ is referred to as UC2.

Other forms of the concentration-response function such as smooth functions represented by cubic or natural cubic splines could also be considered. Additional constraints would be imposed on the parameters to ensure that the estimated concentration-response function is nondecreasing.

Once the form of g(PM_t) is chosen, the model parameters are fit using a constrained maximum likelihood. I obtained the constrained maximum likelihood estimates iteratively in two repeated steps. The first step fixes the parameters corresponding to the PM₁₀ terms at their current values and estimates the parameters corresponding to the confounders. This step can be performed using the GLM function available in S-Plus (Insightful, Seattle, WA, USA). The GLM function allows the parameters corresponding to the confounders to be estimated while keeping the PM₁₀ parameters fixed at their current values. The second step fixes the parameters corresponding to the confounders at their current values and estimates the parameters corresponding to the PM₁₀ terms. This step can be performed using constrained optimization software. These two steps are iterated to obtain the final parameter estimates. GLM estimation was chosen to avoid the recent convergence criteria problems associated with the S-Plus GAM function (Dominici et al. 2002b).

Simulation Study Top

I use simulations to explore the bias and mean squared error properties of the constrained, or biologically plausible, concentration-response functions compared with unconstrained concentration-response functions. Mortality time series are generated using two forms of the PM₁₀ concentration-response function:

• Null concentration-response. For the null concentration-response, PM₁₀ has no effect on mortality. In this situation, the log mean of the simulated Poisson mortality time series is given by log(µ_t) = log(83). The average daily mortality in Cook County was 83.

• Dual change-point concentration-response. This is the concentration-response found from fitting model C2 to the actual Cook County data (see "Application"). In this situation, the log mean of the simulated Poisson mortality time series is given by log(µ_t) = log(83) + f(PM_t), where f(PM_t) is the estimated C2 function for Cook County.

For simplicity, confounders were not included in the simulations. Using the observed Cook County PM₁₀ time series for both the null and dual change-point concentration response, 400 mortality time series were generated. Models C2 and UC2 were estimated from the simulated mortality time series. This was done by fitting the Poisson log linear models log(µ_t) = µ + g(PM_t,C2) and log(µ_t) = µ + g(PM_t,UC2) to each simulated mortality time series, where g(PM_t,C2) and g(PM_t,UC2) correspond to models C2 and UC2, respectively.

Figure 1 gives the results of the simulations where models C2 and UC2 were used to estimate the null concentration response. The curve corresponding to the mean estimate indicates that the bias of model C2 is modest. The curves corresponding to the 97.5 and 2.5 percentiles indicate that model C2 has a lower estimation variance than model UC2. The root mean squared error (RMSE) was used to compare the statistical accuracy of the two methods. Figure 1B gives the RMSE of the 400 concentration-response function estimates at each PM₁₀ concentration. The RMSE for model C2 is always lower than the RMSE for model UC2.

Figure 1.

Results of simulations using the null concentration response to generate mortality. (A) Summary of the 400 concentration-response estimates from model C2 and model UC2 estimation. The solid lines correspond to the mean value and 97.5 and 2.5 percentiles of the 400 concentration-response estimates using model C2; the dotted lines are similarly defined for estimates using model UC2. (B) RMSE at each PM₁₀ concentration.

Figure 2 gives the results of the simulations where models C2 and UC2 were used to estimate the dual change-point concentration response. The curve corresponding to the mean estimate for model C2 showed little bias (the mean estimate for model UC2 is an unbiased estimate of the dual change-point concentration response). The curves corresponding to the 97.5 and 2.5 percentiles indicate that model C2 has a lower estimation variance than model UC2. The RMSE for model C2 was again lower than for model UC2.

Figure 2.

Results of simulations using the dual change-point concentration response to generate mortality. (A) Summary of the 400 concentration-response estimates from the model C2 and model UC2 estimation. The solid lines correspond to the mean value and 97.5 and 2.5 percentiles of the 400 concentration-response estimates using model C2; the dotted lines are similarly defined for estimates using model UC2. (B) RMSE at each PM₁₀ concentration.

The simulations showed that the bias of the constrained concentration-response functions was not large and that the bias is a small part of the RMSE. The bias of the constrained concentration-response functions is more than compensated for by a decrease in estimation variance. The smaller RMSE for the constrained concentration-response functions indicates that they offer an improvement in statistical accuracy over unconstrained concentration-response functions. The constrained concentration-response functions have the added advantage of giving biologically plausible concentration-response function estimates.

Application Top

In this section the Cook County data illustrate the methodology of constrained versus unconstrained concentration-response function estimation. Using data from only one city is adequate for this purpose.

A Poisson log-linear model similar to those used in Dominici et al. (2002b) was fit to the data:

log(µ_t) = µ + St₁(time,8/year)
+ St₂(temp₀,6)
+ St₃(temp_1-3,6)
+ St₄(dew₀,3) 
+ St₅(dew_1-3,3) 
+ γDOW_t + g(PM_t). [3]

Here,

• The subscript t refers to the day of the study

• µ_t is the mean number of deaths on day t

• µ is the intercept term

• St₁(time,8/year) is a smooth function of time with 8 degrees of freedom per year

• St₂(temp₀,6) and St₃(temp_1-3,6) are smooth functions with a total of 6 degrees of freedom. temp₀ is the current day's mean 24-hr temperature, and temp_1-3 is the average of the previous 3 days' 24-hr mean temperatures. St₄(dew₀,3) and St₅(dew_1-3,3) are similar functions for the 24-hr mean dew point temperature

• DOW_t is a set of indicator variables for the day of the week

• γ is a vector of coefficients that contains the mortality adjustments for the day of the week

• g(PM_t) is the PM₁₀ concentration-response function.

The smooth functions, St_i(), were represented using natural cubic splines. The use of GLMs and natural cubic splines was one of the methods used to adjust for confounders in the Revised Analyses of the National Morbidity, Mortality, and Air Pollution Study (NMMAPS), Part II (Health Effects Institute 2003). Another common method to adjust for confounders is the use of GAMs, with smoothing splines (Green and Silverman 1994) used to represent the smooth functions St_i(). In light of the revised NMMAPS study, the GAM models are now fit with stricter convergence criteria. GLM estimation with natural cubic splines is essentially a fully parametric version of GAM estimation with smoothing splines. Currently, there is no firm evidence that favors the use of GLM or GAM estimation (Health Effects Institute 2003).

Single change-point models C1 and UC1 were fit to the Cook County data using a sequence of change-point values (θ) ranging from 5 to 90 µg/m³, using increments of 5 µg/m³. As the maximum PM₁₀ concentration is 89 µg/m³, fitting model C1 and UC1 with a change point of 90 µg/m³ is identical to fitting a constrained linear or unconstrained linear concentration-response function, respectively. For model C1 the optimal change-point value for Cook County occurs at 60 µg/m³. For model UC1 the optimal change-point value for Cook County occurs at 65 µg/m³.

A check of the parameter estimates obtained from the constrained iterative estimation procedure is possible using only the S-Plus GLM function. This can be done by creating two new PM₁₀ time series PML and PMU. PML is equal to PM₁₀ if PM₁₀ is < θ and θ otherwise, and PMU is equal to 0 if PM₁₀ is < θ and (PM₁₀-θ) otherwise. The parameters from the GLM fit corresponding to PML and PMU give the slope of the concentration-response function before and after the change point, respectively. For example, to check the estimate of the constrained iterative estimation procedure for θ = 60, the GLM model can be fit with PML being the only PM₁₀ term included in the model. Not including PMU in the model corresponds to the estimated concentration-response function having zero slope after the change point. If the iterative estimation procedure is giving appropriate estimates, the PM₁₀ parameter estimates from the two methods should be the same. This check was used for each value of the change point θ. In each case the PM₁₀ parameter estimates from the two methods agreed.

Figure 3 contains a plot of the estimated C1 function for Cook County and a 95% pointwise confidence interval. The estimated UC1 function is included for comparison. Model C1 says that the effect of PM₁₀ on mortality increases until the change point of 60 µg/m³, but that there is no incremental effect beyond 60 µg/m³. If an actual concentration-response function were shaped like the estimated UC1 function, the effect of PM₁₀ on mortality would decrease for concentrations above 65 µg/m³. This is not biologically plausible.

Figure 3.

Plot of the estimated C1 function with 95% pointwise confidence interval for Cook County. The estimated UC1 function is included for comparison.

Confidence intervals in Figure 3 were produced via parametric bootstrap. Using the estimated ß₁, ß₂, and θ values from model C1, 100 new mortality time series were generated according to the Poisson model [Equation 3]. Model C1 with θ fixed at 60 µg/m³ was then fit to each of the generated mortality time series. At each PM₁₀ concentration, the 95% confidence interval corresponds to the 97.5 and 2.5 percentiles of the 100 estimated C1 concentration-response functions.

Next, the dual change-point models C2 and UC2 were fit to the Cook County data with θ₁ fixed at 25 µg/m³ and θ₂ fixed at 50 µg/m³. The location of these change points is similar to the knot location used by Daniels et al. (2000). The use of two change points allows extra flexibility in the C2 and UC2 concentration-response functions compared with C1 and UC1. The check described above was performed on the constrained iterative estimation procedure for model C2 estimation. Once again the PM₁₀ parameter estimates from the two methods agreed. Figure 4 contains a plot of the estimated C2 function for Cook County and a 95% pointwise confidence interval. The confidence intervals, as for model C1, were produced via parametric bootstrap. Model C2 says that PM₁₀ has no effect on mortality until the first change point of 25 µg/m³ (an actual threshold), has an increasing effect for concentrations between 25 and 50 µg/m³, and has no incremental effect for concentrations beyond 50 µg/m³. If an actual concentration-response function were shaped like the estimated UC2 function (dotted), it would say that the effect of PM₁₀ on mortality is beneficial at some PM₁₀ levels and that higher concentrations of PM₁₀ can be less harmful than lower concentrations of PM₁₀. This is not biologically plausible.

Figure 4.

Plot of the estimated C2 function with 95% pointwise confidence interval for Cook County. The estimated UC2 function is included for comparison.

The above analyses were repeated using 4 instead of 8 degrees of freedom per year for the smooth function of time. Fewer degrees of freedom for the smooth function of time correspond to less aggressive control of time-related confounders such as seasonality. The results of using 4 degrees of freedom were similar to those obtained using 8 degrees of freedom: The shape of the estimated concentration-response functions were similar.

The shape of the estimated C1 and C2 concentration-response functions could have public health implications. The estimated C1 concentration-response function suggests that a PM₁₀ standard would be inefficient unless it were set below the change point of 60 µg/m³. The reason for this is that according to the estimated C1 concentration-response function, PM₁₀ does not have an incremental effect on mortality beyond 60 µg/m³. Similarly, the estimated C2 concentration-response function suggests that a PM₁₀ standard will reduce PM₁₀-induced mortality only if it reduces PM₁₀ concentrations in the range between 25 and 50 µg/m³.

Table 1.

The estimated UC1 and UC2 concentration-response functions were both biologically implausible. This is not an uncommon occurrence. In many studies that examine the shape of the concentration-response function, the estimated concentration-response functions are decreasing in certain air pollution ranges and/or indicate that air pollution at certain concentrations can be beneficial to health (Daniels et al. 2000; Dominici et al. 2002a; Moolgavkar and Luebeck 1996; Smith et al. 2000b). One might ask whether I am overinterpreting the negative slopes of the unconstrained functions UC1 and UC2: Are the negatives slopes actually statistically different from zero? Although the question of whether the negative slopes are statistically significant is not of particular concern, the reason for using constrained concentration-response functions is to preclude concentration-response functions that are believed implausible, irrespective of statistical significance.

Table 1 contains p-values for change in deviance tests. The tests compare the fitted concentration-response functions with a null model of no PM₁₀ concentration-response function and a linear PM₁₀ concentration-response function. The linear concentration-response function fit to the Cook County data has a positive slope, so the constrained and unconstrained versions of the linear concentration-response function are the same. The "Test" column gives the forms of the concentration-response functions being compared. "Null" indicates a model that does not contain a PM₁₀ concentration-response function. "Linear" indicates a model that contains a linear PM₁₀ concentration-response function. The "p-value constrained" column contains the p-values of the change in deviance test for the constrained forms of the concentration-response functions (C1 and C2). The "p-value unconstrained" column contains the p-values for the unconstrained forms of the concentration-response functions (UC1 and UC2). The p-values would not be small enough to reject the linear concentration-response function at conventional significance levels. Of the two constrained concentration-response functions, model C2 gave the best fit to the data, as measured by the likelihood value. Because models C1 and C2 use the same number of parameters (C1 uses a parameter for estimation of the change-point location), this comparison suggests that model C2 is the preferred constrained concentration-response function.

It is possible that a biologically implausible fitted concentration-response function will show a statistically significant improvement over a linear concentration-response function that is nondecreasing (biologically plausible). For example, if we fit an unconstrained single change-point piecewise linear concentration-response function to Cook County with the change point fixed at 65 µg/m³, this would result in the biologically implausible concentration-response function of Figure 3 (dotted line), which shows mortality decreasing markedly for concentrations > 65 µg/m³. The p-value for the test of whether this concentration-response function is an improvement over a linear concentration-response function is 0.03, suggesting that the biologically implausible concentration-response function is an improvement over a biologically plausible linear concentration-response function. However, rejecting the linear concentration-response function in favor of a biologically implausible nonlinear concentration-response function raises problems of interpretation that could be avoided by the use of biologically plausible nonlinear concentration-response functions. These potential anomalies illustrate the importance of having sufficiently flexible parameterizations of biologically plausible concentration-response functions.

Discussion and Conclusion Top

In this article I introduce a method for fitting PM₁₀ mortality concentration-response functions that constrains the concentration-response functions to be biologically plausible, i.e., a nondecreasing function of the PM₁₀ concentration. Simulations showed an increase in statistical estimation precision that results from constraining the concentration-response function estimates to be biologically plausible. Biases in the constrained concentration-response function are small and are more than compensated for by a decrease in estimation variance. This was the case whether or not PM₁₀ was assumed to have an effect on mortality.

The use of constrained concentration-response functions in contrast to unconstrained concentration-response functions was illustrated using data from Cook County. The estimated unconstrained concentration-response functions were biologically implausible. It is difficult to understand and apply unconstrained concentration-response functions because increases in PM₁₀ could imply decreases in mortality, or PM₁₀ could be beneficial to health. The use of biologically plausible concentration-response functions avoids these problems of interpretation without introducing material bias. However, it is important to keep in mind that a concentration-response function not biologically plausible may indicate a more serious problem, such as inadequate adjustment for covariates. In such instances further investigation into the form of the model being used may be required.

The constrained concentration-response functions fit to the Cook County data did not provide evidence against a linear concentration-response function, and similar results have been reported in previous studies (Dominici et al. 2002a; Schwartz and Zanobetti 2000). Smith et al. (2000a) did not find statistically significant evidence in favor of nonlinear concentration-response functions. However, they question whether testing a nonlinear versus a linear concentration-response function is the appropriate formulation. Schwartz and Zanobetti (2000) also present an argument as to why the concentration-response relationship may look like a logistic curve. Other studies have shown evidence in favor of nonlinear particulate air pollution concentration-response functions (Daniels et al. 2000; Smith et al. 2000b).

The fact that for the Cook County data the constrained concentration-response functions were not statistically different from a linear concentration-response function should not affect the desirability of using constrained concentration-response functions. Constrained concentration-response functions are more accurate than unconstrained concentration-response functions, have low bias, and are biologically plausible.

References Top

1. Daniels MJ, Dominici F, Samet JM, Zeger SL. 2000. Estimating particulate matter-mortality dose-response curves and threshold levels: an analysis of daily time-series for the 20 largest US cities. Am J Epidemiol 152:397–406. Find this article online
2. Dominici F, Daniels M, Zeger SL, Samet JM. 2002. . Air pollution and mortality: estimating regional and national dose-response relationships. J Am Stat Assoc 97:100–111. Find this article online
3. Dominici F, McDermott A, Zeger SL, Samet JM. 2002. . On the use of generalized additive models in time-series studies of air pollution and health. Am J Epidemiol 156:193–203. Find this article online
4. Dominici F, Samet JM, Zeger SL. 2000. Combining evidence on air pollution and daily mortality from the 20 largest US cities: a hierarchical modelling strategy. J R Stat Soc Ser A 163:263–302. Find this article online
5. EarthInfo Environmental Databases. 2003. Available: http://earthinfo.com/databases/databases​.htm [accessed 30 January 2004].
6. Green PJ, Silverman BW. 1994. Nonparametric Regression and Generalized Linear Models: A Roughness Penalty Approach. London:Chapman and Hall.
7. Hastie TJ, Tibshirani RJ. 1990. Generalized Additive Models. London:Chapman and Hall.
8. Health Effects Institute. 2003. Special Report: Revised Analyses of Time-Series Studies of Air Pollution and Health. Cambridge, MA:Health Effects Institute.
9. HEI Perspectives. 2002. Understanding the Health Effects of Components of the Particulate Matter Mix: Progress and Next Steps. Cambridge, MA:Health Effects Institute.
10. Ito K, Kinney PL, Thurston GD. 1995. Variations in PM-10 concentrations within two metropolitan areas and their implications for health effects analyses. Inhal Toxicol 7:735–745. Find this article online
11. McCullagh P, Nelder JA. 1989. Generalized Linear Models. 2nd ed. London:Chapman and Hall.
12. Moolgavkar SH, Luebeck EG. 1996. A critical review of the evidence on particulate air pollution and mortality. Epidemiology 7:420–428. Find this article online
13. Pickle LW, Mungiole M, Jones GK, White A. Atlas of United States Mortality. Hyattsville, MD:National Center for Health Statistics, 1996.
14. Samet JM, Dominici F, Curriero FC, Coursac I, Zeger SL. 2000. Fine particulate air pollution and mortality in 20 U.S. cities, 1987-1994. N Engl J Med 343:1742–1749. Find this article online
15. Schwartz J.. 2000. Assessing confounding, effect modification, and thresholds in the association between ambient particles and daily deaths. Environ Health Perspect 108:563–568. Find this article online
16. Schwartz J, Ballester F, Saez M, Perez-Hoyos S, Bellido J, Cambra K, et al. 2001. The concentration-response relation between air pollution and daily deaths. Environ Health Perspect 109:1001–1006. Find this article online
17. Schwartz J, Zanobetti A.. 2000. Using meta-smoothing to estimate dose-response trends across multiple studies, with application to air pollution and daily death. Epidemiology 11:666–672. Find this article online
18. Smith RL, Davis JM, Sacks J, Speckman P, Styer P. 2000. . Regression models for air pollution and daily mortality: analysis of data from Birmingham, Alabama. Environmetrics 11:719–743. Find this article online
19. Smith RL, Davis JM, Speckman P. 1999. Human health effects of environmental pollution in the atmosphere. In: Statistics for the Environment 4: Statistical Aspects of Health and the Environment (Barnett V, Stein A, Turkman F, eds). Chichester, UK:John Wiley and Sons, 91-116.
20. Smith RL, Spitzner D, Kim Y, Fuentes M. 2000. . Threshold dependence of mortality effects for fine and coarse particles in Phoenix, Arizona. J Air Waste Manag Assoc 50:1367–1379. Find this article online
21. U.S. Environmental Protection Agency. 2004. AirData: Access to Air Pollution Data. Available: http://www.epa.gov/air/data/aqsdb.html [accessed 30 January 2004]. Last Updated: February 12, 2004