Environ. Sci. Technol. 2001, 35, 2554-2561
Estimating the Long-term Trend in the Extreme Values of Tropospheric Ozone Using a Multivariate Approach KIM A. MENEZES AND THOMAS S. SHIVELY* Department of Management Science and Information Systems, University of Texas at Austin, Austin, Texas 78712
This paper deals with the identification and estimation of the long-term trend in the extreme values of tropospheric ozone after allowing for the confounding effects of meteorological conditions. The U. S. Environmental Protection Agency’s (EPA’s) National Ambient Air Quality Standard for ozone is stated in terms of exceedances of a specified threshold level. Therefore, the EPA is concerned with the long-term trend in the probability of an exceedance. A multivariate nonparametric probit regression model estimated within the framework of a hierarchical Bayes model is used to model the probability of an exceedance after allowing for the effects of changing meteorological conditions. There are three advantages to using this model. First, the trends estimated at each site in a region can be separated into a city-wide component that is common to all sites and a site-specific component that is unique to the individual site. Second, the hierarchical Bayes framework allows for combining information across monitoring sites to increase the information available regarding the trend at each individual site. Third, the nonparametric model does not require the a priori specification of the functional forms relating the probability of an exceedance to the meteorological variables. Ozone data from four Houston, Texas monitoring sites for the period 19811997 are analyzed. We find that there is a downward trend in the probability of an exceedance in the 1980s followed by a relatively flat trend in the 1990s.
1. Introduction High levels of tropospheric ozone can cause serious damage to human health. Shortness of breath, nausea, eye and throat irritation, and lung damage are only a few of the long list of health problems caused by high ozone concentrations. The U.S. Environmental Protection Agency’s (EPA’s) National Ambient Air Quality Standard (NAAQS) for ozone is stated in terms of exceedances of a specified threshold level. When this standard is exceeded at a single monitoring site, the EPA classifies the surrounding county or metropolitan area as being in “nonattainment” for ozone. In Texas, the Houston area has been classified by the 1990 Clean Air Act to be a severe nonattainment area that needs to come into compliance by the year 2007. To determine if the Houston area is moving toward attainment, the EPA and the Texas Natural Resource Conservation Commission (TNRCC) are * Corresponding author phone: (512) 471-1753; fax: (512) 4710587; e-mail:
[email protected]. 2554
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 35, NO. 12, 2001
interested in the long-term trend in the probability of an exceedance of the threshold level specified by the NAAQS for ozone. In this paper we analyze daily tropospheric ozone data collected at four monitoring sites in Houston. A multivariate nonparametric probit regression model embedded in a hierarchical Bayes framework is used to model the longterm trend in the probability of an exceedance. This model accounts for the correlation in ozone levels across the four sites and for the effects that meteorological conditions have on the formation of ozone. In addition, the model allows the trend at each site to be separated into city-wide and sitespecific trends. This gives the EPA and TNRCC the ability to get a measure of how the city of Houston is doing overall as well as to determine which sites need specific attention (due, for example, to changing traffic patterns and industrial concentrations). The hierarchical Bayes model also allows for information regarding exceedance probabilities to be combined across monitoring sites in an appropriate way to increase the information available regarding the trend at each individual site. An excellent reference for hierarchical Bayes models is Gelman et al (1). Combining information across sites is important for two reasons. First, in some years certain sites do not have many data available because of the monitoring instruments being inoperable for long periods of time. For example, at one site in 1993 there were only twelve observations available and no exceedances occurred for these twelve observations. It is unrealistic to believe that if all the observations were available there would have been no exceedances recorded in 1993, i.e., it is unrealistic to believe that the probability of an exceedance at the site in 1993 was zero. There were many observations available at the other three sites in 1993 so it is reasonable to use the information from these sites to modify the estimated exceedance probabilities at the site with so few observations. The hierarchical Bayes model combines information across sites in a statistically meaningful way to obtain better estimates of the exceedance probabilities at each individual site. Second, there is evidence (which is difficult to confirm several years later) that some of the ozone monitoring instruments in the Houston network may have underreported ozone levels for scattered months in the period 1981-1997 because of calibration problems. By combining information across sites the Bayesian model will help mitigate potential data problems. The reason is that if one site reports unusually low ozone levels in a given year relative to those of other sites, the model’s ability to combine information across sites will reduce the resulting impact of the underreported ozone on the site’s trend function. Assessing the trend in ozone levels is complicated by the fact that the chemical reactions that form ozone in the ambient air are strongly affected by meteorological conditions. For example, ozone levels tend to be high in hot, sunny conditions because sunlight provides the energy for the chemical reaction that forms ozone. Conversely, windy conditions are typically associated with low levels of ozone because high winds tend to disperse the chemical precursors that react to form ozone. In addition, relative humidity and the direction of the wind, as well as interactions between the meteorological factors, can have an effect on ozone. It is therefore important to properly account for the effects of meteorological conditions when determining the long-term trend in ozone so that one or two years of atypical meteorological conditions do not distort the results. However, there is no chemical theory that specifies the functional form of the relationships between ozone and the meteorological 10.1021/es001838p CCC: $20.00
2001 American Chemical Society Published on Web 04/25/2001
variables. The nonparametric probit regression model used in this paper does not require the a priori specification of the functional forms relating the probability of an ozone exceedance to the meteorological variables. The data are allowed to specify the functional forms, so, given the large amount of data that are available, there is little risk of specifying an incorrect model and thereby invalidating the empirical results. In recent years various statistical techniques that incorporate the influence of meteorological variables have been used to assess the trend in ozone exceedances. One common approach is to use parametric regression models to link some characteristic of ozone, such as the mean level of ozone, to meteorological variables. Bloomfield et al (2) constructed a nonlinear regression model for hourly ozone data in the Chicago area in which meteorological variables, seasonality, and a trend term are included. However, regression models based on the assumption of normally distributed errors often underpredict the extreme ozone values as pointed out in the report of the National Research Council (3). To avoid this problem, Cox and Chu (4) proposed a model for the daily maxima of hourly ozone concentrations based on the Weibull distribution in which the scale parameter is allowed to vary as a function of meteorological conditions. Huang and Smith (5) used regression tree models to analyze high levels of ozone concentrations in the Chicago area, whereas Smith and Shively (6) used nonhomogeneous Poisson processes and generalized Pareto distributions to model the trend in the frequency and size of the exceedances of specified threshold levels after accounting for the effects of meteorological conditions. Niu (7) and Shively and Sager (8) extended the work of Cox and Chu (4) and Bloomfield et al. (2) by using nonparametric regression models to model ozone. Section 2 contains a description of the ozone and meteorological data used in this paper. Section 3 presents the multivariate nonparametric probit regression model used to analyze the ozone data, and Section 4 contains a discussion of the empirical results and their interpretation. The Supporting Information provides the statistical methodology used to estimate the functions nonparametrically.
2. Ozone Data The data used for this study consisted of the daily maximum ozone readings and the daily values of several meteorological variables collected at four monitoring sites in the Houston, TX area: Deer Park (DP), Clinton Drive (CL), Aldine Drive (AL), and Northwest Harris County (NH). The data are for the months May through October in the period 1981-1997. Only data from the “high ozone season” are analyzed because most of the exceedances of the EPA’s NAAQS for ozone occur during this period. The TNRCC believes the primary meteorological variables that affect the chemical reaction that forms ozone are sunlight, wind speed, relative humidity, and wind direction. Hence, these are the explanatory variables we include in our model. Sunlight is an important factor because it provides the energy for the chemical reaction that forms ozone. A direct measure of sunlight is not available at the monitoring sites, but the daily temperature range (daily maximum minus daily minimum) provides a good proxy for sunlight because the more sunlight there is the more the ambient air warms. Wind speed is also an important factor in the formation of ozone because high winds tend to disperse the chemical precursors to ozone so the chemical reaction that forms ozone cannot take place. In addition, there is a strong a priori belief that an interaction exists between sunlight (temperature range) and wind speed because if the wind speed is high enough the chemical precursors will disperse and less reactions will take place no matter how much sunlight there is.
Relative humidity is a measure of the amount of water vapor that is present in the ambient air. High relative humidity readings are typically associated with low levels of ozone because the water vapor tends to wash ozone from the ambient air. Conversely, low relative humidity readings are more likely to be associated with high levels of ozone. As with temperature range and wind speed there is an a priori belief that an interaction exists between humidity and wind speed. There is also reason to believe an interaction effect exists between relative humidity and temperature range. However, including a multiplicative interaction term in the model for each pairwise combination of temperature range, wind speed, and relative humidity induces near perfect multicollinearity so an interaction term for temperature range and humidity is not included. Wind direction is also a potential factor in the level of ozone present in the ambient air because winds from certain directions can blow in the chemical precursors to ozone. Following are the ozone and meteorological variables used in the analysis along with the hypothesized effect of the meteorological variables on the ozone level: (1) ozone The daily ozone value is the maximum of the 13 hourly ozone readings from 6 am through 7 pm; (2) temperature range Temperature range is the difference between the maximum and the minimum hourly temperature readings from 6 am through 7 pm. Large temperature ranges tend to be associated with high ozone levels. (3) wind speed - Wind speed is the average of the hourly wind speed readings from 6 am through 7 pm (in miles per hour). High wind speeds tend to be associated with low ozone levels. (4) relative humidity Relative humidity is the average of the hourly relative humidity readings from 6 am through 7 pm. High relative humidity readings tend to be associated with low ozone levels. The daily wind direction data give the percentage of time from 6 am through 7 pm that the hourly wind direction fell into one of four quadrants. Each quadrant is 90 degrees, with the first quadrant going from 90 to 180 degrees (0 degrees is north), and each subsequent quadrant obtained by adding 90 degrees. For example, if the wind blows from the Southeast (Quadrant 1) for 8 h and the Southwest (Quadrant 2) for 5 h, then the data for the four quadrants would be 0.62, 0.38, 0.00, and 0.00. Quadrant 1 provides the predominant wind direction in Houston relative to major nearby pollution sources, i.e., the TNRCC believes that wind blowing from the first quadrant will have the strongest impact on ozone levels. Only the first three wind quadrants are included in the analysis because the sum of all four wind quadrants add up to 1.00 and hence the fourth quadrant adds no additional explanatory power to the model. Other analyses of the long-term trends in ozone use similar sets of meteorological variables. For example, Bloomfield et al. (2) used temperature, humidity, wind speed, wind direction, visibility, cloud cover, and barometric pressure in an analysis of Chicago ozone data. Smith et al. (6) used temperature, wind speed, and wind direction variables in their analysis of Houston ozone data, whereas Niu (7) uses temperature, wind speed, humidity, and cloud cover variables in an analysis of Chicago ozone data. The missing value convention for the data is the following: if seven or more hourly observations are missing on a given day for ozone or any of the meteorological variables, then the data for that day are considered to be missing. Between 20% and 30% of the observations are missing at any given site. After allowing for missing observations there are 9553 observations available for analysis at the four sites combined. Figure 1a-d plot the percentage of days in each year at each site for which there was an exceedance of the NAAQS for ozone. Figure 1 panels a-d provide an initial feel for the trend in the probability of an exceedance at each site, but VOL. 35, NO. 12, 2001 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2555
FIGURE 1. Panels a-d plot the percentage of days in each year at each site for which there was an exceedance of the NAAQS for ozone. Panels e-h plot the percentage of exceedances at the four sites combined on days when wind speed, temperature range, and relative humidity take on specified values. they can also be somewhat misleading because no attempt is made to control for the meteorological variables. Autocorrelation coefficients for each of the time series of percentage exceedances in Figure 1a-d were computed. However, none of the autocorrelation coefficients are more than 1.3 standard deviations from zero so there is little statistical evidence that the annual exceedance percentages are related through time. Figure 1e graphs the percentage of exceedances (at the four sites combined) on days when wind 2556
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 35, NO. 12, 2001
speed takes on the values 1, 2, ..., 12. Figure 1f and g provide the same graphs for temperature range and relative humidity, although the temperature range and relative humidity values are grouped into bins so that there are multiple observations in each bin. For temperature range, the bin width is two degrees Celsius and for relative humidity the bin width is 0.05. As we would expect, the percentage of exceedances is increasing as the temperature range increases, and decreasing as wind speed increases. For relative humidity,
TABLE 1. Pearson Correlation Coefficients for the Ozone Data Deer Park Clinton Aldine NW Harris
Deer Park
Clinton
Aldine
NW Harris
1.000 0.830 0.659 0.555
1.000 0.742 0.694
1.000 0.778
1.000
the percentages peak at about 0.50. Figure 1h provides a three-dimensional projection graph of the percentage of exceedances for different combinations of wind speed and temperature range. Table 1 contains the pairwise correlations across sites for the ozone data. The correlations for the ozone data vary from 0.56 to 0.83 and indicate the potential importance of taking into account in the formal statistical model the correlation in the data across sites.
3. Statistical Model for Ozone Exceedances This section develops the nonparametric probit regression model used to model the long-term trend in the probability of an exceedance of a specified threshold level at the four sites in Houston. The model is constructed so that the trend at each site can be divided into city-wide and site-specific components. The functions relating the probability of an exceedance to the various meteorological variables are assumed to be the same at each of the four sites, but no assumptions are made regarding their functional form. Ozone levels, and therefore exceedances, are assumed to be correlated across sites. Let Ozj(t) represent the maximum hourly ozone value at site j on day t, L represent a specified threshold level and
Ij(t) )
{
1 if Ozj(t) > L 0 if Ozj(t) e L
17
∑ k)1
8
Rkjskj(t) +
∑ f (w (t))] l
lj
(1)
l)1
where Φ is the standard normal cumulative distribution function, skj(t) is defined for each site j as
skj(t) )
{
j k + ηkj Rkj ) R
1 if day t is in year k 0 otherwise
(2)
wlj(t) represents the value of the l-th meteorological variable at site j on day t, and fl is the unknown function relating the l-th meteorological variable to the probability of an exceedance. f1 - f8 represent the functions for wind speed; temperature range; relative humidity; the interaction terms for wind speed and temperature range, and wind speed and relative humidity; and the three wind direction variables, respectively. Sites 1, 2, 3, and 4 are the Deer Park, Clinton, Aldine, and Northwest Houston County sites, respectively. Note that the Rkj, k ) 1, ..., 17, represents a measure of the long-term trend in the probability of an exceedance at site j for the 17 year period 1981-1997 after allowing for the effects of meteorological conditions. For example, if the Rkj are decreasing, then the probability of an exceedance is decreasing through time at site j for a given set of meteorological conditions. In Section 4, we provide estimates of the Rkj and fl, and then compute and graph the annual probability of an exceedance at each site for a set of meteorological conditions. 3.1. Modeling the Trend Terms Rkj. There is strong reason to believe a priori that the trend terms across sites in a given
(3)
where the ηkj are identically distributed N(0,σ2η) random variables and independent across years and sites (i.e., independent for all k and j). R j k represents the effects of citywide factors on the probability of an exceedance at each site in year k, while ηkj represents the effects of site-specific factors on the probability of an exceedance at site j in year k. σ2η is a measure of the variability in the site-specific influences across the four sites with σ2η ) 0 indicating that there are no site-specific differences. Combining eqs 1 and 3 gives 17
pr(Ij(t) ) 1) ) Φ[
∑
17
R j kskj(t) +
k)1
∑
8
ηkjskj(t) +
k)1
∑ f (w (t))] l
lj
l)1
(4)
The city-wide average probability of an exceedance in year k for a given set of meteorological conditions is computed using
∫
∞
-∞
The probability that ozone exceeds the threshold level L at site j on day t is modeled as
pr(Ij(t) ) 1) ) Φ[
year are related. The reason is that there will be, almost certainly, factors other than meteorological conditions (including the effect of pollution control programs) that affect ozone city-wide and therefore have an impact on ozone levels at all sites in Houston. Conversely, there will be factors that affect the formation of ozone only at specific sites (i.e., traffic patterns and industrial activity) so the Rkj will also have a component that is individual to the site. To capture the citywide and site-specific influences in the trend term, we model Rkj as
8
Φ[Rk +
∑ f (w (t))] f(R )dR l
lj
k
(5)
k
l)1
where f(Rk) is the density function for a N(R j k, σ2η) random variable. The expression in eq 5 represents the average of the 8 exceedance probabilities Φ[Rk + ∑l)1 fl(wlj(t))] with respect to the distribution of the trend term Rk. The change through time in this probability reflects the effect of city-wide factors, other than meteorology, on the long-term trend in the probability of an exceedance at each site. The difference between the city-wide average probability of an exceedance and the probability of an exceedance at a specific site (computed in eq 1 ) represents the effect of site-specific factors on the probability of an exceedance at the site. From a regulatory point of view, a downward trend in the city-wide average probability of an exceedance would provide evidence that pollution control programs implemented city-wide are having a positive (downward) effect on ozone levels. Knowing the effect of site-specific factors on the probability of an exceedance allows the TNRCC to determine which sites need specific attention, and possibly to employ additional resources in the areas surrounding these sites to alleviate the problem. An alternative method for modeling the city-wide average probability of an exceedance is to fit a univariate probit regression model to the exceedances at each site individually and then compute the city-wide average probability of an exceedance for a given set of meteorological conditions by averaging the four individual probabilities. There are three advantages to using the multivariate model outlined in this section to model ozone exceedances. First, the model in this section allows us to account for the correlation in the ozone exceedances across sites and to incorporate this information into the model. Second, using the multivariate model, and more specifically the expression in eq 3 for Rkj, allows for combining information across monitoring sites to increase the information available regarding the trend at each individual site. Also, if the monitoring instrument at one of VOL. 35, NO. 12, 2001 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2557
the sites is incorrectly calibrated, then combining information across sites will help mitigate its effect on the analysis. Third, using the expression in eq 5 to compute the city-wide average probability reduces the impact that an “outlying” estimated value of Rkj will have on the city-wide average probability. This reduces the effect on the city-wide average probability in a given year of an unusually large (or small) number of exceedances that may occur at a specific site because of random variability. 3.2. Modeling the Functions fl. We model the unknown functions fl using the expression
fl(wl) ) βl0 + βl1(wl - wl0) + τl
∫
wl
wl0
(wl - s)Wl(ds) (6)
where wl0 is the minimum value of wl, βl0 is the function value at wl0, i.e., βl0 ) fl(wl0), βl1 is the slope of the function at wl0, i.e., βl1 ) f′(wl0), τl is a smoothing parameter, Wl(s) is a Wiener process with unit dispersion parameter, and Wl(ds) is the change in the process over an interval of length ds. Time is suppressed in the above notation. The functions fl can be estimated only up to a constant because the βl0’s cannot be identified, just as in a multiple linear regression setting. The representation of fl(wl) in eq 6 corresponds to Wahba’s (10) specification of a cubic spline function, whose slope changes stochastically according to a Wiener process. This representation provides a very flexible class of functions with which to model the relationship between ozone exceedances and the meteorological variables. It is also possible in principle to determine whether fl(wl) is a flat function, i.e., whether the meteorological variable wl is related to the probability of an exceedance, using Bayesian variable selection techniques. A Bayesian variable selection technique is developed for nonparametric univariate probit regression models in Shively et al (9). The development of a similar technique for the multivariate probit model in this paper is an issue for further research. 3.3. Modeling the Correlation in Ozone Exceedances Across Sites. To model the correlation in ozone exceedances across sites we use the latent variable approach to probit regression of Chib and Greenberg (11) and Albert and Chib (12). An excellent reference for probit regression models and how latent (unobservable) variables can be used to motivate this type of model is given in Pindyck and Rubinfeld (13, pp 280-287). Let 17
yj(t) )
∑ Rj s
k kj(t)
k)1
17
+
∑η k)1
8
kjskj(t)
+
∑ f (w (t)) + (t) l
lj
j
(7)
l)1
For a given site j, the j(t), t ) 1, ..., Tj, are assumed to be independent N(0,1) random variables. Let Ij(t) ) 0 if yj(t) < 0 and Ij(t) ) 1 if yj(t) g 0. It is readily checked that Ij(t) satisfies eq 1. To model the correlation in ozone levels across sites on day t, and therefore the correlation between ozone exceedances on day t, we assume the j(t), j ) 1, 2, 3, 4, are correlated. Letting (t) ) (1(t), 2(t), 3(t), 4(t))′, we have (t) are independent identically distributed N(0, D) random vectors where D has ones on its diagonal and whose nondiagonal elements need to be estimated. 3.4. Prior Distributions. Diffuse priors are placed on the R j k and βl1. Flat priors on the interval [0,1010] are used for σ2η and the smoothing parameters τ2l . A multivariate N6(0, 0.5I6) prior truncated to the appropriate region C, where C is a convex solid body in the hypercube [-1,1]6, is used for the nondiagonal elements of D (see Chib and Greenberg (11) for the rationale underlying this prior). 3.5. Estimation Procedure. The unknown parameters and functions in the multivariate probit model are estimated using a latent variable approach from Wood and Kohn (14) and Albert and Chib (12). The estimation procedure is imple2558
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 35, NO. 12, 2001
mented using a Markov chain Monte Carlo (MCMC) sampling scheme similar to one proposed by Wood and Kohn (14) for a univariate probit model. We modify Wood and Kohn’s sampling scheme to exploit the structure of the data set in order to provide the computational efficiencies necessary to implement the estimation procedure for a sample containing 9553 observations. Details regarding the implementation of the estimation procedure are given in the Supporting Information.
4. Empirical Results Section 4.1 reports the estimation results for the trend functions that are obtained when the multivariate probit model is applied to the combined data from the Deer Park, Clinton, Aldine, and Northwest Harris County sites. Section 4.2 discusses the differences between the sites in the trend functions. Section 4.3 reports the relationship between the exceedance probabilities and the meteorological variables, and Section 4.4 discusses the estimates of the correlation matrix D. Section 4.5 contains a comparison of the trend function estimates obtained using the multivariate probit model to those obtained using four individual probit regression models. 4.1. Trends in the Probability of an Exceedance. The solid lines in Figure 2a-d provide the annual probability of an exceedance at the four sites of the threshold level 124 ppb (which is the threshold level specified by the EPA’s NAAQS in effect during the period 1981-1997) when the meteorological conditions are conducive to the formation of ozone. Conducive meteorological conditions are defined as wind speed ) 2 mph, temperature range ) 12 °C, relative humidity ) 60%, Quadrant 1 ) 0.46, Quadrant 2 ) 0.23, and Quadrant 3 ) 0.31. The exceedance probabilities are computed using eq 1 with the estimates of the Rkj and the meteorological functions f1 - f8. The city-wide average probability of an exceedance for each year (short-dashed line) is also given in Figure 2a-d for the same meteorological conditions. The city-wide average probabilities are obtained using eq 5 . The long-dashed lines in Figure 2a-d are 90% confidence bands for the city-wide average probability. Figure 2 reveals a downward trend in the early and mid 1980s in the probability of an exceedance at each site and in the city-wide average, and then a flattening out of the trend in the late 1980s and throughout the 1990s. The 90% confidence bands indicate that the changes in the annual exceedance probabilities probably cannot be attributed to random variability in the data. For conditions that are conducive to the formation of ozone the city-wide average decreases from 0.74 in 1981 to 0.56 in 1987 (which is a decline of 24%), and then fluctuates around 0.56 between 1988 and 1997. The interpretation of the decline between 1981 and 1987 is that on days when meteorological conditions are conducive to the formation of ozone there is a 24% decline in the city-wide average probability of an exceedance due to nonmeteorological factors. Thus, the results suggest that citywide pollution control programs implemented in Houston in the early and mid 1980s had a positive impact on ozone. However, the flattening out of the probabilities in the late 1980s and throughout the 1990s suggest that the city may no longer be moving closer to attainment. Similar results are obtained when the meteorological variables are set to their median values (wind speed ) 3 mph, temperature range ) 9 °C, relative humidity ) 67%, Quadrant 1 ) 0.48, Quadrant 2 ) 0.54, and Quadrant 3 ) 0.08). For these meteorological conditions the city-wide probability of an exceedance decreases from 0.15 in 1981 to 0.06 in 1987 (which is a decline of 60%), and then fluctuates around 0.06 between 1988 and 1997. These results can be interpreted to reflect the exceedance probabilities under “typical” meteorological conditions.
FIGURE 2. Panels a-d display the annual probability of an exceedance at each site (solid line) and the city-wide average probability of an exceedance (short-dashed line) estimated using the multivariate probit regression model for a set of meteorological conditions that are conducive to the formation of ozone. The long-dashed lines are 90% confidence bands for the site-specific probabilities. 4.2. Differences in the Trends Between Sites. As discussed in Section 3, the model allows for differences between sites to account for site-specific effects such as changing traffic patterns and industrial concentrations. For example, Figure 2a shows that the probability of an exceedance at Deer Park is as high or higher than the city-wide probability of an exceedance in every year except 1992 and 1993 (there are only twelve observations at Deer Park in 1993). This indicates that Deer Park has an ozone problem that is consistently worse than the other three sites and that it may be useful to focus additional resources for reducing ozone in the Deer Park area. To determine whether the differences we observe graphically between the site-specific probabilities and city-wide average probabilities are real (i.e., cannot be attributed to random variability) we need to determine whether σ2η ) 0. Figure 3 contains a histogram of the σ2η values generated by the MCMC sampling scheme and therefore provides an approximation to the posterior distribution of σ2η. The mean of the generated σ2η values is 0.14 with a standard deviation of 0.04. The location and dispersion of the histogram implies there is strong evidence that σ2η > 0 and therefore that the probability of an exceedance in a specific year differs across sites. 4.3. Estimates of the Meteorological Functions. Figure 4a plots the probabilities of exceeding the threshold level of 124 ppb at the Aldine site in 1987 for different combinations of wind speed (WS) and temperature range (TR) when relative humidity (RH) is set to its 10th percentile (55% humidity). These probabilities are computed assuming the wind blows from Quadrant 1 for 6 h (46% of the time), from Quadrant 2 for 4 h (31% of the time), and from Quadrant 3 for 3 h (23% of the time). The probabilities are obtained from eq 1 by combining the trend function estimate for Aldine in 1987
FIGURE 3. Histogram of the σ2η values generated by the MCMC sampling scheme. (i.e., the estimate of R7,2) and the estimates of the wind direction functions (i.e., f6(0.46) + f7(0.31) + f8(0.23)) with the estimate of f(WS,TR,RH) ) f1(WS) + f2(TR) + f3(RH) + f4(WS × TR) + f5(WS × RH) when RH ) 0.55. Figure 4b and c provide the corresponding probability plots when relative humidity is set to its 50th percentile (67% humidity) and 90th percentile (81% humidity). The estimated probability functions show that the effect of increasing wind speed is to reduce the probability of an exceedance for all temperature range and relative humidity values. They also show a strong interaction effect between wind speed and temperature range because for low wind speed, the probability of an exceedance increases as temperature range increases (except for very large values of temperature range) while for high wind speed increasing temperature range has no effect. Finally, increasing relative VOL. 35, NO. 12, 2001 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2559
TABLE 2. Estimates of the Off-Diagonal Elements of the D Matrix (Standard Errors are in Parentheses) Deer Park Deer Park Clinton Aldine NW Harris
FIGURE 4. Panels (a)-(c) plot the probability of an exceedance of the threshold level 124 ppb in 1987 at the Aldine Drive site for different combinations of wind speed and temperature range when relative humidity is set to its 10th (a), 50th (b), and 90th (c) percentile values (55%, 67%, and 81% humidity), respectively. It is assumed that the wind blows from Quadrant 1 for 6 h, from Quadrant 2 for 4 h, and from Quadrant 3 for 3 h. humidity lowers the probability of an exceedance. The meteorological reasons for the properties of these probability functions have been discussed in Section 2. The dip in the probability function for high values of temperature range and low values of wind speed is probably due to the fact that there are very few data points in this region of the space, rather than being attributable to a meteorological phenomenon. 4.4. Correlation Between Sites. The correlation matrix D defined in Section 3.3 provides a measure of the strength of the relationship between ozone at different sites after the effects of the long-term trend and meteorological conditions are accounted for. The estimates and corresponding standard errors of the off-diagonal elements of D are given in Table 2. The four sites lie very close to a straight line running from southeast to northwest. Ordered from the southeast to northwest the sites are Deer Park, Clinton Drive, Aldine Drive, 2560
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 35, NO. 12, 2001
1.000 0.369 (0.051) 0.100 (0.059) -0.017 (0.065)
Clinton
Aldine
NW Harris
1.000 0.183 (0.053) 0.099 (0.060)
1.000 0.383 (0.046)
1.000
and Northwest Harris County. The three highest correlations are between the adjacent sites of Deer Park and Clinton, Clinton and Aldine, and Aldine and Northwest Harris County, with the smallest correlations between the sites furthest apart: Deer Park and Northwest Harris County, and Clinton and Northwest Harris County. The ratio of each estimated correlation coefficient to its standard error indicates there is strong evidence that the correlations between Deer Park and Clinton, Clinton and Aldine, and Aldine and Northwest Harris County are positive, and moderately strong evidence that the correlations between Deer Park and Aldine, and Clinton and Northwest Harris County are positive. The positive correlations indicate the importance of using the multivariate probit model rather than four univariate probit models that do not account for the relationships between the sites. A comparison between the exceedance probabilities obtained from the multivariate and univariate models is given in Section 4.5. An alternative parametrization of the correlation matrix would account for the spatial orientation of the sites, i.e., allow the correlation between two sites to be a function of the distance between the sites. However, given the large amount of data available for each site we do not believe imposing such a restriction on the correlation matrix is necessary to obtain reliable estimates. 4.5. Comparison of Results with Alternative Models. An alternative method that can be used to model the probability of an exceedance is to use a Gaussian regression model with the natural logarithm of ozone as the dependent variable (see ref 2). However, the Anderson-Darling test for normality performed on the residuals from a Gaussian model applied to the Houston data indicates the assumption of normally distributed errors is violated (p-value < 0.0001). This implies that the natural logarithm is not the correct transformation of ozone to use with the Gaussian model, and any exceedance probabilities computed using this model may be in error. Further, the identity of the appropriate transformation is not clear. It can be shown that the exceedance probabilities computed using a nonparametric probit model are robust to the distribution of ozone. Therefore, we are better off using the nonparametric probit model to model exceedance probabilities because this model does not rely on any distributional assumptions regarding ozone. The trend in the probability of an exceedance for each site was also obtained by fitting a univariate probit regression model to the exceedances at each site individually. The individual univariate models were estimated using the same methodology used for the multivariate model except the Rkj trend terms in eq 1 are assumed to be unrelated across sites and the ozone values across sites are assumed to be independent. The city-wide average probability of an exceedance for a given set of meteorological conditions is computed using the arithmetic average of the four individual probabilities. Figure 5a-d displays the trends obtained from the four univariate models for meteorological conditions that are conducive to ozone formation. The effect of combining
FIGURE 5. Annual probability of an exceedance at each site (solid line) and the city-wide average probability of an exceedance (shortdashed line) estimated using four univariate probit regression models to model ozone exceedances at each site individually: panel a, Deer Park; b, Clinton; c, Aldine; d, NW Harris County. information across sites on the trend at a site in the multivariate model is visible at all four sites by comparing the results of the univariate and multivariate models (i.e., by comparing Figure 2a-d with Figure 5a-d). For example, there is considerably more fluctuation from year to year in the probability of an exceedance at the Deer Park site for the univariate model than for the multivariate model. The effect of combining information because of the lack of data is particularly strong at the Deer Park site in 1993 because it recorded only twelve observations that year. The estimated probability of an exceedance in 1993 for the univariate model (Figure 5a) is zero because there are no exceedances at Deer Park that year whereas the estimated probability in the multivariate probit model (Figure 2a) is very close to the city-wide average probability. A similar effect is seen for the Clinton site in 1991 when there are only 10 observations. The trend in the probability of an exceedance for each site was also estimated using a univariate probit regression model for each individual site that did not include meteorological variables. The trend function estimates are very similar to the plots in Figure 1a-d.
Acknowledgments The authors thank Dr. M. W. Hemphill and three anonymous reviewers for the many helpful comments they made on this manuscript. The views expressed herein are those of the authors and do not necessarily represent the position of the
TNRCC. We are grateful to the TNRCC for providing the data analyzed in this study.
Supporting Information Available Details regarding the implementation of the estimation procedure. This material is available free of charge via the Internet at http://pubs.acs.org.
Literature Cited (1) Gelman, A., Carlin J. B.; Stern H. S.; Rubin, D. B. Bayesian Data Analysis; Chapman and Hall: New York, 1995. (2) Bloomfield, P.; Royle. J. A.; Steinberg, L. J.; Yang, Q. Atmos. Environ. 1996, 30, 3067. (3) National Research Council. Rethinking the Ozone Problem in Urban and Regional Air Pollution; National Academy Press: Washington, DC, 1991. (4) Cox, W. M.; Chu, S. H. Atmos. Environ. 1992, 27B, 425. (5) Huang, L. S.; Smith, R. L. Environmetrics 1999, 10, 103. (6) Smith, R. L.; Shively, T. S. Atmos. Environ. 1995, 29, 3489. (7) Niu, X. J. Am. Stat. Assoc. 1996, 91, 1310. (8) Shively, T. S.; Sager, T. W. Environ. Sci. Technol. 1999, 33, 3873. (9) Shively, T. S.; Kohn, R.; Wood. S. J. Am. Stat. Assoc. 1999, 94, 777. (10) Wahba, G. J. R. Statist. Soc. B. 1978, 40, 364. (11) Chib, S.; Greenberg, E. Biometrika 1998, 86, 221. (12) Albert, J.; Chib S. J. Am. Stat. Assoc. 1993, 88, 669. (13) Pindyck, R. S.; Rubinfeld, D. L. Econometric Models & Economic Forecasts; McGraw-Hill: New York, 1981. (14) Wood, S.; Kohn, R. J. Am. Stat. Assoc. 1998, 93, 203.
Received for review November 2, 2000. Revised manuscript received March 7, 2001. Accepted March 14, 2001. ES001838P
VOL. 35, NO. 12, 2001 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2561
