A simulation study to quantify the impacts of exposure measurement error on air pollution health risk estimates in copollutant time-series models

Background Exposure measurement error in copollutant epidemiologic models has the potential to introduce bias in relative risk (RR) estimates. A simulation study was conducted using empirical data to quantify the impact of correlated measurement errors in time-series analyses of air pollution and health. Methods ZIP-code level estimates of exposure for six pollutants (CO, NOx, EC, PM2.5, SO4, O3) from 1999 to 2002 in the Atlanta metropolitan area were used to calculate spatial, population (i.e. ambient versus personal), and total exposure measurement error. Empirically determined covariance of pollutant concentration pairs and the associated measurement errors were used to simulate true exposure (exposure without error) from observed exposure. Daily emergency department visits for respiratory diseases were simulated using a Poisson time-series model with a main pollutant RR = 1.05 per interquartile range, and a null association for the copollutant (RR = 1). Monte Carlo experiments were used to evaluate the impacts of correlated exposure errors of different copollutant pairs. Results Substantial attenuation of RRs due to exposure error was evident in nearly all copollutant pairs studied, ranging from 10 to 40% attenuation for spatial error, 3–85% for population error, and 31–85% for total error. When CO, NOx or EC is the main pollutant, we demonstrated the possibility of false positives, specifically identifying significant, positive associations for copollutants based on the estimated type I error rate. Conclusions The impact of exposure error must be considered when interpreting results of copollutant epidemiologic models, due to the possibility of attenuation of main pollutant RRs and the increased probability of false positives when measurement error is present. Electronic supplementary material The online version of this article (doi:10.1186/s12940-016-0186-0) contains supplementary material, which is available to authorized users.


Background
Many epidemiologic studies of the health effects associated with ambient air pollution exposure have focused on evaluating the effects associated with a single pollutant. In reality, humans are exposed to a complex mixture of pollutants that can vary spatially and temporally. The challenges of examining the relationship between multipollutant exposures and health effects include high correlations between pollutants preventing the use of standard statistical models, correlation of measurement errors across pollutants, varying degrees of exposure error across pollutants, the potential for interaction between pollutants, and pollutant mixtures varying by location [1][2][3][4][5]. Further, these epidemiologic studies often use data from centrally located measurement sites as the exposure estimate, which may not accurately capture the exposure variability due to the spatial variability of pollutant concentrations. As an example, analysis of pollutants with primarily local sources (CO, NO x , and EC) which exhibit a large degree of spatial heterogeneity [6][7][8] may be more error prone when central-site (CS) monitor measurements are used to estimate exposure, compared to pollutants primarily originating from regional sources (PM 2.5 , SO 4 , O 3 ) which may exhibit spatially homogeneous patterns and thus have a lower degree of measurement error when CS monitor estimates are used. Fixed-location ambient monitors also do not account for human exposure factors such as where people move in time and space (e.g. in-vehicle concentrations impacting exposure), and concentrations in different microenvironments (e.g. infiltration-related factors). There is therefore the potential for exposure measurement error when CS measurements are used as a surrogate for exposure in epidemiologic studies. The presence of exposure measurement error can lead to effect attenuation and reduced statistical power in the resultant health risk estimate [9].
Billionnet et al. summarized the literature on statistical methods to study the effect of multiple pollutants [4], and Oakes et al. published a review of multipollutant exposure metrics [10]. Both identified a comprehensive understanding of exposure measurement error in the context of multipollutant studies as a key issue to address and investigate further. However while these papers detail work on multipollutant exposure metrics and statistical methods for analyzing multipollutant exposures, previous work to quantify the impact of exposure error on health risk estimates has focused on singlepollutant time-series models [6][7][8][11][12][13][14]. While there are inherent difficulties in examining multipollutant exposures, we still do not have a clear understanding of the relationship of exposure measurement error in a more simplistic two pollutant model [15][16][17].
To our knowledge this is the first work which uses empirical pollutant relationships and health data to quantify the impact of correlated exposure measurement error in copollutant time-series models on resultant health risk estimates in a simulation study. We consider Poisson regression and allow for additive, multiplicative, and correlated measurement errors, which are typically not considered in the classical/Berkson error framework [16]. We use previously described, daily exposure metrics ranging from CS measurements to more complex modeling approaches [18] to calculate empirical relationships between pollutants (PM 2.5 , SO 4 , O 3 , CO, NO x , EC). Relationships between exposure metrics were previously used to calculate estimates of exposure measurement error due to spatial variability of pollutant concentrations, and human exposure factors at the ZIP code level [17], but did not include the use of empirical health data nor analysis of an epidemiological model. Our previous work showed the potential for bias in model coefficients for copollutant models, motivating the current work to examine the degree of attenuation of relative risks (RRs) empirically. In the current paper, the previous work is extended in a simulation study by considering a Poisson time-series analysis of emergency department (ED) visits for respiratory disease in the Atlanta, GA metropolitan area to obtain empirical estimates of the attenuation in health risk estimates due to various sources of exposure measurement error.

