A conditional expectation approach for associating ambient air pollutant exposures with health outcomes

doi:10.1002/env.978

Environmetrics. Author manuscript; available in PMC 2010 March 25.

Published in final edited form as:

Environmetrics. 2009 March 25; 20(7): 877–894.

doi: 10.1002/env.978

PMCID: PMC2786090

NIHMSID: NIHMS91921

A conditional expectation approach for associating ambient air pollutant exposures with health outcomes

Kathleen A. Wannemuehler,¹ Robert H. Lyles,² Lance A. Waller,² Robert M. Hoekstra,¹ Mitchel Klein,³ and Paige Tolbert³

¹Division of Foodborne, Bacterial and Mycotic Diseases, National Center for Zoonotic, Vectorborne and Enteric Diseases, Centers for Disease Control and Prevention, The Rollins School of Public Health of Emory University

²Department of Biostatistics, The Rollins School of Public Health of Emory University

³Department Environmental and Occupational Health, The Rollins School of Public Health of Emory University

Abstract

Our research focuses on the association between exposure to an airborne pollutant and counts of emergency department visits attributed to a specific chronic illness. The motivating example for this analysis of measurement error in time series studies of air pollution and acute health outcomes was a study of emergency department visits from a 20-county Atlanta metropolitan statistical area from 1993–1999. The research presented illustrates the impact of using various surrogates for unobserved measurements of ambient concentrations at the zip code level. Simulation results indicate that the impact of measurement error on the association between pollutant exposure and a health outcome can be substantial. The proposed conditional expectation approach provided reliable estimates of the association and exhibited good confidence interval coverage for a variety of magnitudes of association. Use of a single-centrally located monitor, the arithmetic average, the nearest-neighbor monitor, and the inverse-distance weighted average surrogates resulted in biased estimates and poor coverage rates, especially for larger magnitudes of the association. A focus on obtaining reasonable exposure measurements within clearly defined subregions is important when the pollutant exposure of interest exhibits strong spatial variability.

Keywords: Air pollution, Conditional expectation, Environmental epidemiology, Maximum likelihood, Measurement error, Spatial variability

1 Introduction and Background

Epidemiological studies of air pollution and health have a common goal of estimating the effect of exposure to specific airborne pollutants on specific health outcomes. These studies involve inherent measurement error issues, as true exposure, whether personal or ambient, is often not available. In many studies the only exposures available are ambient source concentrations from one or more monitoring stations, which often vary both temporally and spatially. The purpose of this paper is to explore and illustrate the impact of using various surrogates for unknown ambient exposures when estimating the association between local ambient exposure and a health outcome.

Errors-in-variables, or predictor measurement error, can arise when an exposure variable of interest can not be measured directly. Extensive coverage of conceptual and analytic issues related to measurement error can be found in the books by Fuller (1987) and Carroll et al. (1995), as well as in reviews such as those by Carroll (1989), Thomas et al. (1993) and Thurigen et al. (2000).

Dominici et al. (2003) provide details on various epidemiological study designs as well as statistical methods for studying health effects of air pollution. Several papers have explored the impact that exposure measurement error has on the estimation of the association between pollutant exposures and health outcomes (Burnett et al., 1994; Brown et al., 1994; Duddek et al., 1995; Zidek et al., 1998; Dominici et al., 2000; Berhane et al., 2004). Zeger et al. (2000) discuss the difference between the true ambient level and the measured ambient concentration, where the latter is the primary focus of the current study. They suggest that, if the Berkson-type error model holds true when relating ambient to personal exposure, this may have little impact on estimates of the association when the focus is on personal exposure. The Berkson-error model has the property that E[True|Observed] = Observed (Thomas et al., 1993). However, from a regulatory perspective, the error relating measurements of ambient concentration with their true spatially-resolved counterparts, may be of greater interest because it is ambient, rather than personal, pollutant levels that may be regulated.

Sheppard et al. (2005) discuss various study design factors that impact what parameter is being estimated in time-series and panel studies of air pollution effects. Their simulation studies suggest that using a single monitor when there is spatial variation in the concentrations results in small but noticeable attenuation of effect estimates, and that using the average of multiple ambient monitor exposures in time-series studies gives estimates that are less biased. Sheppard (2005) discussed further that when estimation of acute health effects is the goal, the use of ambient concentrations in time-series studies is quite adequate.

Carlin et al. (1998) developed spatio-temporal hierarchical models for the association between ozone and emergency department visits rates in Atlanta. They chose universal kriging methods to obtain exposure estimates and standard errors for each zip code centroid. Using measurements from 10 monitors across the Atlanta metropolitan area, Gelfand et al. (2001) interpolated ozone levels within each zip-code region. They consider a fixed spatial case, and then extend to a spatio-temporal model.

