A Risk Assessment for Ozone Regulation Based on Statistical Rollback

In environmental studies, it is important to assess how regulatory standards for air pollutants affect public health. High ozone levels contribute to harmful air pollutants. The EPA regulates ozone levels by setting ozone standards to protect public health. It is thus crucial to assess how various regulatory ozone standards affect non-accidental mortality related to respiratory deaths during the ozone season. The original rollback approach provides an adjusted ozone process under a new regulation scenario in a deterministic fashion. Herein, we consider a statistical rollback approach to allow for uncertainty in the rollback procedure by adopting the quantile matching method so that it provides flexible rollback sets. Hierarchical Bayesian models are used to predict the potential effects of different ozone standards on human health. We apply the method to epidemiologic data.


Introduction
Regulating high ozone levels is essential as exposure to ozone can affect the risk of respiratory diseases or related deaths. Tropospheric ozone, also known as ground-level ozone, is one of the main harmful air pollutants that can cause adverse health effects. Many areas in the United States have been observed to exceed the current ozone National Ambient Air Quality Standards (NAAQS) . Elevated ozone concentrations have also been a growing concern for rapidly developing nations where emissions of ozone precursors have been risen from expanding transportation networks [1]. Recently, U.S. EPA promulgated a new ozone NAAQS as 0.070 ppm (parts per million) and this change has prompted studies on the effects of the new ozone standard [2][3][4]. To investigate how different ozone NAAQS affect public health and welfare can be of great interest in epidemiologic studies, toxicological studies, or controlled human exposure studies. Based on the reviews of the air quality criteria for ozone (O 3 ) and related photochemical oxidants and the NAAQS for O 3 , a modification in the current ozone regulatory standards provides the required protection for public health and welfare (see also [5,6]). Bell et al. [7,8] disapproved of the EPA's scientific reviews reporting that designations based on air quality data from 2006 to 2008 would be effective in air quality data to take effect in 2010 for the 2007 8-h ozone standard. Therefore, it would be interesting to investigate whether the current regulations are sufficiently stringent to prevent respiratory-related mortality or not.
To assess how changes in the ozone regulations affect mortality, ground-level ozone must be adjusted by strengthening the air quality standards. Although rollback functions, namely, air quality adjustment procedures proposed by the EPA [9], can be useful for adjusting the ozone process, it cannot introduce sufficient variability in the rollback adjustment as the adjustment is deterministic and the EPA regulatory standards are based on the average of three consecutive years' AQI values. Thus, we consider a parametric rollback approach with quantile matching method to allow for uncertainty in the rollback procedure To conduct the risk assessment, we consider hierarchical Bayesian models that provide uncertainty quantification for relevant parameters to predict the potential effects of different regulatory ozone standards. We also describe some variations in the results under different modeling assumptions. We analyze databases from the National Morbidity, Mortality, and Air Pollution Study [10] that is designed to study the public health effects of air pollutants. The Health Effects Institutes began this study with researchers from Johns Hopkins and Harvard University in 1996. The NMMAPS database contains 108 U.S. urban areas from 1987 to 2000 where one can build a multisite time-series model of ozone and mortality simultaneously with meteorology information such as temperature or dew point and air pollution (O 3 , PM 10 , SO 2 , NO 2 , and CO) (see also [11]).
The focus of this paper is on statistical approaches to describe how the new ozone regulations affect mortality. We shall describe our results that mortality decreases as limits of acceptable ozone level get lower through the statistical rollback approach.

Statistical Rollback Approach
The rollback transformation [9] is one of the methods to adjust current ozone processes to follow new ozone regulatory standards. There are two main conditions to consider in air quality data. The first one is baseline conditions characterized by unadjusted air quality data monitored at fixed locations during recent years. The other one is attainment conditions generated by fitting an air quality adjustment procedure (AQAP) to the baseline.