Methods
Estimates of exposure and exposure measurement error Three estimates of daily exposure to ambient PM 2.5 , EC, SO 4 , CO, NO x , and O 3 were derived for 193 ZIP codes in the 20-county Atlanta metropolitan area. Pollutantspecific measured concentrations, modeled exposure estimates, and summary statistics have been described previously in detail [17,18]. Briefly the three exposure assessment approaches include: (1) CS measurements (CS), (2) ambient air quality (AQ) estimates obtained by combining simulations from an air quality model for regional background and a dispersion model for the local contribution, and (3) estimates from a stochastic population exposure model (PE). For PE, the contribution from indoor sources was not included because of the desire to examine the association between the health outcomes and exposure to ambient pollution, and for comparability with the CS measurements and AQ estimates. All three approaches estimate exposures to ambient pollution at each ZIP code centroid in the study area. Daily estimates (8-h maximum for O 3 , 24-h average for other pollutants) for 1999-2002 were generated for each approach.
For each ZIP code, three types of exposure error (δ spatial , δ population , δ total ) were calculated as the difference between two exposure metrics. Detailed summary statistics including the magnitude and variance of δ, and betweenpollutant correlations of δ are presented in [17]. Briefly, the ZIP code-specific exposure error due to a lack of spatial refinement in the exposure estimate is represented by δ spatial = AQ -CS, assuming that the difference between the more (i.e. AQ) and less (i.e. CS) spatially refined metrics gives an estimate of the amount of spatial measurement error that is present in the CS metric. This assumption is made because our air quality models add spatial variability to the AQ metric compared with CS measurements, which lack spatial variability because the same CS measurement was used to represent exposure in each ZIP code. Exposure error introduced when human exposure factors are not included in an exposure estimate is represented by δ population = PE -AQ. Our PE metric includes variability due to human exposure factors such as time-location-activity patterns of individuals, commuting patterns, and infiltration of ambient pollutants to indoor environments. A third type of exposure error, δ total = PE -CS, represents the combined exposure error when both spatial variability and human exposure factors are not accounted for. δ total does not represent all potential sources of exposure error that may be present in a study; instead it represents the total exposure error that we were able to assess in this analysis.

Epidemiologic model
To conduct the simulation, health data from a previously described epidemiologic study were used [19]. Briefly, individual-level data on ED visits for respiratory outcomes from 41 hospitals in the 20-county Atlanta area were aggregated to daily counts for each ZIP code. The same Poisson time-series model in S. Sarnat et al. [19] was used in this simulation study to examine the association between the above described exposure metrics and daily counts of ED visits. The basic form of the model was: log where Y kt is the count of ED visits for respiratory outcomes in ZIP code k on day t. For each pollutant (pollution 1 and pollution 2 ), daily averages (8-h maximum for O 3 , 24-h average for all other pollutants) of same-day concentrations were used. To control for spatial autocorrelation in the baseline ED visits across the ZIP codes, an indicator variable for the k th ZIP code (ZIP k ) was used to represent the areas from which ED counts were spatially aggregated. Dummy variables for day of week and holidays (DOW, indexed by m) and for hospital (hospital, indexed by n) were used; the latter accounted for the differing durations of time each hospital was included in the study. Long-term trends and seasonality of health outcomes (time) were controlled for with parametric cubic splines with monthly knots (g(γ 1 , …, γ N ; x)). Meteorology was controlled for using maximum temperature (IO temp ) with indicators for each degree Celsius, a cubic function for dew point (dewpt), a cubic function for minimum temperature (temp), and dummy variables for seasons. When simulating the health effect we assume zero lag, but expect results to be generalizable to single day lags.