2 Motivating Example

The motivating example for this analysis of measurement error in time series studies of air pollution and acute health outcomes was a study of emergency department (ED) visits from a 20-county Atlanta metropolitan statistical area from 1993–1999. Metzger et al. (2004) and Peel et al. (2005) analyzed the association between ambient concentrations of various pollutants and cardiovascular and respiratory outcomes, respectively. They developed single-pollutant models for particulate matter (PM₁₀), ozone, nitrogen dioxide (NO₂), carbon monoxide (CO), and sulfur dioxide (SO₂). Health outcome data consisted of information on emergency department visits of residents of the Atlanta metropolitan statistical area for 31 area hospitals, including ICD-9 diagnostic codes. Counts of ED visits were aggregated across the region for each day and associated with ambient concentrations of single pollutants measured at a centrally located monitor.

For our illustrative analysis, we focus on the Atlanta area within 30 km of the central downtown monitor operated by the Aerosol Research and Inhalation Epidemiology Study (ARIES). We focus our attention upon the pollutant nitrogen dioxide (NO₂) for the year 1999. NO₂ is a primary pollutant that has been shown to exhibit marked spatial variability (Wade et al., 2006) and to be positively associated with various cardiovascular diseases (Metzger et al., 2004) and respiratory diseases (Peel et al., 2005). NO₂ is measured at 4 locations in the above defined area. Data on particle composition and physical characteristics collected by the ARIES station are available from August 1, 1998 forward. ED visits from 94 zip codes having a centroid within this area were included in the analysis. Visits were determined to be associated with cardiovascular disease using reported ICD-9 codes.

The purpose of this analysis is to illustrate the impact of various surrogates for true ambient concentration on the estimation of the association between ambient exposure and ED visits related to a specific illness or disease of interest. In particular, we propose approaches designed to reduce bias due to spatial variability of the pollutant measurements. In the process, we also explore the feasibility of estimating spatial variability using likelihood-based methods when measurements are available from only a small number of locations.

3 Methods

3.1 Health Outcome Model

A common approach (Metzger et al., 2004; Peel et al., 2005) to modeling daily event counts, such as ED visits, involves fitting a generalized linear model with a log-link and Poisson error (McCullagh and Nelder, 1989) such as the following: R_ht ~ P, where

(1)

Here R_ht and n_ht are the number of visits for a specific illness and total number of patient ED visits from subregion h based on residential zip code on day t, respectively. X_ht is the true, unobserved ambient-level concentration of a specific pollutant at the centroid of subregion h on day t. While we use the notation X_ht for simplicity, often the exposure of interest is the previous day’s concentration or possibly a moving-average of current and prior concentrations. [theta]

(t;Ψ) can include covariates such as day-of-the-week and smooth functions of calendar time and weather to model long-term temporal trends, seasonality, and other relevant events such as influenza epidemics, and to account for any additional temporal correlation in the count time series (Dominici et al., 2000; Metzger et al., 2004; Peel et al., 2005; Peng et al., 2006). In equation 1, β₀ + [theta]

(t;Ψ) represents the log(baseline rate) for the outcome in the absence of exposure. β₁, the parameter of interest, is the log relative risk associated with a unit increase in the daily ambient concentration of the pollutant.

Given the usual assumptions, the above model would be be reasonably straightforward except for the fact that X_ht is unknown for all h and t. Ambient air quality measurements are not available for each subregion, but rather there are a limited number of monitors located throughout the region, each providing a daily measure of ambient concentration (Z_mt) at that specific location. One might think of this as a classical missing data problem (Little and Rubin, 2002), albeit an extreme example, however information provided by the correlation structure of the observed measurements can inform us about the unobserved observations. The problem can also be placed in a measurement error framework (Carroll et al., 1995), by considering the set of daily monitor measurements Z_t as a source of potential surrogates for the unknown exposure measurement X_ht. We use the term surrogate as defined by Carroll et al. (1995). That is, we assume f(R|Z,X,Ψ) = f(R|X,Ψ), i.e., non-differential error whereby Z offers no information about R once X is known.

3.2 Surrogate exposures for X_ht

In this section we investigate a conditional expectation approach designed to adjust for measurement error in ambient exposures due to spatial variability. Assuming known mean and variance, this method is equivalent to the use of the simple kriging predictor. See Schabenberger and Gotway (2005) for details on kriging.

