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 (PM_{2.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 PM_{2.5} over the preceding year in 113 US counties and for the period 2000–2002. We decompose the association between PM_{2.5} and mortality into 2 components: (1) the association between “national trends” in PM_{2.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 PM_{2.5} also have steeper declines in mortality with respect to their national trends.

If monthly mortality rates are caused by average PM_{2.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 PM_{2.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 exposure^{12–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.

#### METHODS

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 PM_{2.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 PM_{10} and mortality and of PM_{2.5} and hospital admissions.^{16,24} These counties have nearly complete PM_{2.5} data (no gaps larger than 3 weeks) for the period over which exposure was averaged, 1999 to 2002.

##### Estimating County-Specific Annual Average PM_{2.5}

For each county and each month, we calculate the average level of PM_{2.5} over the preceding year (denoted by *PMtc*) as follows. First, we estimate the smooth trend in PM_{2.5} using a linear regression model with outcome monitor-specific daily PM_{2.5} level, and as predictor a natural cubic spline of time with 16 degrees of freedom. Second, for each month, we calculate the average PM_{2.5} over the previous year using the fitted values from the regression model described above. This modeled PM_{2.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 PM_{2.5} concentrations across monitors that are online for varying periods of time might induce spurious trends.

##### Measurement Error in County-Specific Annual Average PM_{2.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) (*σw*^{2}); (2) the variability of the slopes between counties (*σb*^{2}); and (3) the intraclass correlation coefficient, ρ=σ*b*^{2}/(σ*b*^{2}+σ*w*^{2}).

##### Analysis of Variance of County-Specific Annual Average PM_{2.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 PM_{2.5} and Mortality

Within each age–sex stratum, we consider the following causal model for the health effects of air pollution:

The parameters δ_{0}*c* 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 PM_{2.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 β_{0}*c* are county-specific intercepts. This model is equivalent to the following:

The term SYMBOL denotes the national trend in annual average PM_{2.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 PM_{2.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 PM_{2.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 PM_{2.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 PM_{2.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 PM_{2.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 PM_{2.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 PM_{2.5} trends from model (3), *w* = (1/*V*_{1})/(1/*V*_{1} + 1/*V*_{2}), and *V*_{1} and *V*_{2} are the statistical variances of η̂_{1} and *η̂2*. That is, δ̂_{1} is a weighted average of the association between the national PM_{2.5} and mortality trends and the association between the local PM_{2.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 α_{0}*cs* are county and age–sex stratum-specific intercepts, *s*^{rs}(*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 PM_{2.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 PM_{2.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 PM_{2.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.

##### Sensitivity Analyses

We assess the sensitivity of the results to different choices of *d*, from *d* = 0 to *d* = 32. We vary *d* on the log_{2} 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 PM_{2.5} is averaged. We fit the same models, using average PM_{2.5} over the previous 2 years as exposure.

#### RESULTS

In the measurement error analysis of PM_{2.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 PM_{2.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 PM_{2.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 PM_{2.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 PM_{2.5} is decreasing faster than the national trend also have mortality rates decreasing faster than the national trend. In LA county, for example, PM_{2.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 PM_{2.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 PM_{2.5} slope is −0.048 (interquartile range [IQR] = 0.056), which corresponds to an average decrease of 0.58 μg/m^{3} PM_{2.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 (*e*^{12 × −}^{2.112} ^{× 10−3} = 0.9750). We evaluate the association between the PM_{2.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 PM_{2.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 PM_{2.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 PM_{2.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 PM_{2.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/*V*_{1})/(1/*V*_{1}+1/*V*_{2})]. 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 PM_{2.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/m^{3} increase in PM_{2.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/m^{3} increase in PM_{2.5}, which is remarkably similar to the PM_{2.5} effect estimated in previous cohort studies.^{12–14} However, as *d* > 0 (that is, as we start to control for smooth trends in PM_{2.5} and mortality), the evidence changes. For *d* ≥ 8, we find no evidence of an association between local trends in PM_{2.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 PM_{2.5} effects are very different, but all statistically significant. The change in mortality associated with a 1 μg/m^{3} increase in PM_{2.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 PM_{2.5} concentrations are averaged, by using PM_{2.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.

#### DISCUSSION

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 PM_{2.5} and mortality rates. We show how this association can be decomposed into 2 components: the association between national trends in PM_{2.5} and mortality, and the association between local trends in PM_{2.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 PM_{2.5} on mortality.

Chay et al^{9} 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 colleagues^{28} used extended follow-up data from the Harvard Six Cities Study^{14} to examine trends in average PM_{2.5} and mortality rates in 6 US cities. They partitioned time into 2 periods, 1974–1989 and 1990–1998. Controlling for average PM_{2.5} in the first time period, they found that a reduction in average PM_{2.5} in the second period was associated with a reduction in the mortality rate.

In our analysis, we define long-term exposure as average PM_{2.5} over the preceding year. National monitoring data for PM_{2.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 PM_{2.5} exposure measure is an important scientific question that deserves further research.

Our analysis focuses on 113 counties with relatively complete PM_{2.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 PM_{2.5} measurement criteria) and using as exposure the annual average PM_{2.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 PM_{2.5} trends is within-county variability (measurement error). Using a regression calibration correction,^{29} we estimate that our PM_{2.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 PM_{2.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 PM_{2.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 PM_{2.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.