In environmental epidemiology, we often conduct observational studies in which exposures to environmental agents cannot be controlled by the investigator. Inference about the health effects of the exposures is generally drawn from a statistical model that controls for potential confounders by including these factors as covariates. Confounding bias caused by omitting important confounders from the regression model is the most common threat to the validity of the exposure effect estimates.1–7
This paper illustrates an approach to diagnosing confounding bias under a causal model linking an environmental exposure and health outcome, estimated using spatiotemporal data. To test the model, we decompose the association between the exposure and health outcome into orthogonal components, corresponding to distinct scales of spatial and temporal variation. If the model adequately controls for confounding, then the exposure effect estimates should be similar at the different spatial and temporal scales. We show that the overall exposure effect estimate is a weighted average of the scale-specific exposure effect estimates. Differences among the scale-specific estimates indicate confounding by omitted covariates.
We illustrate our approach in a study of the mortality effect of 12-month exposure to fine particulate matter (PM2.5). We develop a log-linear regression model for multisite time-series data to estimate the association between month-to-month variation in mortality rates and month-to-month variation in average PM2.5 over the preceding year in 113 US counties and for the period 2000–2002. We decompose the association between PM2.5 and mortality into 2 components: (1) the association between “national trends” in PM2.5 and mortality; and (2) the association between county-specific deviations from the national trends, that is, between “local trends.” This second component provides evidence as to whether counties having steeper declines in PM2.5 also have steeper declines in mortality with respect to their national trends.
If monthly mortality rates are caused by average PM2.5 concentration in the previous year, the associations between the national and local trends should be the same, absent confounding and measurement error. Our proposed approach allows us to assess the validity of this causal hypothesis.
We hypothesize that the association between the national trends in PM2.5 and mortality is likely to be confounded by slowly time-varying factors, such as changes in industrial activities and the economy, improving health care, and large scale weather events.8–11 Our approach can be used to focus on the component of association that is least likely to be confounded, the association between the local trends.
The statistical framework proposed in this paper draws from both cohort studies of long-term exposure12–15 and multisite time-series studies of short-term exposure.16–23 As in cohort studies, we focus on long-term average exposure (averaged over the previous year). As in time-series studies, we estimate associations between temporal changes in exposure and outcome within counties, to guard against bias due to county-specific characteristics that do not vary with time.
We construct mortality counts (Ytc) and number of people at risk (Ntc) for each county c and month t for 6 strata (2 sexes and 3 age groups: 65–74 years; 75–84 years; and >85 years), using Medicare enrollment files. Our study population includes 8.2 million Medicare enrollees living on average 6 miles from an EPA PM2.5 monitor.
The locations of the 113 US counties included in the study are shown in Figure 1. The counties are categorized into 7 geographic regions. The regions are based on our previous national multisite time-series studies of PM10 and mortality and of PM2.5 and hospital admissions.16,24 These counties have nearly complete PM2.5 data (no gaps larger than 3 weeks) for the period over which exposure was averaged, 1999 to 2002.
Estimating County-Specific Annual Average PM2.5
For each county and each month, we calculate the average level of PM2.5 over the preceding year (denoted by PMtc) as follows. First, we estimate the smooth trend in PM2.5 using a linear regression model with outcome monitor-specific daily PM2.5 level, and as predictor a natural cubic spline of time with 16 degrees of freedom. Second, for each month, we calculate the average PM2.5 over the previous year using the fitted values from the regression model described above. This modeled PM2.5 allows us to impute small gaps of missing data when calculating annual averages. For counties with multiple monitors, we use the one with the most complete data, that is, the one with the smallest maximum and average gap in observations and with the longest observation period. We use data from a single monitor rather than from all the monitors within a county because averaging ambient PM2.5 concentrations across monitors that are online for varying periods of time might induce spurious trends.
Measurement Error in County-Specific Annual Average PM2.5 Trends
To investigate whether the observed variation in PMtc trends across counties represents true between-county variability in long-term exposure, rather than differences between monitors within a county (“measurement error”), we perform the following analysis. First, for each monitor, we linearly regress PMtc on t and estimate the slope. Here, we use all monitors with at least 80% of the data available and with no gaps longer than 1 month. Second, we fit a 1-way random effects model to the monitor-specific estimated slopes and calculate: (1) the variability of the slopes within county (measurement error) (σw2); (2) the variability of the slopes between counties (σb2); and (3) the intraclass correlation coefficient, ρ=σb2/(σb2+σw2).
Analysis of Variance of County-Specific Annual Average PM2.5 Trends
To quantify the variability of PMtc in space and time, we conduct the following analysis of variance. We fit a linear model with PMtc as the dependent variable, and with the following predictors: (1) county-specific intercepts (the spatial dimension); (2) a natural cubic spline of month with 16 degrees of freedom (the time dimension); and (3) an interaction between the county-specific indicators and the smooth function of time (the space-by-time interaction).
Causal Model for Annual Average PM2.5 and Mortality
Within each age–sex stratum, we consider the following causal model for the health effects of air pollution:
The parameters δ0c are county-specific intercepts, which are included in the model to control for unmeasured county-specific characteristics that do not vary with time. The parameter δ1 denotes the association between month-to-month variation in PMtc and month-to-month variation in mortality.
Estimates from model (1) are likely to be confounded by factors that cause trends in PM2.5 and mortality. Examples of such confounders are policy changes affecting the economy, industrial activity, and health care and large scale weather events.8–11 A popular approach to controlling for unmeasured temporal confounding at the national level is to add to the model a smooth function of time:
where s(t; d) is a smooth function of time modeled using a natural cubic spline with d degrees of freedom. We emphasize that this model controls for temporal trends at the national level, since s(t; d) is common to all counties. The parameters β0c are county-specific intercepts. This model is equivalent to the following:
The term SYMBOL denotes the national trend in annual average PM2.5, calculated as the fitted values of a linear regression model having PMtc as dependent variable (for all counties) and a natural cubic spline of time with d degrees of freedom (s(t; d)) as predictor. The term s*(t; d − 1) is a smooth function of time modeled using a natural cubic with d − 1 degrees of freedom, orthogonal to SYMBOL and PMtc.
Models (2) and (3) yield the same predicted values. The only difference between the 2 models is in parametrization: model (3) takes the smooth function s(t; d) in model (2), which is represented by a set of d basis functions, and breaks it into: (1) SYMBOL, which is a linear combination of the d basis functions; and (2) the remaining smooth function, s * (t; d − 1). The parameters η2 in model (3) and β1 in model (2) are exactly the same.
Model (3) allows us to estimate the association between PM2.5 and mortality trends at 2 different scales: national and local. The parameter η1 denotes the association between month-to-month variation in the national trend in SYMBOL and month-to-month variation in the national trend in mortality rates. The parameter η2 denotes the association between month-to-month variation in county-specific deviations in PMtc from the national trend, and month-to-month variation in county-specific deviations in mortality from the national trend. In other words, η2 provides evidence as to whether counties having steeper declines in SYMBOL also have steeper declines in mortality relative to the national trend.
If model (1) describes the causal link between annual average PM2.5 and mortality, then the estimates of η1 and η2 in model (3) should be equal, absent confounding and measurement error. Therefore, a comparison of η̂1 ofη̂2 provides important evidence on the causal hypothesis formulated in model (1).
In model (3), the term SYMBOL controls for the national trend in annual average PM2.5, and s*(t;d − 1) controls for the remaining national trend in mortality. This implies that the effect of SYMBOL (η1), which represents the association between trends in PM2.5 and mortality at the national scale, is potentially confounded by time-varying factors such as changes in the economy and health care. We focus on η2, the association between trends in PM2.5 and mortality at the local scale, because we believe that this exposure effect is less likely to be confounded. To bias the estimation of η2, a confounder must cause county-specific deviations in PMtc and mortality from their national trends. An example of such a factor is “health consciousness,” a characteristic of counties that relates to their aggressiveness in implementing national air pollution regulatory standards and in improving health care.
It can be shown that the PM2.5-mortality association as measured by model (1) is a composite of 2 pieces of information:
where η̂1 and η̂2 are the estimated coefficients of the national and local PM2.5 trends from model (3), w = (1/V1)/(1/V1 + 1/V2), and V1 and V2 are the statistical variances of η̂1 and η̂2. That is, δ̂1 is a weighted average of the association between the national PM2.5 and mortality trends and the association between the local PM2.5 and mortality trends.
We also consider a pooled model that combines information across age–sex strata and allows for stratum- and region-specific smooth functions of time:
where α0cs are county and age–sex stratum-specific intercepts, srs(t;d) is a stratum- and region-specific smooth function of time modeled using a natural cubic spline with d degrees of freedom, and α1 is the PM2.5 effect common to all age–sex strata. When d = 0, model (4) is an age–sex stratum pooled version of model (1), and α1 is the association between PM2.5 and mortality without control for trends. When d > 0, model (4) is a pooled version of model (2), or equivalently of model (3). The parameter α1 is the association between month-to-month deviations in PM2.5 and mortality from their respective stratum- and region-specific trends, ie, the association between local trends.
In all log-linear models, we use a negative binomial variance model,25
We fit the models by iterating between fitting the log-linear model for fixed φ, and estimating φ using a method of moments estimator.26
We report results for all models when d = 16 degrees of freedom are used to model the national trends over 3 years.
We assess the sensitivity of the results to different choices of d, from d = 0 to d = 32. We vary d on the log2 scale so as to maintain the same knots as d increases. We also calculate robust standard errors,27 which account for residual autocorrelation in monthly mortality rates. Robust and model-based standard errors are similar, and hence we report only the results using model-based standard errors. We also explore the sensitivity of our results to the time period over which PM2.5 is averaged. We fit the same models, using average PM2.5 over the previous 2 years as exposure.
In the measurement error analysis of PM2.5 trends, we find that 80% of the total variability in monitor-specific trends is attributable to variability among counties.
Table 1 summarizes the results of the analysis of variance of PMtc. We find that 91% of the total variance in PMtc can be attributed to the space component, and 5% to the space-by-time component. Note that the space-by-time variance of PMtc, which provides the main source of information for estimating η2 in model (3), is larger than the variance due to the time component, and accounts for 57% of the temporal variance.
Figure 2A displays regional and national linear trends in annual average PM2.5 concentrations. We estimate these trends by linearly regressing PMtc on t. Figure 2B shows regional and national trends in log mortality rates. These trends are estimated by log-linearly regressing Ytc on t with offset Ntc. The log-linear models are fit separately for each age–sex stratum, and the fitted values are averaged across strata. Annual average PM2.5 concentrations are decreasing over time in all regions except in the Northeast and Central regions. Mortality rates are decreasing in all regions. This information is used to estimate the association between the national trends in PM2.5 and mortality in model (3).
Figure 3A shows how county-specific linear trends in PMtc deviate from the national linear trend. County-specific PMtc trends are calculated by linearly regressing PMtc on t. The deviations are the differences between these county-specific trends and the national trend. The deviations are centered at zero to draw attention to the trends, rather than to the levels. Figure 3B shows how county-specific linear log mortality rate trends deviate from the national linear trend. For each county and age–sex stratum, we calculate the trend in the log mortality rate by log-linearly regressing Ytc on t with offset Ntc. Deviations are the differences between the county- and stratum-specific trends and the national stratum-specific trend. The deviations are centered at zero and averaged across age–sex strata. Three counties with very different trends—Los Angeles county (CA), Peoria county (IL), and De Kalb county (GA)—are identified. This plot examines whether counties in which PM2.5 is decreasing faster than the national trend also have mortality rates decreasing faster than the national trend. In LA county, for example, PM2.5 is increasing relative to the national trend, but mortality is decreasing relative to the national trend. Observe the substantial variability in the county-specific deviations from the national trend. This information is used to estimate the association between local trends in PM2.5 and mortality in model (3).
Figure 4 shows a scatterplot of the slopes estimated by linearly regressing PMtc on t versus the slopes estimated by log-linearly regressing Ytc on t with offset Ntc. The mortality rate slopes are averaged across age–sex strata. Los Angeles county (CA), Peoria county (IL), and De Kalb county (GA) are again highlighted. The median PM2.5 slope is −0.048 (interquartile range [IQR] = 0.056), which corresponds to an average decrease of 0.58 μg/m3 PM2.5 concentration per year (12 × 0.048 = 0.58). The median log mortality rate slope is −2.112 × 10−3 (IQR = 1.904 × 10−3), which corresponds to a 2.50% decrease in the mortality rate each year on average (e12 × −2.112 × 10−3 = 0.9750). We evaluate the association between the PM2.5 slopes and the mortality slopes using a weighted linear regression model, where the weights are the inverse variances of the mortality slope estimates. The regression line is superimposed. There is no evidence of a positive association between the rates of change in PM2.5 and log mortality rates (slope estimate = −0.001; 95% CI = −0.006 to 0.003).
Table 2 displays the results of models (1) and (3), separately for each age–sex stratum. We report results for model (3) when d = 16, but note that any d ≥ 8 provides qualitatively similar results. The first column contains estimates of δ1 from model (1), and the second and third columns show estimates of η1 and η2 from model (3). As expected from Figure 2, we find a strong evidence of an association between national trends in PM2.5 and mortality (second column). However, there is no evidence of an association between local trends in any of the strata (third column). This is consistent with the data displayed in Figure 3 and the exploratory analysis shown in Figure 4.
The first column of Table 2 contains results from model (1). These estimates quantify the association between annual average PM2.5 and mortality without control for temporal confounding. In each age–sex stratum, δ̂1 lies between η̂1 (second column) and η̂2 (third column). This follows from the weighted average result, equation (4). Observe that the positive association between PM2.5 and mortality estimated based on model (1) (δ1) is a combination of a very strong positive association between national trends (η1) and a null association between local trends (η2). The large difference between these 2 effects (η1 and η2) suggests that they should not be combined in a weighted average. In the fourth column of Table 2, we show the weight that is given to the national trend component [(1/V1)/(1/V1+1/V2)]. We find that the national trend component accounts for about 40% of the information contained in δ1.
Figure 5 shows estimates of the association between annual average PM2.5 and mortality based on the pooled model (5), as a function of the degrees of freedom allowed in each stratum- and region-specific trend term per year. When d = 0, we estimate the association without control for temporal confounding. We estimate that a 1 μg/m3 increase in PM2.5 is associated with an 0.86% increase in mortality (95% CI = 0.64% to 1.09%). This corresponds to an 8.96% increase in mortality for each 10 μg/m3 increase in PM2.5, which is remarkably similar to the PM2.5 effect estimated in previous cohort studies.12–14 However, as d > 0 (that is, as we start to control for smooth trends in PM2.5 and mortality), the evidence changes. For d ≥ 8, we find no evidence of an association between local trends in PM2.5 and mortality.
Figure 5 also displays the results of model (4) separately for each year. Again if there is a causal association between exposure and outcome, the estimated association should be similar in different subsets of the data. When d = 0, the 3 year-specific PM2.5 effects are very different, but all statistically significant. The change in mortality associated with a 1 μg/m3 increase in PM2.5 ranges from a 4.02% decrease in 2001 (95% CI = 3.25% to 4.79%) to a 5.30% increase in 2002 (95% CI = 4.41% to 6.19%). As d increases, the 3 year-specific estimates become more similar, and settle around a null effect.
We explore the sensitivity of our results to the time period over which PM2.5 concentrations are averaged, by using PM2.5 averaged over the previous 2 years as exposure (and using mortality data for 2001 and 2002). The results of the age–sex stratum-specific models are shown in Table 3. For model (3), using now just 2 years (24 months) of mortality data, we report results when d = 8 degrees of freedom are used to model the national trends. Results are qualitatively similar for all d ≥ 4. The results shown in Table 3 are qualitatively similar to those in Table 2. We find an association between national trends for most strata, but no association between local trends.
This paper illustrates an approach to the assessment of confounding bias in observational studies where environmental exposures and health outcomes vary in time and space. We introduce a causal model for the association between monthly variations in annual average PM2.5 and mortality rates. We show how this association can be decomposed into 2 components: the association between national trends in PM2.5 and mortality, and the association between local trends in PM2.5 and mortality. We find a very large association at the national scale, and no evidence of association at the local scale. We believe that the national trend component is more likely to be confounded than the local trend component. If we set aside the association between national trends, we are left with no evidence of an effect of PM2.5 on mortality.
Chay et al9 estimated the association between trends in air pollution and adult mortality in the United States using an instrumental-variables approach. Following the Clean Air Act of 1970, counties were designated as “attainment” or “nonattainment” according to their levels of total suspended particulates. These authors compared changes in total suspended particulates levels and mortality rates across attainment and nonattainment counties. They found that, while nonattainment status was associated with large reductions in total suspended particulates in the years 1971–1972, nonattainment status was not significantly associated with reductions in adult or elderly mortality.
In another recent paper, Laden and colleagues28 used extended follow-up data from the Harvard Six Cities Study14 to examine trends in average PM2.5 and mortality rates in 6 US cities. They partitioned time into 2 periods, 1974–1989 and 1990–1998. Controlling for average PM2.5 in the first time period, they found that a reduction in average PM2.5 in the second period was associated with a reduction in the mortality rate.
In our analysis, we define long-term exposure as average PM2.5 over the preceding year. National monitoring data for PM2.5 started in 1999 and therefore we do not have data to estimate exposures for longer time periods. Our sensitivity analysis suggests that, when a different exposure averaging period is used, results do not change qualitatively. Determining the appropriate long-term PM2.5 exposure measure is an important scientific question that deserves further research.
Our analysis focuses on 113 counties with relatively complete PM2.5 data over the study period, and uses data from the best single monitor for each county. We conducted the same analysis using a larger set of 250 US counties (meeting less strict PM2.5 measurement criteria) and using as exposure the annual average PM2.5 concentration averaged across all monitors in each county. This produced very similar results.
In these data, we estimate that 20% of the total variability in PM2.5 trends is within-county variability (measurement error). Using a regression calibration correction,29 we estimate that our PM2.5 local trends coefficient is attenuated by 20% (1 − 0.80, where 0.80 is the intraclass correlation). In contrast, we assume that the national trend in PM2.5 is estimated without error, since it is based on data from 113 counties. We conclude that the attenuation of the local trends coefficient is not enough to explain the discrepancy between the effects of the local and national PM2.5 trends.
Our study, as with most air pollution studies, is potentially affected by various sources of bias. This bias comes from 3 sources. First, we use county-level exposure to represent individual-level exposure. Previous studies have shown that this tends to bias exposure effects towards the null.30,31 The second source of bias is due to the lack of information on area-level time-varying confounders that affect both PM2.5 and mortality trends. We control for such factors by including a smooth function of time in the regression models. The third source of bias is due to the lack of adjustment for individual-level covariates beyond age and sex. However, previous cohort studies have found the air pollution-mortality association to be robust to the adjustment for both time-varying and time-invariant individual-level confounders.32
Our proposed methods can be used more generally to diagnose unmeasured confounding in observational studies where the exposure and outcome vary in time and space. We decompose the exposure variable into orthogonal components and allow each component to have a unique effect on the outcome. If there is a causal link between exposure and outcome, then the exposure components must affect the outcome equally, assuming there is no confounding or covariate measurement error. Therefore, differences in these scale-specific effects are a useful diagnostic tool for assessing confounding and its magnitude. If the exposure effects differ, we suggest focusing on the exposure effects that are thought least likely to be confounded. A priori knowledge about the potential confounders can guide the partitioning: the least confounded exposure effects are those corresponding to scales of variation at which the confounders are approximately constant.
1. Greenland S, Morgenstern H. Confounding in health research. Annu Rev Public Health
2. Christenfeld NJ, Sloan RP, Carroll D, Greenland S. Risk factors, confounding, and the illusion of statistical control. Psychosom Med
3. Dominici F, McDermott A, Hastie T. Improved semiparametric time series models of air pollution and mortality. J Am Stat Assoc
4. Peng RD, Dominici F, Louis TA. Model choice in time series studies of air pollution and mortality. J Royal Stat Soc, Ser A
5. Burnett R, Ma R, Jerrett M, et al. The spatial association between community air pollution and mortality: a new method analyzing correlated geographic data. Environ Health Perspect
6. Jerrett M, Burnett RT, Willis A, et al. Spatial analysis of the air pollution-mortality relationship in the context of ecologic confounders. J Toxicol Environ Health A
7. Touloumi G, Samoli E, Pipikou M, Le Tertre A, Atkinson R, Katsouyanni K. Seasonal confounding in air pollution and health time-series studies: effect on air pollution effect estimates. Stat Med
8. Curriero FC, Heiner KS, Samet JM, et al. Temperature and mortality in 11 cities of the eastern United States. Am J Epidemiol
9. Chay K, Dobkin C, Greenstone M. The Clean Air Act of 1970 and adult mortality. J Risk Uncertainty
10. Chay K, Greenstone M. Air quality, infant mortality, and the Clean Air Act of 1970. Tech. Rep. Working Paper No. 04–08. MIT Department of Economics; 2003.
11. Greenstone M. The impacts of environmental regulations on industrial activity: evidence from the 1970 and 1977 Clean Air Act Amendments and the Census of Manufacturers. J Political Econ
12. Pope CA, Thun MJ, Namboodiri MM, et al. Particulate air pollution as a predictor of mortality in a prospective study of US adults. Am J Respir Crit Care Med
13. Pope CA, Burnett RT, Thun MJ, et al. Lung cancer, cardiopulmonary mortality, and long-term exposure to fine particulate air pollution. J Am Med Assoc
14. Dockery DW, Pope CA, Xu X, et al. An association between air pollution and mortality in six US cities. New Engl J Med
15. Jerrett M, Burnett RT, Ma R, et al. Spatial analysis of air pollution and mortality in Los Angeles. Epidemiology
16. Dominici F, Peng RD, Bell ML, et al. Fine particulate air pollution and hospital admission for cardiovascular and respiratory diseases. J Am Med Assoc
17. Bell ML, Dominici F, Samet J. Time series studies of particulate matter. Annu Rev Public Health
18. Samet JM, Dominici F, Curriero FC, et al. Fine particulate air pollution and mortality in 20 US cities, 1987–1994. New Engl J Med
19. Samet JM, Zeger SL, Dominici F, et al. The National Morbidity, Mortality, and Air Pollution Study, Part I: Methods and Methodological Issues
. Cambridge, MA: Health Effects Institute; 2000. Available at: http://pubs.healtheffects.org/view.php?id=228
20. Samet JM, Zeger SL, Dominici F, et al. The National Morbidity, Mortality, and Air Pollution Study, Part II: Morbidity and Mortality from Air Pollution in the United States
. Cambridge, MA: Health Effects Institute; 2000. Available at: http://pubs.healtheffects.org/view.php?id=118
21. Samoli E, Schwartz J, Wojtyniak B, et al. Investigating regional differences in short-term effects of air pollution on daily mortality in the APHEA project: a sensitivity analysis for controlling long-term trends and seasonality. Environ Health Perspect
22. Samoli E, Touloumi G, Zanobetti A, et al. Investigating the dose-response relation between air pollution and total mortality in the APHEA-2 multicity project. Occup Environ Med
23. Katsouyanni K, Toulomi G, Samoli E, et al. Confounding and effect modification in the short-term effects of ambient particles on total mortality: results from 29 European cities within the APHEA2 project. Epidemiology
24. Peng RD, Dominici F, Pastor-Barriuso R, et al. Seasonal analyses of air pollution and mortality in 100 US cities. Am J Epidemiol
25. McCullagh P, Nelder JA. Generalized Linear Models
. Boca Raton, FL: Chapman and Hall; 1989.
26. Agresti A. Categorical Data Analysis
. Hoboken, NJ: Wiley; 2002.
27. Liang KY, Zeger SL. Longitudinal data analysis using generalized linear models. Biometrika
28. Laden F, Schwartz J, Speizer FE, et al. Reduction in fine particulate air pollution and mortality: extended follow-up of the Harvard Six Cities Study. Am J Respir Crit Care Med
29. Carroll RJ, Ruppert D, Stefanski LA. Measurement Error in Nonlinear Models
. Boca Raton, FL: Chapman and Hall; 1995.
30. Zeger SL, Thomas D, Dominici F, et al. Measurement error in time-series studies of air pollution: concepts and consequences. Environ Health Perspect
31. Dominici F, Zeger SL, Samet JM. A measurement error model for time-series studies of air pollution and mortality. Biostatistics
32. Krewski D, Burnett RT, Goldberg MS, et al. Reanalysis of the Harvard Six Cities Study and the American Cancer Society Study of Particulate Air Pollution and Mortality, Part II: Sensitivity Analyses
. Cambridge, MA: Health Effects Institute; 2000.