Currently published studies have utilized a variety of surrogate measures for ambient exposure. Metzger et al. (2004) and Peel et al. (2005) fit models using a single, centrally-located monitor. Zeger et al. (2000) and Wade et al. (2006) recommend averaging, possibly using a spatial average, over monitors within a specific region. Burnett et al. (1994) divided a geographical region of interest into subregions and assigned each subregion a single monitor.

Wong et al. (2004) compared four different weighted average interpolation methods: spatial averaging, nearest neighbor, inverse distance weighting, and ordinary kriging, to assess each method’s influence on the estimated exposure measurement. Other commonly used surrogates are to utilize a single, centrally located monitor or the arithmetic average of all or a subset of exposure measurements. We take the comparison of surrogates a step further by assessing the impact of each surrogate measure on the estimation of the association between exposure and a health outcome of interest.

We propose the following method where we model the temporal and spatial variability exhibited by the observed ambient concentrations through the mean and covariance structures, respectively. Utilizing estimated parameters, one can then interpolate concentration levels at each location. This method is equivalent to kriging. Modeling the spatial variability is based on describing a region’s air pollution measurements as a spatial random field in which the spatial dependence of measurements at different locations can be expressed through the variance-covariance matrix.

More specifically, let t = 1,…, T denote the day of observation, m = 1,…,M the monitor, h = 1,…,H the subregion, and L = M +H. We can then express Z, the vector of observed concentrations, as:

and X, the vector of unobserved ambient exposure levels, as:

We then assume the following model:

(2)

where

and e ~

(0,∑).

ζ(t;Θ) is a function that defines the regional average exposure for each day. As in model 1, long-term temporal trends are accounted for via smooth functions of calendar time with pre-specified knots. Variables such as day-of-the-week, precipitation, wind direction and magnitude can be included to further account for weather, traffic patterns and any additional temporal correlation in the time series. Let Θ be the vector of parameters corresponding to the time-specific covariates included in the model.

Conditional on the fixed covariates indexed by Θ, we assume that each day’s set of pollutant measurements is an independent realization of the underlying spatial process of the pollutant measurements across a defined region. We seek to explore the impact of various surrogates for the unknown subregion ambient exposures. For this purpose take advantage of parametric isotropic models (Waller and Gotway, 2004) to estimate the spatial covariance parameters. Conditional on the covariates, we assume the spatial correlation between exposure measurements is the same for all pairs of equally distant locations and is not dependent on direction, i.e., we assume stationarity and isotropy, respectively.

Under the defined structure, ∑ is a block diagonal matrix where there are identical blocks of size L = M + H for each day, t. For example, if we assume a spatial exponential covariance structure with a nugget effect given 2 monitors (M = 2) and a single subregion centroid (H = 1), i.e. L = 3, then for day t,

(3)

This model gives rise to the following variance-covariance structure:

where, under an exponential model, equation M8

and d_ll′ is the distance between location l and l′.

Once a structure is chosen for ∑, maximum- or restricted maximum-likelihood methods are available to estimate the covariance parameters. We take advantage of the MIXED procedure’s (SAS, 2004c) ability to fit models with spatial covariance structures. We emphasize that the estimation of the covariance parameters from a parametric model using likelihood methods with a handful of monitors can be challenging, and untenable without the assumption of independence across days. Nevertheless, the following analysis and simulation study illustrate that estimation of the covariance parameters is still possible given only a handful of monitors.

Under the model framework described above, utilization of multivariate normal ( x2133

) theory enables calculation of the conditional expectation (CE) of the unobserved exposure measurements (X) given the observed exposure measurements (Z):

(4)

where Z_a = (z₁₁,…, z_1T, z₂₁,…, z_2T ,…, z_M1,…, z_MT)′ is a vector of observed concentrations from the M monitors on T days, and X_a = (x₁₁,…, x_1T, x₂₁,…, x_2T ,…, x_H1,…, x_HT)′ is a vector of the unknown concentrations. In practice, unknown parameters are replaced by estimates; for example ζ_{Z_a}(t;Θ) and ζ_{X_a}(t;Θ) are both replaced by [zeta]

_{Z_a}(t;Θ), the vector containing estimates of the average exposure for each day across the region based on the observed exposures, z_a, and day-specific covariates. Liang and Liu (1991) suggest that Whittemore’s approach (Whittemore, 1989), which replaces the true value X_a by E[X_a|Z_a], may offer a valid approximation in a non-linear setting if β₁ is small. They also suggest the method produces consistent estimates of the exposure effect β₁, our primary focus of interest, in models with a log-link. The method does not ensure consistent estimation of the intercept β₀.

