The relationship between daily cardiovascular mortality and daily ambient concentrations of particulate pollutants (sulfur, arsenic, selenium, and mercury) and daily source contributions from coal power plants and smelters (individually, combined, and with interaction) in Phoenix, AZ, 1995–1998: A multipollutant approach to acute, time-series air pollution epidemiology: I

Abstract

The objective of this paper is to estimate the increase in risk of daily cardiovascular mortality due to an increase in the daily ambient concentration of the individual particulate pollutants sulfur (S), arsenic (As), selenium (Se), and mercury (Hg) using single-pollutant models (SPMs) and to compare this risk to the combined increase in risk due to an increase in all four pollutants by including all four pollutants in the same model (multipollutant model, MPM) and to the risks from source contributions from power plants and smelters. Individual betas in a multipollutant model (MPM) were summed to give a combined beta. Interaction was investigated with a pollutant product term. SPMs (controlling for time trends, temperature, and relative humidity), for an interquartile range (IQR) increase in the pollutant concentration on lag day 0, gave these percent excess risks (±95% confidence levels): S, 6.9% (1.3–12%); As, 2.9% (0.4–5.5%); Se, 1.4% (–1.7 to 4.6); Hg, 9.6% (4.8–14.6%). The SPM beta for S (as sulfate) was higher than found in other studies. The SPM beta for Hg gave the largest t-statistic and beta per unit mass of any pollutant studied. An (IQR) increase in all four pollutants gave an excess risk of 15.4% (7.5–23.8%), slightly smaller than the combination of S and Hg, 16.7% (9.1–24.9%). The combined beta was 71% of the sum of the four individual SPM betas, indicating a reduction in confounding among pollutants in the combined model. As and Se were shown to be noncausal; their SPM betas could be explained as confounding by S.

Implications: The combined effect of several pollutants can be estimated by including the appropriate pollutants in the same statistical model, summing their individual betas to give a combined beta, and using a variance–covariance matrix to obtain the standard error. This approach identifies and reduces confounding among the species in the multipollutant model and can be used to identify confounded species that have no independent relationship with mortality. The effect of several pollutants acting together may be higher than that of one pollutant. Further work is needed to understand the strong relationship of mortality with particulate mercury and sulfate.

Introduction

Over the past decade there has been an effort to move air quality management from a single-pollutant approach to a “multipollutant” approach (National Research Council, 2004; Hidy, 2010; Chow, 2010; Mauderly et al., 2010). This has been accompanied by an effort to develop an air pollution research program suitable to provide a sound basis for multipollutant air quality management (Vedal, 2011; Dominici, 2010; Greenbaum and Shaikh, 2010). It has been fairly common to include two pollutants in the same statistical model and compare the individual betas and confidence intervals in an attempt to identify the more toxic pollutant or the more likely causal pollutant (Burnett, 2004; Roberts, 2006; Brook et al., 2007). However, it is now generally recognized that this approach is not reliable (Vedal, 2011; Tolbert et al., 2007; Ito et al., 2007). Ito et al. (2007) argue that unless the pollutant variables are uncorrelated, pollutants in a multipollutant model (MPM) cannot be considered as independent variables. Therefore, multipollution approaches to air pollution epidemiology that use betas or t-statistics (or derived confidence intervals) as an indicator of the relative toxicity of the individual pollutants may provide indications but are not necessarily reliable (Burnett, 2004; Roberts, 2006; Brook et al., 2007). While the individual betas in an MPM may not be considered representative of the effects of the individual pollutants, the sum of the betas (combined beta) may be considered representative of the combined effect of the multiple pollutants. Thus, by putting several pollutants in the same statistical model, we can obtain the combined health effect of the pollutants in the model. Thurston et al. (2013) have used a similar approach, combining many pollutants in an MPM to obtain a total risk index for long term exposure to air pollutants. This paper is apparently the first to apply this multipollutant approach to acute, time-series air pollution epidemiology. Source contributions obtained from source apportionment models also provide an approach to investigating the combined effects of multicomponent mixtures.

Methods

Air pollution data

Chemical composition data were obtained from the U.S. Environmental Protection Agency (EPA) National Exposure Research Laboratory (NERL) research monitoring site in central Phoenix. The site and measurements program are described in detail in Mar et al. (2000, 2003) and in a EPA report (EPA, 1998). In brief, NERL investigators collected daily, gravimetric, integrated 24-hr (starting at 0700), fine particle samples (PM_2.5) on Teflon filters. The concentrations of x-ray fluorescence (XRF) elements (including S, Hg, As, and Se) were determined at the EPA (Research Triangle Park, NC) with energy-dispersive x-ray fluorescence. Temperature and relative humidity measurements (1238 days from Feburary 9, 1995, until June 30, 1998) were obtained from measurements at the monitoring site, or in case of missing data, from the National Oceanic and Atmospheric Administration weather station at the Phoenix (AZ) Sky Harbor airport. Concentrations of XRF elements were available for 993 days. Pearson correlation coefficients between pollutants, temperature, and relative humidity are shown in Table 1. Except for the correlations between Se and As (r = 0.687), As and S (r = 0.400), and Se and S (r = 0.220), all correlations between pollutants were weak. Since the species differed greatly in concentrations, nanograms and micrograms per cubic centimeter (ng/cm³ and µg/cm³) concentration units were converted to interquartile range concentration units (IQR) prior to analysis. IQR values are given in Table 2. Betas are calculated in terms of IQR units; one IQR unit equals the IQR concentration. Plots of the pollutants versus time are provided in the supplemental data. Plots of S, As, and Se are similar, with higher concentrations during August to December. Hg concentrations did not show any seasonal variation.

Table 1. Pearson’s correlation coefficient, r