Simulation of true exposures and health outcomes
True exposure was defined as an exposure estimate without spatial and/or population measurement error.
To simulate the true daily exposure for each pollutant, linear models for each type of daily exposure were fitted separately for each ZIP code and for each pollutant; for example, for spatial measurement error, the model is given by AQ kt = θ k,1 + θ k,2 * CS kt + ε k,t for ZIP code k and day t, such that both additive (θ 1 ) and multiplicative bias (θ 2 ) between the two exposure metrics were accounted for. See Additional file 1: Table S1 for the mean and standard deviation across ZIP codes of estimates of θ k,1 and θ k,2 . For each Monte Carlo iteration, we first calculated a 'true' exposure mean for each day by using the unrefined measurement as the predictor (i.e., CS for spatial or total measurement error, AQ for population measurement error) and a realization of the measurement error coefficients (θ 1 and θ 2 ) drawn from their asymptotic bivariate normal distribution. Across-ZIP code spatial heterogeneity in the measurement error coefficients was allowed, however for each ZIP code we assumed that additive and multiplicative bias remains constant across days. Finally, the residual of the true exposure (ε t ) was subsequently drawn from a bivariate normal distribution with covariance equal to the ZIP code specific empirical measurement error covariance (see Additional file 1: Table  S2 for a summary of the median across ZIP codes of the correlation of measurement error covariance) for the two pollutants of interest, and then combined with the true exposure mean to obtain the full true exposure value (true_exp). We chose to simulate true exposure from the less refined exposure because our measurement error framework assumes that the true exposure is more variable than the error-prone (unrefined) exposure. The resultant measurement error between the true and errorprone exposures, derived here from empirical data, do not strictly follow the standard classical or Berkson measurement error models (Zeger et al. [16]). Using observed health data from the study described in S. Sarnat et al. [19], a design matrix X (without the two pollutants) was constructed for the Poisson model in Equation (1). The Poisson mean μ for each day and each ZIP code was calculated as: where B is a vector of the estimated regression coefficients from Equation (1), subscripts 1 and 2 denote pollutants 1 and 2, RR 1 and RR 2 are the hypothetical RRs per interquartile range (IQR) increase to allow for comparison across pollutants, true_exp 1 and true_exp 2 are the simulated true exposures, and IQR 1 and IQR 2 are the interquartile ranges for the corresponding true exposure metrics. We assumed RR 1 = 1.05 and RR 2 = 1, and in subsequent text refer to pollutant 1 as the 'main pollutant' (i.e., the pollutant with an assumed health effect) and to pollutant 2 as the 'copollutant' (i.e., the pollutant assumed to have no effect). The assumed RR 1 = 1.05 was chosen based on previously established single-pollutant models for the same data set, specifically the O 3 -respiratory analysis [19]. Poisson counts were drawn from a Poisson distribution of the simulated Poisson mean μ. Two versions of the Poisson time-series model (Equation (1)) were then fitted: one using the simulated Poisson counts with the true exposure, and another using the same simulated Poisson counts, together with the unrefined exposure values (standardized by their corresponding IQRs). This resulted in two sets of estimated log RR:β 1; true and β 2; true for the 'true exposure' scenario, andβ 1; noisy and β 2; noisy for the 'noisy exposure' scenario (i.e., using the unrefined exposure metric). The above simulation was run 1000 times for each pollutant pair and measurement error scenario. The mean over the 1000 estimates for each RR was calculated and used as the overall point-estimate associated with each pollutant pair and measurement error scenario. As an example, the RR for the main pollutant, true exposure scenario, was calculated as: where n indexes the simulation iteration and N = 1000. Standard error of the RR was calculated as the standard error of the 1000 βs. A student's two-sided t-test was used to determine significant differences between RR noisy;a and RR noisy;b for two different copollutants a and b, when paired with the same main pollutant, across Monte Carlo simulations.