To obtain the CE estimates, we first estimate the parameters of the mean structure and covariance structure in 2 stages using the MIXED procedure (SAS, 2004c), utilizing only the observed measurements (Z). The first stage uses restricted maximum likelihood (REML) to estimate the vector of parameters, Θ, defining the mean structure. In this first stage, we assume observations of the daily ambient concentration measurements are independent after removing the temporal trends.

Stage-1 Model:

(5)

where ϵ ~

(0, σ²I). The fitted values, = [zeta]

_z(t;

) and the residuals, r = Z − , are obtained. The default optimization algorithm in MIXED is ridge-stabilized Newton-Raphson. We utilized Fisher’s scoring up to the first 50 iterations. In the second stage, a no intercept model is fit to the residuals to obtain REML estimates of the covariance parameters.

Stage-2 Model:

(6)

where e ~

(0,∑), and ∑ takes a spatial exponential form as described in equation 3. Using the estimates of Θ and ∑, the conditional expectation estimates (Eq. 4) are then calculated. See appendix A for details.

4 Simulation Study

4.1 Simulating ambient concentration

True ambient concentrations were simulated for each day of two calendar years for each of five monitored locations and 94 centroid locations, for a total of 99 locations. The actual locations of the 94 Atlanta area ZIP code centroids were used. The five monitor locations were distributed across the region such that one was at the center and the other four towards the outer perimeter. We used a publicly available Excel spreadsheet (Dutch, 2005) to convert monitor and ZIP code centroid locations from latitude-longitude to northing-easting scale based on the Universal Transverse Mercator (UTM) System. This allowed for a scale in kilometers (km) for graphing purposes as well as the unit of measure for the effective range.

Exposure data were simulated as the sum of three sine curves to loosely mimic the short-and long-term trends of NO₂:

where

for all monitor and centroid locations. We chose this method for defining the smooth trend of average daily exposure for its relative ease, as well as for the realism and assessment of using cubic splines to smooth temporal trends when the underlying functional forms of the temporal trends are unknown. Further day-to-day variability is incorporated via the error term, e. Here, e ~ x2133

(0,∑), where ∑ takes a spatial exponential form with a nugget effect as described in equation 3. Covariance parameter values were chosen loosely based on the 1999 NO₂ data from Atlanta: partial sill equation M14

, effective range (ρ_b) = 50 km, and nugget equation M15

. The simulated pollutant measurements had a mean and standard deviation of approximately 50 and 17 ppb, respectively. Figure 1 shows a time-series plot of the simulated exposure measurements with corresponding cubic smoothing splines fit through the data for two of the monitors.

Figure 1

Simulated Ambient Concentration at 2 Monitors

4.2 Simulating Health Outcomes

Approximately 1,000,000 ED visits were dispersed across days and ZIP codes for a 2-year period. Using the observed daily counts and total annual count of ED visits, a rate for each day of the year was calculated. Using the CALL RANTBL routine (SAS, 2004a) and the daily rates, each of the 500678 observed visits were randomly assigned to one of 365 days. Within each day, each visit was then assigned to one of the 94 zip codes using the same routine. Distribution of the probability mass across the 94 zip codes loosely replicated the observed; probabilities ranged from 0.01% to 2.6%. The 2-stage process was repeated using different seed values to obtain counts for a second year. These daily subregion counts of all ED visits, ranging from 0 to 65 with a median of 12, were fixed for all subsequent simulation runs. Daily counts of an event of interest, i.e. ED visits related to a specific cardiac or respiratory event, were simulated from a binomial distribution for each of the 94 subregions (i.e. ZIP codes) using the same day exposure of the true ambient concentrations. That is, R_ht ~ B(n_ht, λ_ht) where n_ht is the total number of ED visits from ZIP code h on day t. ln(λ_ht), the natural log of the probability of the visit being attributed to the health outcome of interest, was defined by the following:

As some subregions had very small total numbers of ED visits, simulating from a binomial distribution insured that the number of events would always be less than or equal to the total number of ED visits for each day in any particular subregion.

The coefficient of the current-day’s exposure, β₁, was set at: 0.001, 0.005, 0.01 and 0.05. We incorporate smooth functions of calendar time via cubic splines with knots on the 21^st of March, June, September and December. See Green and Silverman (1994, Ch.2) or Seber (1977, Ch.8) for details. Indicator variables for day-of-week (DOW) and federal holiday (HOLIDAY) were included in an attempt to more accurately reflect the true nature of ED visits. All terms on the right side of the above equation, except β₁X_ht, were fixed across all simulations.

4.3 Surrogate Exposures and Health Outcome Models

The monitor located near the center of the region was chosen as the central monitor. Utilizing all five monitors’ observed measurements, predictions using each of the four methods were calculated: daily average, inverse-distance weighted, nearest-neighbor, and the conditional expectation (CE) estimates.