Quantile Matching Approach
Suppose that current ozone level, denoted by x, is a random variable having a distribution function F θ . We assume that adjusted ozone level corresponding to x under a new scenario follows the same distribution with different parameters, F θ * , and there exists a mapping m such that a random variable m(x) has the distribution function F θ * . Thus, the adjusted ozone level can be obtained using m(x). That is, θ * is the adjusted model parameter based on the relationship between current levels and adjusted levels under the new scenario. Moreover, let F(x; θ) and G(z; γ) be the ozone distribution functions based on the baseline condition and attainment condition, respectively, (e.g., G ≺ F). That is, ). Note that F and G can be estimated by baseline and attainment years' ozone data, respectively. Attainment year based on its AQI values or the condition G ≺ F is required. G may require adjustment based on ρ(s) = x(b) x(s) , where x(b) is current design values and x(s) is design values under scenario s.
Let f θ be a density function of current ozone levels x. Then adjusted rollback values z corresponding to x are obtained as follows: 1.

2.
Compute q x i , the quantile of x i such that for i = 1, . . . , n.

5.
{z i } n i=1 are adjusted (rollback) values of x. We extend this approach through the Bayesian framework by putting prior distributions of parameters, π(θ). By generating samples from posterior distribution, π(θ|x) ∝ f θ π(θ), we can generate rollback ensembles in the same fashion.
In general, the distribution of x can be expressed as a mixture of the form where u is a threshold, h can be any parametric distribution such as Weibull, truncated Gamma or Gaussian distribution, and g is a Generalized Pareto distribution (GPD) with parameters θ 2 = (σ, ξ). The quantile matching approach is simple and straightforward but causes various problems in adjusting multiple ozone processes in large dimensions. Now we consider more feasible approaches in the next section.

Weibull Approach in Rollback
Let y t (c) denote hourly (or daily average, daily maximum, daily 8-h maximum) ozone level for each time t and city c and y * t (c) denote rollbacked process corresponding to baseline process y t (c). That is, y * t = γ(y t ) η , where we omit city index for simplicity. We fit y t to a parametric distribution such as Weibull distribution.
We obtain following information for each year i and each city c: where κ and δ are weibull parameters. Note that γ and η are functions of κ and δ.

Now we consider a model of the form
For example, we consider a simple linear model: For new standard regulation AQI * , we can obtain κ * (c) and δ * (c) for each city c using It is observed that Σ(c) introduces some uncertainty in κ * (c) and δ * (c). In general, hierarchical structures for (α 1 (c), β 1 (c), α 2 (c), β 2 (c)) can be considered. Usually, δ provides a better fit based on ACLV1. Thus, we may need to adopt the relationship between κ and δ: Otherwise, we may select an average value or baseline value of κ or δ.
Let PAQI s (c) be the 4th largest 8-h daily maximum concentration of z t , which is Then,

Log-Normal Approach in Rollback
Let y t be a baseline process and y * t be rollbacked process. Assume that Similarily, we fit y t to a log-normal distribution for each time t and city c. It is observed that y t can be any covariate in the risk model (e.g., daily average, daily maximum, and daily 8-h maximum). Then, can obtain following information for each year i and each city c:

Now we consider a simple linear model:
For new standard regulation AQI * , we can determine µ * (c) and σ * (c) for each city c using where µ and σ are estimated based on baseline year. In general, the following relationship is considered: . Furthermore, well-fitted models of µ(c) and δ(c) can provide a good estimation based on following relationships We can estimate one parameter (e.g., δ * or µ * ) first and then tune other parameter (e.g., κ * or σ * ) to make rollbacked AQI set to be new standard regulation AQI * based on the range of the parameter. In the 8-h daily maximum rollback, we can directly solve the following equation by plugging in estimated parameters of δ * or µ * : Furthermore, the following empirical relationships can be useful.
As previous approaches (except for the 8-h maximum covariate) do not guarantee rollbacked values satisfying exact new standard regulation AQI * , hourly scale rollback will be more appropriate for adjustment after rollback. Notice that EPA standard regulation is based on the average of three consecutive years' AQI values.