r	Temp	RH	Se	As	Hg
RH	0.540	1
Se	–0.279	0.322	1
As	0.210	0.254	0.687	1
Hg	0.084	0.006	–0.028	–0.066	1
S	0.476	0.017	0.220	0.400	0.024

Table 2. Relationship of cardiovascular mortality with single and multiple pollutants (for an IQR increase in pollutant concentration)

	IQR^a	Beta^b	SE^c	t^d	%ER^e	%LCI^f	%UCL^g
Hg	1.1	0.091	0.023	3.99	9.60	4.77	14.64
S as Sulfate	844	0.066	0.027	2.45	6.86	1.34	12.67
As	1.7	0.028	0.013	2.30	2.93	0.42	5.50
Se	0.5	0.014	0.016	0.88	1.41	−1.70	4.63
	Sum Betas^h	Comb Beta^i	Comb^j SE	Comb t	Comb %ER	Comb %LCI	Comb %UCL
Hg + S + As + Se	0.201	0.143	0.036	3.98	15.38	7.53	23.80
Hg + S	0.158	0.154	0.034	4.48	16.69	9.06	24.86

Notes: ^aIQR = interquartile concentration in ng/m³ (IQR =75%tile – 25%tile).

^bBeta = the change in the natural logarithm of mortality for one IQR unit change in the pollutant concentration [β = ln(rate ratio)].

^cSE = standard error of beta from GAM-Exact.

^dt = t-statistic (beta/SE).

^e%ER = Increase in relative risk as a percent of the baseline risk.

^f%UCI = upper 95% confidence limit.

^g%LCI = lower 95% confidence limit.

^hSum Betas = sum of betas from single-pollutant models.

ⁱComb Beta = sum of individual betas in a co- or multipollutant model.

^jComb (parameter) = the value of that parameter for the combined beta.

Mortality data

Mortality data for all of Maricopa County from January 1, 1995, to June 30, 1998, were obtained from the Arizona Center for Health Statistics in Phoenix, AZ. Death certificate data included residence ZIP Code and the primary cause of death as identified by the International Classification of Diseases, Ninth Revision (ICD-9; World Health Organization, Geneva, Switzerland). In this study, relationships were investigated with cardiovascular mortality (CVM) only (ICD-9 codes 390–448.9) for persons 65 years and older. For the time period covered in the analysis, there were 4847 cardiovascular deaths or an average of 3.91 deaths per day. Autocorrelation in the natural logarithm of the daily death count was low; R² never exceeded 0.01 out to lag 10.

Study area

Only the deaths of residents in ZIP Codes located near the monitoring site were included in this study. This ZIP Code region, hereafter called Central Phoenix, was recommended by the Arizona Department of Environmental Quality (Phoenix, AZ) as being a area in which the average community pollutant concentrations would be well represented by the concentrations measured at the monitoring site (Mar et al., 2000, Wilson et al., 2007). A map of Phoenix showing the ZIP Code areas used in this study and the location of the monitoring site may be found in Wilson et al. (2007). Maps and additional information on the socioeconomic status of each ZIP Code are given in the supplemental data.

Statistical model

Mortality outcomes were modeled using Poisson log-link regression allowing for overdispersion in a generalized additive model (GAM) (Dominici et al., 2004) with a stringent convergence criterion (10e-10 with 1000 iteration steps). Mortality count data are distributed lognormally, so transformation from counts to the natural logarithm was appropriate. Autocorrelation in the residuals was checked for S and found to be low (R² < 0.007) out to lag 10. The distribution of residuals appeared to be adequately random (plot shown in supplemental material). Analyses were conducted using R (R Core Team, 2012). The GAM-Exact technique, described by Dominici et al. (2004), was used to obtain the final betas and the standard errors in betas. The base model used smoothing splines to control for time trends, temperature, and relative humidity. The number of degrees of freedom (df) for the time trend (df = 12) was selected by averaging the df obtained by identifying the minimum Akaike information criterion (AIC) (df = 15) and the df obtained by identifying the minimum absolute sum of the partial autocorrelation function (PACF) (df = 9). The potential confounding effects of weather were accounted for by including smooths for temperature (lag 2 days, df = 5) and relative humidity (lag 0 days, df = 2); lag days and df were based on minimization of the AIC. Indicator variables were used to adjust for day of the week. An extensive analysis of the sensitivity of the results for particulate matter (PM) to temperature and relative humidity controls is given in Mar et al. (2000). In the earlier Phoenix studies (Mar et al., 2000, 2003), missing values were replaced with the average of all measured values. For this study, the Source Apportionment Workshop (Mar et al., 2006), and the more recent Phoenix study (Wilson et al., 2007), days with missing values were omitted from the analyses.

Single-pollutant models (SPMs) had only one pollutant variable in the model. In the multipollutant approach, variables representing two or more pollutants were included in the same statistical model. The individual betas calculated by the model were added to give a combined beta. The standard error for the combined beta was estimated using a variance/covariance matrix (T. Lumley, personal communication, 2008). Since the analysis used IQR units, using the combined beta in the equation [β = ln(rate ratio)] yields the increase in risk for one IQR increase in concentration of each pollutant in the model. To convert to gravimetric units, divide beta by the IQR (in ng/m³) and multiply by the number of units; that is, [(β × 10)/IQR = ln(rate ratio)] gives the relative risk for a 10-ng/m³ increase in the pollutant. With an SPM, we may think of the increase in risk being for a daily or a yearly average increase in the pollutant. For a combined beta, from an MPM, we must think in terms of a yearly increase in risk for an increase in the yearly average of each pollutant in the model.

Results and Discussion

Single-pollutant models

The beta, standard error, t-statistic, probability, percent increase in risk for an IQR increase in the pollutant concentration, and ±95% confidence intervals (CI) for the relationship of each pollutant with CVM are given in Table 2 for lag day 0, the only lag day out to lag day 10 for which the relationship between mortality and pollutant was statistically significant.