Temporal trends were simulated using a mixture of several sine curves (see section 4.1), however in the stage-1 model we remove the long-term temporal trends by fitting smooth functions of calendar time via cubic splines with knots on the 21^st day of each month. Indicator variables for day-of-the-week (DOW) were also included. In the stage-2 model we assume a spatial exponential structure with a nugget effect to estimate the spatial covariance parameters based on modeling the residuals from the stage-1 model. Using the estimated parameters of Θ and ∑, we then calculate the estimated conditional expectations for each day at each location.

Using each of the proposed surrogates, a time-series model was then fit using the GEN-MOD procedure (SAS, 2004b) under an assumed Poisson distribution with a log link as described in equation 1 with covariates HOLIDAY, DOW, and cubic spline with knots on the 21^st day of each month. Invoking the Poisson approximation to a binomial (Dudewicz and Mishra, 1988), the Poisson model for event counts remains reasonable given that events are rare (i.e., small R_ht). For comparison purposes in the simulation study, we also fit the health outcome model using the true CE, based on the known Θ and ∑, and the true centroid subregion measurements, X.

4.4 Simulation Results

Table 1 gives the average and standard deviation of 1000 simulation runs for the estimated covariance parameters. Note the true parameters chosen for the simulation reflect both the spatial correlation as well as the large amount of unexplained variability, either due to unknown explanatory variables or to small-scale variation (e.g. local sources, instrument error) that exist in the real data. Even with only a few monitors, reliable estimation of these spatial covariance parameters is achieved due to the many independent replicates of the spatial process over time.

Table 1

Simulation: Spatial Covariance Parameters

Table 2 gives the average and standard deviation of 1000 simulation runs for the estimated β₁’s in the four different association scenarios, along with the average of the estimated standard errors. Coverage rates, average relative risk (RR) and average corresponding lower and upper confidence bounds are also given for a 20 ppb increase in exposure. Figure 2 and Figure 3 provide histograms of the 1000 estimates of β₁ for a visual comparison of the methods’ performances for β₁ = 0.001 and 0.05.

Table 2

Simulation: Response- Current-day Exposure Model

Figure 2

Simulation (β₁ = 0.001): Histogram of 1000 estimated β₁’s for each method.

Figure 3

Simulation (β₁ = 0.05): Histogram of 1000 estimated β₁’s for each method.

Table 2 indicates that even when β₁ is small, there is considerable bias and sub-optimal coverage of β₁ when using the central monitor and nearest neighbor surrogates. The arithmetic and inverse-distance weighted averages exhibit less bias than both the CM and NN; however, both are clearly more biased than the estimated CE approach. Results also show an increase in the amount of bias and a decrease in the coverage rates for the CM, NN, AVE, and IDW surrogates when the magnitude of the association increases, i.e., as β₁ gets larger. In contrast, the CE approach accurately estimates β₁ for all four magnitudes of β₁, and provides near-optimal coverage for the three smallest magnitudes of association. The bias remains small for 20β₁ = 1 (i.e., β₁ = 0.05). However, the coverage rate suffers in this case, indicating the possible need for a variance adjustment when effect sizes are larger than those seen in real world ambient pollution and health studies. The performance of the estimated CE method is quite respectable, making it the favored approach.

5 Real Example: Atlanta 1999

To further illustrate the utility of the CE approach, we apply the method to 1-hour daily maximum NO₂ data from the Atlanta metropolitan area from the year 1999.

5.1 Exposure Model: NO₂

Prior to estimating the spatial covariance parameters, we fit the following stage-1 model to ambient NO₂ to remove the large scale trend as well as any temporal correlation:

where w_j(t) = (t − τ_j)³ if t ≥ τ_j and w_j(t) = 0 otherwise. ϵ_mt ~ [mathematical script N]

(0, σ²) and Θ = (α₀,…, α₂₃,γ_e,1,…,γ_e,15).

The exposure model includes indicator variables for day-of-the-week (DOW) as well as same-day (RAIN₀) and previous day (RAIN₋₁) precipitation. Also included are average wind direction (WINDDIR), average wind vector magnitude (WINDV EC) and the corresponding interaction between the two. Average wind direction was categorized into four 90° quadrants, measured from 0°N, and wind magnitude was categorized using the quartiles as cut-points. Weather variables (precipitation, wind vector and wind magnitude), measured at Hartsfield-Atlanta International Airport, were obtained from the National Climatic Data Center network.