Statistical Modeling
We apply an overdispersed Poisson model in generalized linear models (see [8,12]) to the NMMAPS data. Denote Y c t by the number of daily non-accidental deaths in community c on day t. The Poisson process with intensity function µ c t can be expressed as where the parameter φ c describes overdispersion for community c. Note that all overdispersed Poisson models are assumed to be mutually independent over time. We can also model the intensity function µ c t with some essential covariates such as ozone levels at different lags, seasonality, long-term trends, weather, and co-pollutants for three age groups (<65, 65-74, and ≥75 year). Natural cubic splines are useful tools for getting smoothing functions of time to account for seasonality and long-term trends in which influenza epidemics, for example, can affect mortality. The interaction term between smoothing functions of time and age-specific indicators (<65, 65-74, and ≥75 years) is considered as it can adjust the possible seasonal mortality patterns by age group. We also control for some potential covariates related to weather by smoothing functions of dew point, average dew points of the previous three days, temperature, the average temperature of the previous three days as follows: + interaction terms for age groups and time, where µ c t is the expected number of deaths, x c t is the average daily O 3 concentrations of the current and the previous day in community c on the day t, DOW t is the days of the week (categorical) on day t. We define ns(time, 7/year), ns(T c t , 6), ns(T c t−1,t−3 , 6), ns(D c t , 3) and ns(D c t−1,t−3 , 3) as the natural cubic spline function of time with 7 degrees of freedom per year, temperature with 6 degrees of freedom, the average temperature of the previous three days with 6 degrees of freedom, dew point with 3 degrees of freedom and the average dew points of the previous three days with 3 degree of freedom, respectively. The last term in (4) indicates the interaction terms of age-specific indicators and natural cubic spline functions of time.
Relative mortality rates associated with exposure to ozone over the past few days can be estimated in a specific community by constrained or unconstrained distributed-lag models as daily ozone levels are readily available. Those models can be more flexible in the sense that they can be more suitable for exploring the time lag between exposures to ozone and deaths than single-lag models (see [7]).
For example, the constrained distributed-lag models (CDL) and unconstrained distributed-lag models (UDL) can be expressed as respectively. The national average relative mortality rate caused by ozone effects can be estimated through a Bayesian hierarchical model in which both variability within-community and across-community can be accounted. The national average relative rate will be estimated by integrating estimates for the relative rate from the distributed-lag models for each specific community. Through the two-stage model, first, variation across communities is considered over the short-term ozone effects. The national average relative rate is then estimated. The ozone effects O 3 on mortality is modeled for each community c as follows: where β c is the true relative rate in a specific community,β c is its estimate, andΣ c is the estimated covariance matrix corresponding toβ c . We put the multivariate normal distribution on β c , β c |µ, Ω ∼ MV N(µ, Ω), where µ is the true national average relative rate and Ω is the true covariance matrix of relative rate in a specific community, β c . Note that there are still several sensitivity issues related to the modeling: (1) co-pollutant such as PM 10 can be included as a potential confounder; (2) days with high temperatures can be excluded to control for the effects of heat waves; (3) the degrees of freedom (df) in the smooth functions of time needs to be specified to control for long-term trends and seasonality; and (4) various ozone exposure metrics such as daily average, 1-h maximum, and 8-h maximum can be considered (see [7]).

Inferences
Estimating the impact of new ozone regulatory standards on the total nonaccidental deaths is not achieved through the national average relative rate µ as the rollback adjustment is applied to each community separately. Instead, we directly calculate the expected total nonaccidental deaths of the original and adjusted (rollbacked) ozone process, respectively.
Denote g c by a rollback function for community c, which adjusts observed ozone concentrations to meet new ozone regulation standards. Our focus is then on the difference in expected nonaccidental deaths between before and after rollback transformation. Let , and x c t is the design matrix of unconstrained (or constrained) distributed-lag ozones for community c (see [7]). If the amount of ozone reduction r is the same for each community c, then the reduction ratio of the expected total death is 1 − exp −rβ c for each community c at time t. Thus, the reduction rate of the national expected total death is 1 − exp(−rμ) at time t. As the µ is obscure, the reduction rate of the national expected total death can be better estimated by weighting mean effects across the 98 cities, where weights are proportional to populations. For simplicity, the unconstrained distributed-lag model is only considered here although the constrained distributed-lag models can be also readily applied via reparameterization. Let z t be rollbacked ozone concentrations of x t for community c (i.e., z c t = g c (x c t )). After the rollback transformation, the expected total deaths for community c during the ozone season is as follows where T O 3 is the set of time indices for the ozone season. The log ratio of national expected death before and after the rollback transformation during ozone season can be expressed as where C is the set of indices of communities in NMMAPS data and w c,t is the weight of community c at time t, We can estimate the posterior distribution of the log ratio in (9) through the MCMC approach. As π(β|Y) can be approximated by π(β|β) and a prior distribution of (µ, Ω) well, we can generate β from the posterior distribution π(β|Y). One of the efficient algorithms in application is TLNISE (two level Normal independent sampling estimation) algorithm (see [13]). Without a Bayesian approach, restricted maximum likelihood (REML) method can provide similar estimates.
Finally, with samples from π(β|Y) the posterior distribution of the log ratio of national expected death in (9) can be estimated as follows where β (i) indicates each sample from π(β|Y) or π(β|β) for i = 1, . . . , N.