Sulfate

In this study, the percent excess risk (95% confidence intervals) for an increase of 1 µg/m³ of sulfate (measured as XRF-S) was 8.0% (1.4–15%), higher than reported in most previous studies (EPA, 2009). Most studies of air pollution epidemiology use county-level mortality data. The use of ZIP-Code-level mortality data makes it possible to focus on a population in an area where the ambient concentration is well represented by the monitor. The exposure error for this population should be lower than the exposure error for the entire county (Mar et al., 2000). For Maricopa County, the value for sulfate is 1.8% (–1.6 to 5.2%), more in line with other results for county-level mortality (EPA, 2009). In Phoenix the selected ZIP-Code population also has a lower and more homogeneous socioeconomic status (SES) than the larger county population. Further work is needed to determine which population property is most important, a lower exposure error or a more homogeneous and lower SES.

Arsenic and selenium

As and Se were included in the NPACT study (Ito et al., 2013) and Se has been used as an indicator of coal combustion (Laden et al., 2000).

Mercury

The extremely strong relationship of CVM with particulate Hg was unexpected and surprising. The neurotoxic effects of Hg are well known. However, there is some evidence from toxicology to indicate that Hg may also have cardiovascular effects. Mozaffarian (2009) lists the following “Experimentally-observed effects of mercury which may increase CVD risk” (p. 1899):

Systemic Effects

Promotion of free radicals and reactive oxygen species

Inhibition of antioxidant systems (glutathione peroxidase, catalase, superoxide dismutase)

Increased lipid (e.g., LDL cholesterol) peroxidation

Promotion of blood coagulation (clotting)

Inhibition of endothelial cell migration

Direct Cardiovascular Effects

Reduction in myocardial contractile force

Increased calcium release from myocardial sarcoplasmic reticulum

Reduction in left ventricular myosin ATPase activity

Decreased heart rate variability and increased blood pressure

Hence, it would not be possible, on the basis of toxicity alone, to rule out the possibility that particulate Hg, like other components of air pollution, may increase CVM.

However, it is not unreasonable to ask whether the relationship of CVM with Hg is real or spurious. There are several reasons to suspect that it may be spurious. Many of the values are negative, and the average is –0.13 ng/m³. However, this may reflect too large a correction for the filter blank or the XRF background. The confidence limit is about the same size as the IQR. Hg does not show any seasonal variation. Hg has an unusually large standard deviation for an air pollutant. Perhaps Hg comes from a local source that impacts the sampling site on some days and not on others. The large percent excess risk of 9.6% (4.8–14.6%) per increase in one IQR of 1.1 ng/m³, if true, would mean that for the same increase in mass, Hg would be 1000 times as toxic as sulfate. This suggests that Hg may be a surrogate for some other pollutant emitted from the same source in much larger quantities, perhaps total Hg (particulate + gaseous + semivolatile). (The semivolatile Hg may have evaporated from the filter by the time the filter was analyzed.) However, there is one reason to feel that Hg deserves further study. Hg’s t-statistic value of 3.99 gives a p of < 10⁻⁴, or, statistically speaking, there is less than 1 chance in 10,000 that the relationship is spurious.

As far as we are aware, this is the first study to look for a relationship between the concentrations of particulate Hg collected on a filter and CVM. Although the Hg concentration time series is not of the highest quality, the high beta and t-statistic suggest that the relationship might be real and that additional epidemiologic studies using particulate Hg would be useful. State monitoring agencies should report particulate Hg levels (readily available from the XRF analyses used to determine metal concentrations) to the EPA monitoring database. Studies of particulate Hg relationship with mortality should include populations near monitoring sites and populations with different SES in order to investigate exposure error and SES effects. In addition, some studies should be undertaken to obtain better measurements of both particulate and gaseous Hg over sufficiently long time periods to provide data suitable for epidemiologic analyses.

Combined model with all four pollutants

We would like to know the combined effect of S, As, Se, and Hg. However, we cannot just add the SPM betas because each SPM beta will represent not only the effect of the single pollutant in the model but also, due to confounding, part of the effect of other correlated pollutants. Including several pollutants in an MPM prevents confounding among the pollutants in the model. However, the combined beta may be confounded by other pollutants, not included in the model, that have an effect independent of the pollutants in the model and are correlated with one or more of the pollutants in the model. Thus, the combined beta may overestimate the effect of the pollutants in the model, but, if so, it will underestimate the effect of all pollutants present in the atmosphere. Table 2 shows the results for a combined model containing S, As, Se, and Hg. The combined beta is almost identical to the combined beta for S and Hg. This suggests that the SPM relationship of As and Se with CVM can be attributed to their correlation with S. The combined beta for the four pollutants is only 71% of the sum of the betas from the four SPMs. Thus, the use of the combined model reduces confounding among the pollutants in the combined model and gives a more accurate estimate of the combined effects of the four pollutants.

Causal, or correlated but noncausal; confounder or confuser?

Confounding presents a major problem in interpreting the results of epidemiologic analyses. We would like to know the combined effect of S, As, Se, and Hg. However, because of confounding, the sum of the SPM betas would misrepresent the result. An extensive literature treats confounding in clinical studies, for example, Steyer and Schmitt (1994 and references therein). Tests have been developed for using the SPM and co-pollutant model (CPM) to identify confounding in clinical epidemiology. First will describe these tests and show that they are not applicable to acute, time-series epidemiology (ATSE). Then we develop some tests to identify and quantify confounding in ATSE. More importantly, we show how to identify species that are noncausal but are correlated with a second pollutant. If the second pollutant is causal, or has an independent relationship with mortality, the noncausal pollutant, when included in a CPM, will change the relationship of second pollutant with mortality. We call such species confusers.

The test developed in Steyer and Schmitt (1994) uses the following notation:

(1)

(2)

The test was: If α₁ β₂ = 0, not confounded; if α₁ β₂ ≠ 0, confounded. The test was applied to examples from clinical epidemiology. Thus, if α₁ = 0, or near 0, X has no effect and will not confound; if β₂ = 0, or near 0, W has no effect and X will not confound W. If α₁ β₂ ≠ 0, there will be confounding. This test does not account for what happens when X and W are correlated. In ATSE it would be rare to have an α or β near zero. Thus, this test is not adequate for ATSE.

Quantitative criteria are needed for ATSE. Lecture notes from a Boston University (2015) online course give a different but quantitative test, also for clinical epidemiology.

A simple, direct way to determine whether a given risk factor caused confounding is to compare the estimated measure of association before and after adjusting for confounding. In other words, compute the measure of association both before and after adjusting for a potential confounding factor. If the difference between the two measures of association is 10% or more, then confounding was present. If it is less than 10%, then there was little, if any, confounding.

Using the preceding notation, the criterion is given as (α₁ − β₂)/β₂. This test is also illustrated with examples from clinical epidemiology. This test would work if X and W had a 0 correlation with each other. But in that case the correlation coefficient would provide the test. It might also work if X and W were correlated and had the same correlation with their respective dose to the body. However, this is generally not the case for air pollutants. Thus, this test is also not adequate for ATSE.

Perhaps it will be useful here to recall that the definition of confounding is different for clinical epidemiology and ATSE.

Conditions necessary for confounding in clinical studies

1. The confounding factor must be associated with both the risk factor of interest and the outcome.

2. The confounding factor must be distributed unequally among the groups being compared.

3. A confounder cannot be an intermediary step in the causal pathway from the exposure of interest to the outcome of interest.

Conditions necessary for confounding in ATSE

1. The confounding variable must be associated (inversely or directly) with the outcome of interest.

2. The confounding variable must be associated (inversely or directly) with the outcome of interest, independent of the variable of interest.

3. The confounding variable must be associated (inversely or directly) with the variable of interest.

4. The confounding variable is not in the causal pathway of the variable of interest to the outcome of interest.

As stated in condition 2, a confounder must have an association with the variable of interest and also with the outcome. In the following section, we develop a logical foundation for identifying and quantifying confounding in ATSE. We also show how to identify a variable that is not associated with the outcome and therefore is noncausal, but is associated with the variable of interest and therefore can lead to confusion.

Let us consider an ideal situation in which a population is exposed to only two pollutants with daily concentrations of P₁ and P₂ (given in the same units) and we have a time series of the daily pollutant concentrations and the number of deaths per day, Y. We assume a linear dependence of ln Y on P₁ and P₂ and express the expected values E(ln Y│P₁) of variable ln Y as a function of P₁ for a SPM or E(ln Y│P₁, P₂) for a CPM. In order to simplify notation, we follow the example of Steyer and Schmitt (1994) and use α for the coefficient of pollutant concentration in SPMs and reserve β for CPMs and MPMs. Thus, we have for our ideal situation

(3)

(4)

If P₁ and P₂ are correlated, that is, R² (P₁│P₂) does not equal zero, P₁ and P₂ will confound each other. Part of the relationship of P₁ with ln Y will be transferred to P₂ and part of the relationship of P₂ with ln Y will be transferred to P_1. The sign of the quantity transferred is given by the sign, s, of the Pearson’s correlation coefficient, r. The amount of relationship depends on R² (P₁│P₂). Thus, we have for single-pollutant models

(5)

(6)

where α_x^T is the true, or unconfounded alpha.

For CPMs or MPMs we have

(7)

It will be convenient to define α^S as the sum of the α terms from the SPMs and β^C as the sum of the β terms from the CPM, which we refer to as the combined β, as

(8)

(9)

We can consider several cases. Since we can rarely prove causality with a single ATSE study, instead of referring to pollutants as causal, or possibly causal, we merely assume that the pollutant has a relationship with mortality that is independent of the relationship of the other pollutants in the model. We call such species independently related to mortality (IRM)—that is, pollutant P₁ is IRM. An IRM pollutant could be causal or it could be noncausal but highly correlated with another pollutant not included in the model.

Case 1

P₁ and P₂ are both IRM but uncorrelated, that is, R² (P₁│P₂) = 0.

P₁ will not confound P₂ and P₁ will not confound P₁.

(10)

(11)

(12)

This important case occurs when a source apportionment model such as factor analysis or principal component analysis is used to obtain source contributions that are uncorrelated.

Case 2

P₁ and P₂ are correlated and both P₁ and P₂ are IRM. In SPMs, P₁ and P₂ are said to confound each other. We assume that mortality is not caused by ambient concentrations of pollution but by the dose of pollutant to the lung. Different pollutants may differ widely in the correlation between ambient concentration and dose to the lung. In our ideal situation, however, P₁ and P₂ both have the same degree of correlation with the amount of P₁ and P₂ deposited in the lung.

In the combined model, the statistical analysis will remove the confounding so that

(13)

(14)

(15)

(16)

If r is positive, β_X will be less than α_X. If r is negative, β_X will be more than α_X. We can use ((α_X – β_X)/β_X) × 100 as a criterion to give some idea of how much confounding contributes to α_X. Call this Criterion C₁.

Case 3

This is like Case 2, but P₁ and P₂ have different degrees of correlation with their respective dose. Generally, one pollutant will be better correlated with the dose than the other. The statistical analysis will distribute more of the effect to the pollutant that has the better correlation with the dose. Therefore,

(17)

However, the sum of the β terms will still equal the sum of the true α^T terms, that is,

(18)

Therefore, β^C, the sum of the β terms in the CPM, will give the true effect on mortality of a unit increase in both pollutants.