Cubic splines were included to model the temporal trends with knots on the 21^st of each month. We also included the previous day’s average exposure (_t−1) to remove any residual autocorrelation not handled by the cubic splines. The residuals obtained from the above model were then modeled in turn, assuming a spatial exponential covariance structure with a nugget effect, to estimate the spatial covariance parameters.

5.2 Health Outcome Model

Using the parameter estimates from the stage-1 and stage-2 exposure models above, the CE estimate was calculated at each centroid for each day. For comparisons, the outcome data was modeled in five ways: the ARIES monitor as the central monitor (CM), the nearest monitor’s exposure for each of the 94 ZIP code centroids, the inverse-distance weighted average, the arithmetic average, and the CE estimate. We assume the following health effects model:

The exposure of interest, equation M19

is the previous day’s ambient concentration. Covariates included are federal holiday and day-of-the-week, as well as smooth functions of calendar time, g₁(time,quarterly), with knots on the 21^st of March, June, September and December and smooth functions for average daily temperature and dew point with knots fixed at the first, second and third quartiles.

5.3 Atlanta Results

Table 3 gives the results from the five models for the model of cardiovascular (CVD) disease visits, which for 1999 had an average daily observed rate of ≈ 0.021. The estimated covariance parameters for the 1999 NO₂ measurements can be found in the footnote of Table 3. As in the simulation study, results show that the point estimates of the association between previous day’s ambient concentration of NO₂ with cardiovascular disease using the conditional expectation surrogate are considerably larger than those using either the central monitor or arithmetic average surrogates, both of which ignore spatial variability. In contrast to the simulation trends, in this particular case the inverse-distance weighted average and the nearest-neighbor surrogate estimates are larger than the CE estimate. We note that in the simulation study where 20β₁ = 0.02, ≈ 30% of the simulations resulted in NN and/or IDW estimates that were greater than the estimates from the CE approach. Because the estimation of β₁ is dependent upon the choice of covariance structure, in practice we recommend a sensitivity analysis to assess the impact of different covariance structures.

Table 3

Atlanta: Health Outcome Model - 1999

6 Discussion

Our results demonstrate that when exposure concentrations exhibit spatial variation across a defined region of interest, using a single centrally-located monitor, the nearest neighbor monitor, or the arithmetic or inverse-distance weighted average of several monitors can lead to substantial bias, generally, underestimation of the effect, even when the true association between exposure and outcome is small and the outcome event is rare. The impact of measurement error on the association between pollutant exposure and a health outcome can be substantial. A focus on obtaining reasonable exposure measurements within clearly defined subregions is important when the pollutant exposure of interest exhibits strong spatial variability. Although there are some limitations and computational pitfalls associated with staying within a likelihood framework, we have shown that having many independent replicates of a defined spatial process over time allows one to estimate spatial covariance parameters reasonably well. Having these estimates then allows one to impute subregion exposures that can be incorporated effectively into the health outcome model. This method provides reliable estimates of the association and exhibits good CI coverage for associations of typical magnitudes. With the real data, we estimated the increase in relative risk for a corresponding 20 ppb increase in exposure to be ≈ 1.06 for cardiovascular disease. Results for similar magnitudes of association from the simulation indicate that we can expect near-nominal CI coverage rates in such a setting.

In fitting the health outcome models discussed here, we assumed that the cubic splines adequately addressed the serial correlation. We refer to Metzger et al. (2004), where results of sensitivity analyses which fit GEE models with a stationary 4-dependent correlation structure indicated minimal serial correlation for various cardiovascular diseases. This assumption may not hold for some health outcomes of interest and thus one would need to explore models that account for the additional correlation, regardless of the surrogate used.

The focus of future work could involve incorporating temporal correlation into the spatial fields in an effort to bridge the likelihood approach presented here and the hierarchical Bayesian spatio-temporal models described in Banerjee et al. (2003).

Acknowledgements

R.H.L. was supported in part by and R01 from the National Institute of Environmental Health Sciences (ES012458), and by an Early Career Award from the American Chemistry Council and the American College of Epidemiology. This work was also supported by grants from: National Institute for Environmental Health Sciences (R01ES11294), U.S. Environmental Protection Agency (R82921301-0), and the Electric Power Research Institute (EP-P27723/C13172). We would like to thank Kristi Metzger, James A. Mulholland, Jennifer Peel, Stefanie E. Sarnat, and Jing Yu for all their helpful input and discussions.

Appendix A

To estimate the conditional expectation we first need to rearrange Z* (eq. 2) such that daily observed measurements are grouped within monitor, followed by the daily unobserved measurements grouped within centroid (see description following Eq. 4). That is,

(7)

ζ(t;Θ) and ∑ are then rearranged in the following way:

The rearrangement of ∑ results in a matrix that is no longer block diagonal, but rather a matrix that consists of L² diagonal matrices of size T × T. The matrices on the diagonal of ∑_a take the form equation M22

I_T, where the diagonal elements of each matrix correspond to equation M23

The off-diagonal matrices of ∑_a take the form equation M24

I_T, where the diagonal elements of each matrix correspond to equation M25

For example, assuming the function f(·) takes a spatial-exponential structure, then for M=2, H=1, and T=2, equation M26

and the components of ∑_a can be expressed as:

If we assumed that ζ_{X_a}(t;Θ) and ∑_a were known, the conditional expectation under a multivariate normal framework, as defined above, would be equivalent to the simple kriging predictor p_sk(Z_a, s₀) (Schabenberger and Gotway, 2005). That is, given the random field Z(s) : s [set membership]

^d, where Z(s) = [Z(s₁),…,Z(s_n)]′ is the vector of observed data at locations s₁,…, s_n, the simple kriging predictor, assuming µ(s) and ∑ are known, is: p_sk(Z; s₀) = E[Z(s₀)|Z(s)] = µ(s₀)+σ′∑⁻¹(Z(s)−µ(s)), at location s₀. This optimal linear predictor is the best linear predictor under squared-error loss (Schabenberger and Gotway, 2005).

In what follows, we provide a transformation using Kronecker products that makes calculation of the conditional expectation feasible, even for the purpose of repeated simulation studies. First, recall that equation M28

Let ζ_{X_a}(t;Θ) = [x with macron]

= (

₁,

₂,…,

_H)′ be a vector of length H × T, where each [x with macron]

_h is a vector, length T, of the expected measurements for each day. And let z^d = z_a −ζ_{Z_a}(t;Θ), where z^d is a vector, length M × T, of each day’s difference between the observed measurement, z_mt, on day t at monitor m and the expected measurement for each day, t. Then, equation M29

. Using Kronecker (direct) product notation, Σ_{X_az_a}=S_{X_az_a} [multiply sign in circle]

I_T, where

and Σ_{Z_a} = S_{Z_a} [multiply sign in circle]

I_T, where

and I_T is an identity matrix of size T×T. Using the properties of Kronecker products, it then follows that equation M32

. Now, if we let equation M33

where each S_h is a row vector of length M, then

To reduce the computation time for calculating the conditional expectation for each day at each location, we can iterate through each location one at a time by calculating E[X_a(h)|Z_a] = _h + (S_h [multiply sign in circle]

I_T)z^d. This enables efficient calculation of the CE estimates.

References