Results
We investigate how the rollback transformation predicts that ozone series change under new ozone regulation standards possibly affecting public health. The transformation is based on the current AQI and new regulation standards in the attainment years. Figure 2 shows Q-Q plots of ozone levles under new regulation standards against the current AQI using various rallback functions and statistical rollbacks. A common AQI is applied to all cities (common rollback) or a different AQI is applied to each city (city-specific rollback).
The regular meteorological model in the NMMAPS is based on nonlinear functions mainly consisting of temperature and dewpoint at lag 0 and the average of them at lags 13. However, meteorological confounding in a distributed-lag model for ozone may exist at lags 46 so that we account for temperature and dewpoint at lag 0 through 6 in the "distributed-lag" meteorological model through nonlinear splines with 4 and 3 df, respectively.
We use the reduction rate of the expected total death to assess the effect of the proposed new ozone regulatory. Posterior means and 95% credible intervals for total mortality per 1000 deaths are shown in Tables 1 and 2. They are based on unconstrained distributed-lag models with a common rollback and city-specific rollback approach, respectively. Most credible intervals have positive lower limits in the common rollback approach except for Weibull rollback at level 75 regularization, providing good evidence that there are reduced rates. As expected, the mortality rate decreases as lower regulation increases. Similar results are observed in CDL. The ranges of the estimated reduction rates are approximately 0.94-1.9 for 75 ppb, 1.47-2.9 for 70 ppb and 4.9-5.0 for 60 ppb for all models and rollback functions. Statistical rollbacks tend to show similar results as quadratic rollbacks for overall regularization levels. The ranges of the estimated reduction rates are roughly 1.0-2.0 for 75 ppb, 1.4-2.8 for 70 ppb and 2.4-4.8 for 60 ppb for all rollback functions in the city-specific rollback approach. Posterior variances of the relative risks for each specific city vary with rollback functions. Rather than we include daily average ozone concentration as a covariate, we include daily maximum ozone concentration and daily 8-h maximum ozone concentration so that they reduce the mortality rate more in Table 3. We also conduct statistical inference with simple MLE and pooled MLE and then compare them with Bayesian inference. Both MLE and pooled MLE give a slightly higher reduction rate than one based on the Bayesian inference in Table 4. Several issues associated with high temperature on the mortality reduction. These are investigated by presenting mortality reduction with and without high temperaturein Table 5.
Several issues associated with high temperature on the mortality reduction. These are investigated by presenting mortality reduction with and without high temperature in Table 5.

Concluding Remarks
As the rollback approach provides an adjusted ozone process in a deterministic manner, we introduce a statistical rollback approach based on the quantile matching method to allow for uncertainty and thus provide flexible rollback sets. The proposed methods are applied to epidemiologic data (NMMAPs) and we assess the impact of new ozone regulation standards on public health under various settings. Another possible model is to consider different time lags in the Poisson process model with constrained or unconstrained distributed lags together. Using respiratory deaths rather than nonaccidental deaths may also provide useful results. One of the issues that need to be addressed is collinearity between ozone and temperature, for example. Other studies revealed strong effects in the Northeast and Industrial Midwest, less strong but still significant effects in the Southeast and possibly Southern California . Thus, we can also consider a regional or spatial structure in the method. Finally, the parametric rollback approach using quantile matching may allow uncertainty in ozone process itself, which can be extended to a fully Bayesian framework.