It is not uncommon for analysts to control for confounding by P₂ by putting P₂ in a CPM model with P₁ and either assuming that β₁ will equal α₁^T or assuming that β₁ and β₂ will give the relative toxicities of P₁ and P₂. However, this will only be true in the unlikely situation described earlier when P₁ and P₂ have the same correlation with their dose.

For this case we cannot use C₁, but we can use C₂, ((α^S – β^C) × 100/β^C), as a criterion for confounding or confusing. This gives the difference between the combined beta, β^C, from the CPM, and the sum of the alpha’s from the SPMs. If C₂ is positive we have positive confounding (r is positive) and the individual α terms are too large; if C₂ is negative we have negative confounding (r is negative) and the individual α terms are too small. The individual β terms cannot be used in C₁ because their values will depend on their correlation with each other and the correlation of the pollutant’s concentration with its dose.

Case 4

P₁ and P₂ are correlated but only P₁ is IRM, that is, α₂^T = 0. α₂ will be positive because P₂ will be confounded by P₁.

(17)

In this case, α₂ is due solely to the correlation of P₁ with P₂ and P₂ has no independent relationship with ln Y. In the CPM, all the relationship due to P₁ will be captured in β₁, and hence β₂ will equal 0.

(18)

(19)

(20)

(21)

In the ideal case, β^C would equal α₁ + α_{2 .} In a nonideal case, relationships will not be as exact but the criterion C₃ =100 × (β^C – α₁₎/ β^C) will indicate how much β^C has been increased and provide information regarding how likely P₂ is to be causal. If C₃ is near zero or negative, P₂ is likely noncausal.

Case 5

P₁ and P₂ are correlated and P₁ is IRM but P₂ is noncausal. This is similar to Case 4. However, in this case, P₁ is a multicomponent mixture and P₂ is highly correlated with the causal component of P₁. For example, in the case of S and As, As is highly correlated with the SO₂ emitted from power plants and therefore likely to be highly correlated with the toxic metals neutralized and made soluble by reaction with sulfuric acid generated from the SO₂ emissions. In this case, since both P₁ and P₂ are correlated with the toxic component of P₁, the statistical analysis will distribute a portion of the effect to both P₁ and P₂. Since including P₂ causes β₁ to be less than α₁, as happened in Case 4, we might think that P₂ is a confounder of P₁. However, since P₂ is noncausal, it would not be considered a confounder so we call it a confuser. Since including P₂ in the model will not increase the overall effect given by β^C, β^C will not be significantly greater than α₁ and in the ideal case β^C would equal α₁. Again criterion C₃ provides a useful test of whether P₂ is causal or correlated, confounder or confuser. If the difference between β^C and α₁ relative to α₁, given by C₃ is very small, we can consider P₂ to be noncausal and a confuser. Since C₃ equal –6.26% for the S + As CPM, we can identify As as noncausal since it has no relationship with mortality except that due to its correlation with S. Again we cannot use C₁ since the βs depend on the correlation of P₁ and P₂ with the causal component of P₁.

Case 6

Let us add a third pollutant, P₃. P₂ and P₃ are correlated with each other but not with P₁. P₁ and P₂ are causal but P₃ is not. Then the SPM for P₃ is

(22)

(23)

Or since α₃^T = 0,

(24)

In a CPM with P₁, P₃ will appear to be causal and a confounder. However, including P₃ in a CPM with P₂ will show (Case 5) that P₃ is noncausal and a confuser. The relationship between P₃ and mortality is due to its correlation with the causal species P₂.

Relationships in the real world are not as clean. We need to remember that

(25)

where the summation is over all P_X’s known and measured, known but not measured, and unknown. (For P_X= P_1, s is positive and R² = 1, so the first term in the summation is α₁^T∙P₁.) Therefore, while the relationships we discuss are useful, they may not be exact because of the many additional possibilities for confounding.

Co- and multipollutant models

We are now ready to interpret the results of the CPMs and MPMs given in Table 3.

Table 3. Analysis of confounding and identification of noncausal species

	α	β	t-Statistic		C1^a (as %)	C2^b (as %)	C3^c (as %)
	SPMs	MPM	SPMs	MPM	α^x – β^x/βx	α^s – β^c/β^c	β^c – α^1/β^c
r = 0.024
Hg	0.0916	0.0903	3.99	3.93	1.47%
S	0.0663	0.0641	2.45	2.39	3.46%		40.65%
Hg + S	0.1579	0.1544	NA	4.48		2.30%
r = 0.400
S	0.0663	0.0435	2.45	1.37	52.42%
As	0.0289	0.0189	2.30	1.27	52.90%		–6.26%
S + As	0.0952	0.0624	NA	2.32		52.57%
r = 0.220
S	0.0663	0.0653	2.45	2.28	1.54%		1.01%
Se	0.0140	0.0017	0.88	0.10	734.51
S + Se	0.0803	0.0670	NA	2.39		19.91%
r = –0.066
Hg	0.0916	0.0941	3.99	4.09	−2.68
As	0.0289	0.0312	2.30	2.46	−7.29		26.88%
Hg + As	0.1205	0.1253	NA	4.74		–3.82%
r = –0.028
Hg	0.0916	0.0915	3.99	3.98	0.15%
Se	0.0140	0.0138	0.88	0.86	1.55%		12.97%
Hg + Se	0.1056	0.1053	NA	3.79		0.34%
S	0.0663	0.0434	2.45	1.38	52.78%
As	0.0289	0.0270	2.30	1.50	7.09%
Se	0.0140	–0.0161	0.88	–0.79	–186.97%		–22.18%
S + As + Se	0.1092	0.0543	NA	1.89		101.24%
Hg	0.0916	0.0927	3.99	4.03	–1.15%
S	0.0663	0.0365	2.45	1.16	81.51%
As	0.0289	0.0227	2.30	1.55	27.07%		–1.59%
Hg + S + As	0.1868	0.1520	NA	4.38		22.94%
Hg	0.0916	0.0903	3.99	3.95	1.47%
S	0.0663	0.0629	2.45	2.23	5.38%
Se	0.0140	0.0020	0.88	0.12	615.71%		0.52%
Hg + S + Se	0.1719	0.1552	NA	4.36		10.80%