Banerjee S, Carlin BP, Gelfand AE. Hierarchical Modeling and Analysis of Spatial Data. Boca Raton: Chapman and Hall; 2003.
Berhane K, Gauderman WJ, Stram DO, Thomas DC. Statistical issues in studies of the long-term effects of air pollution: The southern California children’s health study. Statistical Science. 2004;19:414–449.
Brown PJ, Le ND, Zidek JV. Multivariate spatial interpolation and exposure to air pollutants. The Canadian Journal of Statistics. 1994;22:489–509.
Burnett RT, Dales RE, Raizenne ME, Krewski D, Summers PW, Roberts GR, Raad-Young M, Dann T, Brook J. Effects of low ambient levels of ozone and sulfates on the frequency of respiratory admissions to Ontario hospitals. Environment Research. 1994;65:172–194. [PubMed]
Carlin BP, Xia H, Devine O, Tolbert P, Mulholland J. Spatio-temoporal hierarchical models for analyzing Atlanta pediatric asthma er visit rates. In: Gastonis C, Kass RE, Carriquiry A, Gelman A, Higdon D, Pauler DK, Verdinelli I, editors. Case Studies in Bayesian Statistics. New York: Springer-Verlag; 1998. pp. 303–320.
Carroll RJ. Covariance analysis in generalized linear measurement error models. Statistics in Medicine. 1989;8:1075–1093. [PubMed]
Carroll RJ, Ruppert D, Stefanski LA. Measurement Error in Nonlinear Models. New York: Chapman and Hall; 1995.
Dominici F, Sheppard L, Clyde M. Health effects of air pollution: A statistical review. International Statistical Review. 2003;71:243–276.
Dominici F, Zeger SL, Samet JM. A measurement error model for time-series studies of air pollution and mortality. Biostatistics. 2000;1:157–175. [PubMed]
Duddek C, Le ND, Zidek JV, Burnett RT. Multivariate imputation in cross-sectinal analysis of health effects associated with air pollution. Environmental and ecological statistics. 1995;2:191–212.
Dudewicz EJ, Mishra SN. Modern Mathematical Statistics. New York: John Wiley & Sons; 1988.
Dutch S. 2005. [accessed 03/04/2008]. http://www.uwgb.edu/dutchs/FieldMethods/UTMSystem.htm.
Fuller WA. Measurement Error Models. New York: Wiley; 1987.
Gelfand AE, Zhu L, Carlin BP. On the change of support problem for spatio-temporal data. Biostatistics. 2001;2:31–45. [PubMed]
Green PJ, Silverman BW. Nonparametric Regression and Generalized Linear Models: A roughness penalty approach. London: Chapman and Hall; 1994.
Liang K, Liu X. Estimating equations in generalized linear models with measurement error. In: Godambe V, editor. Estimating Functions. USA: Oxford University Press; 1991. pp. 47–63.
Little RJA, Rubin DB. Statistical Analysis with Missing Data. ”Second” ed. Hoboken, NJ: John Wiley and Sons, Inc; 2002.
McCullagh P, Nelder JA. Generalized Linear Models. Second ed. Boca Raton: Chapman and Hall/CRC; 1989.
Metzger KB, Tolbert PE, Klein M, Peel JL, Flanders WD, Todd K, Mulholland JA, Ryan PB, Frumkin H. Ambient air pollution and cardiovascular emergency department visits. Epidemiology. 2004;15:46–56. [PubMed]
Peel JL, Tolbert PE, Klein M, Metzger KB, Flanders WD, Todd K, Mulholland JA, Ryan PB, Frumkin H. Ambient air pollution and respiratory emergency department visits. Epidemiology. 2005;16:164–174. [PubMed]
Peng RD, Dominici F, Louis TA. Model choice in time series studies of air pollution and mortality. Journal of the Royal Statistical Society A. 2006;169:179–203.
SAS. SAS OnlineDoc 9.1.3. Cary, NC: SAS Institute Inc.; 2004a. SAS/STAT: CALL RANTBL Procedure.
SAS. SAS OnlineDoc 9.1.3. Cary, NC: SAS Institute Inc.; 2004b. SAS/STAT: GENMOD Procedure.
SAS. SAS OnlineDoc 9.1.3. Cary, NC: SAS Institute Inc.; 2004c. SAS/STAT: MIXED Procedure.
Schabenberger O, Gotway CA. Statistical Methods for Spatial Data Analysis. Boca Raton: Chapman and Hall; 2005.
Seber GAF. Linear Regression Analysis. New York, NY: John Wiley & Sons; 1977.
Sheppard L. Acute air pollution effects: consequences of exposure distribution and measurements. Journal of Toxicology and Environmental Health. 2005;68:1127–1135. [PubMed]
Sheppard L, Slaughter JC, Schildcrout J, Liu LJ, Lumley T. Exposure and measurement contributions to estimates of acute air pollution effects. Journal of Exposure Analysis and Enviromental Epidemiology. 2005;15:366–376. [PubMed]
Thomas D, Stram D, Dwyer J. Exposure measurement error: Influence on exposure-disease relationships and methods of correction. Annual Review of Public Health. 1993;14:69–93. [PubMed]
Thurigen D, Spiegelman D, Blettner M, Heuer C, Brenner H. Measurement error correction using validation data: a review of methods and their applicability in case-control studies. Statistical Methods in Medical Research. 2000;9(5):447–474. [PubMed]
Wade KS, Mulholland JA, Marmur A, Russell AG, Hartsell B, Edgerton E, Klein M, Waller L, Peel JL, Tolbert PE. Effects of instrument precision and spatial variability on the assessment of the temporal variation of ambient air pollution in Atlanta, Georgia. Journal of the Air and Waste Management Association. 2006;56:876–888. [PubMed]
Waller LA, Gotway CA. Applied Spatial Statistics for Public Health Data. Hoboken, NJ: John Wiley & Sons; 2004.
Whittemore AS. Errors-in-variables regression using stein estimates. The American Statistician. 1989;43:226–228.
Wong DW, Yuan L, Perlin SA. Comparison of spatial interpolation methods for the estimation of air quality data. Journal of Exposure Analysis and Environmental Epidemiology. 2004;14:404–415. [PubMed]
Zeger SL, Thomas D, Dominici F, Samet JM, Schwartz J, Dockery D, Cohen A. Exposure measurement error in time-series studies of air pollution: Concepts and consequences. Environmental Health Perspectives. 2000;108:419–426. [PMC free article] [PubMed]
Zidek JV, White R, Le ND, Sun W, Burnett RT. Imputing unmeasured explanatory variables in environmental epidemiology with application to health impact analysis of air pollution. Environmental and Ecological Statistics. 1998;5:99–115.