Statistical analyses
Mean RR estimates from using the true and noisy exposures were compared to the true RR in order to determine if the presence of measurement error induces bias and a loss of power. Various summary statistics are presented, including percent attenuation, root mean-square error (RMSE), and power/type I error.
Percent attenuation of the RR for the main pollutant was calculated to determine the impact of inclusion of measurement error on attenuation of the RR. RMSE was calculated to assess both bias and loss of precision in the estimates. Lastly, statistical power/type I error was calculated to assess false positive associations. As an example, for pollutant 1, noisy scenario, we define: All statistical analyses and model simulations were completed in R, version 3.0.2 (R Foundation for Statistical Computing; http://www.r-project.org/).

Results
Mean RRs for the main pollutant over 1000 runs of the health model for both the noisy and true exposure metrics are shown in Fig. 1. For all types of measurement error, resultant RR of the main pollutants (pollutant 1) when the health model used the true exposure metric (i.e., the simulated exposure metric without measurement error) is close to 1.05. Similarly, the resultant RR of the copollutants (pollutant 2) using the true exposure metric is close to 1, as expected due to the assumed RR 1 = 1.05 and RR 2 = 1. In contrast the presence of measurement error in the noisy exposure metric results in considerable attenuation for all types of measurement error (Fig. 1, Table 1). For spatial measurement error, we see greater attenuation for RR 1;noisy when a local pollutant (CO, NO x , EC) is the main pollutant (29-40%) compared to a regional pollutant (PM 2.5 , SO 4 , O 3 ; 10-15%) (ranges represent the attenuation across all copollutants models for a specific main pollutant). For population measurement error, there is little attenuation when CO is the main pollutant (3-4%), and substantial attenuation (82-85%) of RR 1;noisy when NO x is the main pollutant. Similarly for total measurement error, the highest level of attenuation of RR 1;noisy is seen when NO x is the main pollutant (85%), and the least attenuation for CO (31-32%). For all other pollutants, for both population and total measurement error, attenuation of RR 1;noisy is moderate-high (Table 1). For all pollutants and all types of measurement error, attenuation of RR 1;noisy mirrors the patterns of attenuation  [17] . The scenario with an assumed RR 1 = 1.05 and RR 2 = 1.05 was run as a sensitivity analysis. Resultant RR from the sensitivity analysis are presented in Additional file 1: Figure S1. Two-sided t-tests were used to compare differences in attenuation of the main pollutant RR by copollutant; there were minimal differences seen across copollutants for the case of RR 1 = 1.05 and RR 2 = 1, with more differences seen for the sensitivity analysis (see Additional file 1: Supplemental Text, Table S3, and Table  S4). The coverage of the 95% confidence intervals was calculated for the main pollutant and copollutant for each of the different scenarios analyzed, and with one exception was in the range of 0.93-0.98.
The RMSE ratios comparing the RMSE for the noisy and true model estimates are presented in Table 2, reflecting both bias and loss of precision. The RMSE ratio for the copollutants is typically small (≤2.0), while the RMSE ratio for the main pollutants ranges from 1.3 to 26.1. While the main pollutant RMSE ratio for regional pollutants (PM 2.5 , SO 4 , O 3 ) ranges from 1.8 to 13.7, there is little difference seen across copollutants and measurement error type for the same main pollutant. In contrast, the copollutant can have a substantial impact on the main pollutant RMSE ratio when local pollutants (CO, NO x , EC) are the main pollutant. For example, for δ spatial for CO, NO x , and EC as the main pollutant, the RMSE ratio of the main pollutant ranges from 6.6-11.5, 6.0-10.4, and 6.6 -12.1 respectively, dependent upon the copollutant ( Table 2).
Power estimates for all main pollutants were equal to 1, indicating we will always detect the health effect association for the main pollutant (results not shown) given the simulation setup. Copollutant type I error for the true scenarios (i.e., no measurement error) were approximately 0.05 (results not shown). We present the type I error for the copollutant, noisy scenario (i.e., with measurement error), when the copollutant has no effect (RR 2 = 1) (Fig. 2). Figure 2 shows that for CO, NO x , or EC as the main pollutant there is a substantial increase a b c Fig. 1 Attenuation of RR due to measurement error in a copollutant model (RR 1 = 1.05, RR 2 = 1). For x-axis labels, top row indicates the main pollutant (pollutant 1), bottom row indicates the copollutant (pollutant 2). Overall point estimates shown are the mean over 1000 estimates; error bars indicate the 95 th confidence interval for the 1000 estimates (i.e., the 2.5 th and 97.5 th percentiles of simulated effect estimates; note that extremely narrow confidence intervals result in some non-visible error bars). Note that when copollutant true and noisy model results do not differ substantially, plotted data points overlap. a Spatial measurement error (δ spatial ). b Population measurement error (δ population ). c Total measurement error (δ total )

Discussion
While previous work has examined the impact of exposure measurement error in single pollutant models [13] and studies have proposed statistical methods for handling multipollutant effects [15,[20][21][22], little work has been done to quantify the impact of measurement error on health risk estimates in multipollutant models, including copollutant models. This study builds on previous work examining the potential attenuation of model coefficients in copollutant epidemiologic models based on empirical covariance structures [17]. While the previous work built the foundation for the current analysis In simulations, all main pollutants had RR = 1.05, and all co-pollutants had RR = 1. For each pollutant pair, 1000 simulations were run and results averaged through development of exposure metrics and calculation of empirical relationships between pollutants, the current manuscript extends the work with the addition of empirical health data applied in an epidemiological model with the previously developed exposure metrics. Using the empirical measurement error covariance structures, estimates of additive and multiplicative bias between exposures, and between-pollutant relationships, we extended the previous work by conducting a simulation study to assess the impact of measurement error on health risk estimates obtained from a copollutant timeseries model. We show the substantial attenuation of RRs resulting from the presence of exposure measurement error in exposure estimates for most pollutants and measurement error types examined, an empirical result which agrees with the predictions in previous work [17]. Lastly, we have shown evidence for the detection of false positive associations between the copollutant and the health outcome, particularly when CO, NO x , or EC are the main pollutant. The RMSE ratio comparing the RMSE of the risk estimates from the noisy and true exposures reflects both the bias and standard error of the estimates. When we compare the RMSE ratio for a main pollutant-measurement error type, we see that for regional pollutants as the main pollutant, RMSE ratio changes very little across copollutants. However the RMSE ratio of the local pollutants as main pollutants is impacted substantially depending upon whether regional or local pollutants are the copollutants. For all local pollutant-(as main pollutant) measurement error pairs except for CO-population measurement error, we see a substantial increase in RMSE ratio for the main pollutant when the copollutant is a regional pollutant, compared to when the copollutant is a local pollutant. Thus from examination of the RMSEs (Table 2) we conclude that the inclusion of a regional copollutant in a model with a local main pollutant substantially increases the degree of bias and standard error of the main pollutant estimates, compared to having a local copollutant. However Fig. 1 demonstrates that the copollutant (whether regional or local) has little effect on the degree of bias in the main pollutant's effect estimate. The two results taken together allow us to conclude that the differences in RMSE across copollutant type are due to differences in standard error, rather than bias, with the spatial structure of the copollutant influencing the standard error but not the bias of the main pollutant's effect estimate. We also see differences in the RMSE ratio of local pollutants as main pollutants when the copollutant is varied, even when the copollutants have an assumed null effect (RR = 1). This indicates that relationship of these copollutants and their error with the main pollutant does have an impact on the bias and standard error of main pollutant RR estimates. Type I error estimates in this study show the likelihood of detecting a false positive effect for a copollutant with an assumed null effect (i.e., RR = 1) is often reduced when measurement error is not present, especially when CO, NO x , and EC are the main pollutants. It is likely that this impact is more evident for CO, NO x , and EC than it is for PM 2.5 , SO 4 , and O 3 because the former are dominated by local sources, thus tend to be more spatially variable [17] and have higher degrees of correlated exposure measurement error. This can be important when trying to understand the independent effect of a pollutant that is a part of a complex mixture. For example, the EPA's most recent Integrated Science Assessments (ISAs) for CO [23] and NO x [24] both state that while there is evidence of consistent positive associations between short-term exposure to CO and NO 2 and effects on the respiratory system in epidemiologic studies, the challenge remains to disentangle the independent effect of CO or NO 2 related health effects from the larger air pollution mixture. There is the possibility that CO and/or NO x are serving as indicators of combustionrelated emissions, particularly from traffic, for some health outcomes. The moderate to high correlations between CO, NO x , and other pollutants generated from combustion processes (e.g., EC, PM 2.5 ) further complicates interpretation of epidemiologic studies of these pollutants. This provides motivation for the use of exposure metrics which minimize measurement error in epidemiologic analyses.
While the simulation presented here has many strengths, including accounting for multiplicative bias and correlated measurement errors, and the use of empirical data, there are inherent limitations. First, we assume that the refined exposure metric in each of our three constructed measurement error types provides an accurate representation of measurement error. That is, we assume that the AQ exposure metric accurately and completely accounts for spatial variability ignored in the CS metric, and we assume the PE exposure metric appropriately accounts for population variability ignored in the AQ metric. The authors acknowledge that the AQ metric and PE metrics themselves include some error, but in the absence of a 'true' spatially refined metric (e.g., a fine-scale measurement/monitor network in the study area) and a set of population-based personal exposure measurements, we operate under the above assumptions. Any error present in the refined metrics will remain an influence on the calculated RRs. Further, though it is common to use an individual's exposure to ambient-generated pollution in an epidemiologic study, the inclusion of exposure from indoor sources would allow for analysis of the impact of exclusion of this exposure source on commonly calculated RRs. Though we expect results presented here to have general applicability regarding the impact on bias of RRs in a copollutant time-series model, we acknowledge that results presented here are dependent on empirical data from the Atlanta metropolitan area. Multipollutant relationships may be different in other cities depending on the city-specific source profile and the magnitude of measurement errors may depend on the location of central monitor and population characteristics.
Simulation studies have examined the impact of measurement error in hypothetical multipollutant epidemiologic studies [16,25]. Exposure measurement error is often categorized into two classes: Berkson error and classical error, which are often assumed to be independent of the true exposure [16]. In reality, exposure measurement error is a combination of the two. Adding to the complexity, the consequences of each error type are different [5,26]. In most time-series studies, Berkson error will not bias the effect estimates, but will tend to increase the variance of the regression coefficients (increasing the width of the CIs and decreasing power), while classical measurement error tends to bias the true effect towards the null, with the magnitude of the effect attenuation depending on the error variance of the exposure estimate relative to the variance of the true exposure [16]. Empirical analyses have shown impacts of both types of error in single pollutant models [14,[27][28][29]. When two pollutants are measured with error, the correlation between the pollutants, and the correlation between their errors, predicts the magnitude of the bias, while the sign of the correlation between the pollutants predicts the direction of the bias [16,20]. In this analysis, we do not follow standard classical or Berkson measurement error models, but utilize the empirical relationships between different exposure metrics from a real-world epidemiologic study. By assuming that the more refined exposure metric is the true exposure while treating the less refined metric as the error-prone exposure, we observe effect attenuation and bias away from the null, depending on the copollutant pair considered. In addition, we note that time-series analyses typically rely on a small number of monitors to derive population-averaged exposure metrics. Recently there has been increasing interest in conducting spatiotemporal modeling of air pollutants via land-use regression models [30] and data fusion methods [31] to capture spatial heterogeneity within the study region. However, measurement errors are also associated with spatial predictions due to spatial smoothing and estimation uncertainty [21,32].

Conclusions
To the best knowledge of the authors, this study is the first to quantify the impact of measurement error on health risk estimates in copollutant models using empirical data. In addition, this study directly addresses a