Notes: ^aC1 indicates the amount of confounding of pollutant X in a SPM by pollutant Y, i.e., first entry, Hg confounded by S; second entry, S confounded by Hg.

^bC2 indicates the amount of confounding of pollutants 1 and 2 by each other in SPMs.

^cC3 indicates the likelihood that P1 is causal or noncausal. If, C3 is near zero or negative, P2 is likely noncausal.

Hg + S: An example of Case 1 with possible contributions from Case 2 and 3. Confounding is low, 1.47% for α₁ (Hg), 3.46% for α₂ (S), and 2.30% for β^C. The contribution of S to β^C is 40.65%, indicating that S is likely to be IRM. R² is small and r is positive, so confounding slightly increases α₁ over α₁^T, α₂ over α₂^T, and α^S over β^C.

S + As: Example of Case 5. Confounding is high, >50% for α₁ (S), α₂ (As), and β^C (S + As). However, as discussed under Case 5 earlier, the contribution of As to β^C is a minus 6.26%, indicating that As is unlikely to be causal. The reduction in β^C occurs because the As time series provides only error, that is, it is all noise with no signal, so it biases the β^C to a lower value.

Hg + As: Example of Case 6. Since r is negative, β₁ will be larger than α₁ (Hg), and β₂ will be larger than α₂ (As). The contribution of As to β^C is 26.88%, suggesting that As is likely to be IRM. However, we know from our analysis of S + As that As is not likely to be causal and that its contribution to β^C is due to its correlation with S.

Hg + S + As: Extension of Case 5 to an MPM. In this three-pollutant model, the addition of As to Hg + S adds only 3.93% to β^C, confirming that As is noncausal. Therefore, we would consider that As is a confuser, not a true confounder. The confounding of α₂ (S), α₃ (As), and β^C (Hg + S + As) β^C is high due to confounding of S and As by each other. The confounding of α₁ (Hg) (–1.15%) is a combination of positive confounding by S and negative confounding by As. For testing Criterion 3 in a MPM, we would use β₁ + β₂ instead of α₁.

The Steyer and Schmitt (1994) test would indicate confounding in all these examples, incorrectly in the case of Hg + S. The Boston University (2015) test would correctly indicate nonconfounding for Hg and S and, in the MPM case, correctly indicate that S and As were confounded but Hg was not. It would not work for the S + As example or another situation in which P₁ and P₂ were correlated. However, neither test would provide any indication that As and Se were likely noncausal.

Is sulfate causal?

It is unlikely that the sulfate anion (SO₄^2-) is the biologic agent causing the toxicity associated with exposure to sulfate (Lippman et al., 1996). It is more likely that the cations (positive ions) associated with the sulfate anion (SO₄^2-), determine the toxicity of sulfate containing particles. Sulfuric acid, formed from atmospheric oxidation of SO₂, can react with basic metal oxides of toxic metals to convert insoluble oxides into soluble ions (Cd, Fe, Hg, Ni, Pb, V, Zn). These heavy metal ions are thought to be toxic (Amdur et al., 1978; Costa and Dreher, 1997; Ghio et al, 1999). Sulfate may also exist in the atmosphere as sulfuric acid aerosol (H₂SO₄) or, if neutralized by ammonia, as (NH4)₂SO₄ particles. Both types of particles yield potentially toxic acidic solutions in water (Lippmann et al., 2000; Gwynn et al., 2000).

In a review of the available information on the epidemiology and toxicology of air pollutants, the British Committee on the Medical Effects of Air Pollutants (COMEAP, 2009) concluded that

Sulphates could play a role in the capacity of PM_2.5 to drive pulmonary inflammation via either increasing metal mobilization from PM, and hence increasing the ability to drive oxidative stress leading to inflammation, or decreasing the pH of lung-lining fluid, leading to decreased pathogen clearance and hence increased inflammation. Inflammation in the lung could then initiate changes in blood clotting, and/or activation of macrophages in atherosclerotic plaques increasing their instability.

While sulfate is likely not the biologically causal agent, sulfuric acid is in the causal pathway, for example:

Which link to call causal is perhaps more of a semantic than a scientific question. The key concern, however, is: “Will control of SO₂ from power plants, and hence a reduction of sulfate in PM_2.5, lead to a reduction in CVM?” This study, and many other studies that find an independent relationship of sulfate with CVM, suggest that the answer to the question is yes.

Source contributions

Sulfate is actually a multicomponent pollutant, having several fine particulate sources (power plants [oil and coal], smelters, possibly traffic [S in gasoline]) as well as coarse particulate sources (CaSO₄ from wallboard and possibly other sulfates in resuspended soil). Since Hg is not correlated with sulfate and As and Se have been shown to be noncausal, we are left with sulfate as the only measured pollutant available to represent pollution from both coal power plants and smelters.

Source apportionment techniques may help separate sulfate into its various sources. Thurston et al. (2011) have suggested that sulfate, or a group of associated species, may serve as a surrogate for particles from coal combustion. Using factor analysis in a study of 64 cities, they find a coal combustion factor including contributions from SO₂, As, Se, and Hg. They did not include sulfate in their source apportionment analysis. It is likely that Hg, Se, and As are not the PM_2.5 components causing the health effects associated with sulfate, but instead are serving, either individually or as a source signature, as a marker for other gaseous or particulate by-products of coal combustion.

Several types of source apportionment models (SAMs) have used the Phoenix pollutant database (Mar et al., 2000; Hopke et al., 2006; Mar et al., 2006). We next use factor analysis (FA), positive matrix factorization (PMF), and principal component analysis (PCA). In Phoenix, As, Se, and S are correlated but they are not correlated with SO₂ or Hg. FA yielded one factor for regional sulfate. PMF splits regional sulfate into one coal combustion source and one smelter source. PCA finds one coal combustion and two smelter sources. SAMs appear to reduce the contribution from coarse sulfate in soil in PM_2.5 as shown by the reduced correlation between of coal power plant and smelter sources and SOIL compared to the correlation between XRF-S and SOIL (Table 4).

Table 4. Pearson’s correlation coefficient, r, for source contributions with SOIL and S and for the coal power plant (PP) sources with the smelter sources

r	SOIL	S as SO4	RegSO4 FA	Coal PP PMF
S	0.263	1
RegSO4 (FA)	0.084	0.963	1
CoalPP (PMF)	0.067	0.872	0.854	1
CoalPP (PCA)	0.072	0.953	0.918	0.914
r-PMF	Coal PP	r-PCA	Coal PP	Smelter 1
Smelter	–0.209	Smelter 1	1.9 × 10^−8	1
		Smelter 2	–2.1 × 10^−8	2.0 × 10^−8

Note: SOIL is the sum of crustal element oxides minus nonsoil potassium.

We can also combine sources into the same model in order to estimate the health effects of coal power plants and smelters together. Table 5 gives results of SPM and MPM for these source contributions. It is interesting to see how the t-statistic and the percent excess risk both increase as we go from a single species sulfate (2.45 and 6.86%), to regional sulfate from FA (2.90 and 9.40%), to a coal power plant and a smelter from PMF (2.87 and 11.20%). Further splitting the smelter factor into two smelters (PCA) does not further improve the statistical significance (t-statistic) or increase the risk calculated for the combination of the two smelters and the coal power plant (PCA) (2.48 and 9.17%). This may be because PCA splits out Ni and V as having separate sources. Note that the negative correlation between the coal power plant and smelter 2 sources causes the β terms to be larger than the α terms. The increase in the t-statistic and the percent excess risk in going from XRF-S to sources may be due to removing the S in soil, as shown by the decrease in correlation in the source contributions with soil compared to that of XRF-S with SOIL (Table 4).

Table 5. Relationship of cardiovascular mortality with source contributions (for an IQR increase in the contribution of that source)

	α^a	β^a	t^b	t^b
	SPM	MPM	SPM	MPM	%ER^c	%LCI^d	%UCL^e
S as sulfate	0.066		2.45		6.86	1.34	12.67
FA RegSO4	0.090		2.90		9.40	16.24	2.96
PMF CoalPP	0.063	0.065	2.12	2.12	6.54	12.97	0.49
PMF Smelter	0.040	0.041	1.92	1.92	4.05	8.37	–0.09
PMF CoalPP + Smelter	0.103	0.106	NA	2.87	11.20	19.57	3.41
PCA CoalPP	0.055	0.052	2.04	1.92	5.62	11.31	0.21
PCA Smelter 1	0.004	0.005	0.22	0.26	0.41	4.14	–3.19
PCA Smelter 2	0.033	0.031	2.38	2.24	3.33	6.16	0.57
PCA CoalPP + Smelters 1 and 2	0.091	0.088	NA	2.48	9.17	17.02	1.86
PCA CoalPP + Smelter 2	0.087	0.083	NA	2.78	8.64	15.18	2.47

Notes: ^aα (SPM) or β (MPM) = the change in the natural logarithm of mortality for one IQR unit change in the source contribution [α or β = ln(rate ratio)].

^bt = t-statistic (beta/SE).

^c%ER = Increase in relative risk as a percent of the baseline risk.

^d%UCI = upper 95% confidence limit.

^e%LCI = lower 95% confidence limit.

It would appear that an MPM containing species from coal power plants and smelters or a combination of source contributions from coal power plants and smelters could be used to estimate the health effects from these sources. Source contributions have the advantage of separating the two types of sources, while MPMs have the advantage of providing information on which species are noncausal.

Interaction

Hg and S

To investigate a possible interaction involving Hg and S, we use a variable obtained by multiplying Hg times sulfate and transforming the product into IQR units. (Hg and S are also in IQR units.) Results are shown in Table 6. In a SPM, Hg*S has a positive, statistically significant (SS) beta. However, when Hg*S is included with Hg and S in a MPM, its beta is negative and not SS. However, the individual betas for Hg and S in the MPM are both greater than in an SPM or a CPM. The combined beta from all three variables is higher than that from the CPM. This indicates that when both Hg and S concentrations are high, neither is as effective as when only one is high. Apparently they interfere with each other in some way (antagonism?). When Hg*S is included in the MPM, it controls for this interaction and gives the correct betas for Hg and S if only one species was present in the atmosphere. AICs were calculated but were not different enough to suggest that one model was better than another.

Table 6. Interaction of Hg and S and power plant (PP) and smelter sources

	α^a	β^a	t^b	t^b
	SPM	CPM	SPM	CPM
Hg and S
Hg	0.092	0.13	3.99	2.61
S	0.066	0.133	2.45	1.65
Hg*S	0.087	–0.072	3.78	–0.89
Hg + S + Hg*S	0.245	0.19	NA	3.70
Hg + S	0.245	0.154	NA	4.48
PMF
PMF CoalPP	0.063	0.068	2.12	1.80
PMF smelter	0.040	0.043	1.92	1.54
PMF product	0.030	–0.003	2.07	–0.12
CoalPP + smelter + product	0.133	0.109	NA	2.60
PMF CoalPP + smelter	0.103	0.106	NA	2.87
PCA
PCA CoalPP	0.055	0.062	2.04	1.35
PCA Smelter2	0.033	0.038	2.38	1.39
PCAProduct	0.035	–0.009	2.59	–0.26
CoalPP + Smelter2 + product	0.122	0.091	NA	2.28
CoalPP + smelter	0.087	0.083	NA	2.78

Notes: ^aα (SPM) or ^aβ (MPM) = the change in the natural logarithm of mortality for one IQR unit change in the source contribution [α or β = ln (rate ratio)].

^bt = t-statistic (beta/SE).

Power plant and smelter sources

As shown in Table 6, the interactions between coal power plant sources and smelter sources for both the PMF and PMC SAMs give results very similar to those found for Hg and S. When both source contributions are high, neither is as effective as when only one is high.

Sensitivity Analysis

Additional analyses were conducted to investigate how changes in the statistical model might affect the resulting betas and t-statistics.

Controls for temperature and relative humidity

In the statistical model, controls for temperature and relative humidity were based on minimization of AIC. We calculated beta for the relationship of the various pollutants with CVM with no controls, with controls used in the model based on minimization of AIC, and with controls based on inclusion of temperature and relative humidity as explanatory variables in a distributed lag model, for lag day 0 and lag days 0 + 1, 0 + 1 + 2, and 0 + 1 + 2 + 3. The betas were not very sensitive to changes in the temperature and relative humidity controls, suggesting that confounding by temperature or relative humidity was not a problem. Results are shown in the supplemental data.

Temporal control

In the statistical model, variations in beta with time were controlled by use of a smoothing spline with 12 df (about 4 df per year of study data) based on the average of the minimization of the AIC (df = 15) and the minimization of the PACF (df = 9). We calculated beta as a function of time df from 1 to 25. Using more degrees of freedom than 12 does not appreciably affect the beta for Hg but slightly increases the beta for XRF-S and slightly decreases the betas for As and Se. Using less than 9 degrees of freedom would decrease all species except Hg. Results are shown in the supplemental data.

Conclusions

The multipollutant approach to air pollution epidemiology, whether by use of MPMs or by use of source contributions, allows us to investigate the combined health effects of several pollutants. Since the use of an MPM reduces the amount of confounding among the species in the model, the combined beta gives a better indication of the combined effect than could be obtained if only SPM betas were available. Source contributions from SAMs also allow us to investigate the combined effect of the several pollutants included in the source signature. Further studies should examine other multipollutant groups, such as vehicle emissions and multipollutant groups based on other source categories.

The strong relationship of sulfate (measured as XRF S) with cardiovascular mortality (CVM) observed in this study indicates that the use of ZIP-Code-level mortality data provides higher and more statistically significant relative risks than obtained with county-level mortality data. Further work is needed to determine which population property is most important, a lower exposure error or a more homogeneous and lower SES.

The strong relationship of particulate Hg (collected on a filter) with CVM suggests that Hg may affect cardiovascular as well as neurological health. This relationship of particulate Hg with CVM is suspicious and needs to be confirmed by other studies. Therefore, further studies of the relationship of Hg with mortality are needed.

Comparison of the combined beta (β^c) from a CPM containing concentrations of two pollutants with the sum of alphas (α^s) from two SPMs of the same pollutants allows us to identify confounding and estimate how much the excess risk would be increased if we used α^s instead of β^c to calculate the relative risk. The difference between α_x from a SPM and the corresponding β_x allows us to identify confounding and estimate how much one pollutant is confounded by the other pollutant in the CPM. The difference between α¹ and β^c provides information on how much pollutant 2 contributes to β^c from a CPM with pollutants 1 and 2. If the contribution to β^c is small or negative, it is likely that pollutant 2 is noncausal. As and Se were found to be noncausal. Their SPM relationship with CVM may be attributed to their correlation with S. This study provides no support for an independent effect of As or Se on CVM.

Since As and Se are noncausal and Hg was not correlated with sulfate, sulfate is the only species left to represent emissions from both coal power plants and smelters. Using FA to yield a regional sulfate source reduces the error due to sulfate in soil and yields a higher t-statistic and percent excess risk. Further splitting regional sulfate into a coal power plant source and smelter sources further increases the t-statistic and percent excess risk when these sources are combined in a CPM (PMF) or an MPM (PCA).

Acknowledgments

Jane Koenig originated the technique of reducing exposure error by limiting the population in an epidemiologic study to a small area whose air pollutant concentrations were well represented by the monitoring site. Dr. Koenig also obtained from the state of Arizona the ZIP Code mortality data used in this study. Therese Mar conducted the statistical analyses for the earlier Phoenix studies and performed some preliminary analyses for this study. Larry Purdue, Robert Stevens, and Roy Zweidinger organized the EPA monitoring study that provided the pollutant and weather database used in this study. Thanks to Gary Norris for providing source contributions from factor analysis and positive matrix factorization and to Quingyu Meng for providing source contributions from principal component analysis. Thanks to Tom Lumly for writing a variance–covariance matrix code to calculate correct standard errors in the combined model. Thanks to David Svendsgaard, Lucas Neas, and Jason Sacks for helpful reviews. An early version of this paper was presented at the 11th International Congress on Combustion By-Products and their Health Effects, Research Triangle Park, NC, May 31–June 3, 2009.

Disclaimer

Some early work on this paper was done while the author was an employee or volunteer (Emeritus Program) at the U.S. Environmental Protection Agency. However, much of the work was done after the author retired. The views expressed in this paper are those of the author and do not necessarily reflect the views or policies of the U.S. Environmental Protection Agency.

Supplemental Materials

Supplemental data for this article can be accessed at http://dx.doi.org/10.1080/10962247.2015.1033067.

Additional information

Notes on contributors

William E. Wilson

William E. Wilson retired in 2011 from the National Center for Environmental Assessment of the U.S. Environmental Protection Agency.