A Hierarchical Modeling Approach for Estimating ... - ACS Publications

Suite 180, Chapel Hill, North Carolina 27517, Office of. Enforcement and Compliance Assurance, U.S. Environmental. Protection Agency, 1200 Pennsylvani...
0 downloads 0 Views 124KB Size
Download PDF
Environ. Sci. Technol. 2004, 38, 1176-1182

A Hierarchical Modeling Approach for Estimating National Distributions of Chemicals in Public Drinking Water Systems S O N G S . Q I A N , * ,† A N D R E W S C H U L M A N , ‡ JONATHAN KOPLOS,§ ALISON KOTROS,§ AND PENNY KELLAR† The Cadmus Group, Inc., 6330 Quadrangle Drive, Suite 180, Chapel Hill, North Carolina 27517, Office of Enforcement and Compliance Assurance, U.S. Environmental Protection Agency, 1200 Pennsylvania Avenue NW, MS 2222A, Washington, DC 20460, and The Cadmus Group, Inc., 57 Water Street, Watertown, Massachusetts 02472

Water quality studies often include the analytical challenge of incorporating censored data and quantifying error of estimation. Many analytical methods exist for estimating distribution parameters when censored data are present. This paper presents a Bayesian-based hierarchical model for estimating the national distribution of the mean concentrations of chemicals occurring in U.S. public drinking water systems using fluoride and thallium as examples. The data used are Safe Drinking Water Act compliance monitoring data (with a significant proportion of left-censored data). The model, which assumes lognormality, was evaluated using simulated data sets generated from a series of Weibull distributions to illustrate the robustness of the model. The hierarchical model is easily implemented using the Markov chain Monte Carlo simulation method. In addition, the Bayesian method is able to quantify the uncertainty in the estimated cumulative density function. The estimated fluoride and thallium national distributions are presented. Results from this study can be used to develop prior distributions for future U.S. drinking water regulatory studies of contaminant occurrence.

1. Introduction Under the Safe Drinking Water Act (SDWA), the United States Environmental Protection Agency (EPA) develops drinking water regulations to protect public health. These regulations can be expressed as the “maximum contaminant level” (MCL), the maximum concentration of a contaminant deemed safe and allowable in drinking water. The SDWA requires that the EPA, at least once every 6 years, review existing drinking water regulations. Potential revisions to the regulations depend in part on the reevaluation of the public exposure to contaminants in drinking water as based on the extent and degree of contaminant occurrence. Currently, contaminant occurrence is measured by the frequency or number of drinking water systems in the nation with mean contaminant concentrations exceeding MCLs. This paper * Corresponding author phone: (919)403-5105; fax: (919)401-5867; e-mail: [email protected]. † The Cadmus Group, Inc., Chapel Hill, NC. ‡ U.S. Environmental Protection Agency. § The Cadmus Group, Inc., Watertown, MA. 1176

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 38, NO. 4, 2004

describes the hierarchical modeling approach developed in support of the EPA’s first Six-Year Review of National Primary Drinking Water Regulations. The model can be used to generate estimates of the national distribution of system mean contaminant concentrations using a very large number of existing compliance monitoring data from nearly all drinking water systems in 16 states. An important feature of the model is its capability to handle censored data. In this section, we briefly discuss some statistical methods that have been used by the EPA to handle censored data. This is followed by a presentation of the data collection process for the Six-Year Review in Section 2.1. The statistical model is presented in Section 2.2. The paper concludes with a presentation of occurrence analysis results from two example chemicals and a discussion in Section 3. Assessing contaminant occurrence in drinking water is difficult because concentrations are usually low, and monitoring results often include concentration values below the method detection limit (MDL). When a concentration value is below the MDL, the exact concentration value is only known to be less than a specific value. Data points with this characteristic are referred to as “nondetections” or “left censored”. When the underlying probability distribution of the contaminant concentration is of interest, a censored value contributes less information than an exactly measured value does. Nevertheless, the censored value has information that should not be disregarded when estimating the probabilistic distribution parameters. There are many definitions of detection or quantification limits (e.g., ref 1). In this paper, we use the lowest analytical concentration (or level) that can be reliably achieved within specified limits of precision and accuracy under routine laboratory operating conditions. This concentration is often referred to as the Minimum Reporting Level (MRL). Many statistical methods are available for estimating distribution parameters when censored data are present (e.g., ref 2). The simple substitution method replaces the leftcensored value by a specific number (e.g., 0, one-half of the MRL value, or the MRL value), which may lead to biased estimates. Alternatively, a regression on ordered statistics (ROS) method for estimating the mean and variance of a log-normal distribution has been proposed (3, 4). Suppose we have observations from a concentration variable C: (c1, c2,‚‚‚, cn), in which k observations are known to be less than the MRL. Let c1, c2,‚‚‚, ck be censored and ck+1,‚‚‚, cn be accurately observed. Assuming C ∼ LN(µ,σ2) and letting Y ) {yi ) log(ci), i ) 1,‚‚‚, n}, we have Y ∼ N(µ,σ2). The percentiles of the uncensored log concentration data {yi, i ) k + 1,‚‚‚, n} can be estimated unambiguously. For a given percentile p, yp ) µ + σxp, where yp is the p percentile of Y and xp is the p percentile of the standard normal distribution. Consequently, by fitting a simple linear regression of yi (i ) k + 1,‚‚‚, n) against their corresponding standard norm quantiles, the resulting slope and intercept are the estimates of σ and µ, respectively (5). Because ROS is a rank based method, its validity may be affected when ties are present. Another frequently used method for estimating parameters of normal or other known distributions is the maximum likelihood estimator (MLE). We divide the complete data Y into the observed (Yobs) and the censored (Ycen) parts, i.e., Y ) (Yobs, Ycen). The likelihood function is defined as L(θ| Yobs) ∫Rf (Yobs, Ycen|θ)dYcen, where R is the censoring region, f is a density function, and θ is the distribution parameter (6). The MLE of the distribution parameters (µ and σ2) is the estimator that maximizes the likelihood function. However, the presence of censored values makes computation of the 10.1021/es020686q CCC: $27.50

 2004 American Chemical Society Published on Web 01/06/2004

MLE more difficult. A commonly used computational method is the Expectation-Maximization (EM) algorithm (6, 7). Because both the ROS and MLE methods assume observations are independent random samples from the same distribution and the contaminant concentration distribution varies from system to system, it is necessary to apply these two methods separately for each drinking water system. Those separately estimated system means are then pooled to form the national distribution of system means. For example, ROS was used in the EPA’s arsenic analysis (8) to calculate means of arsenic concentrations in drinking water for individual public water systems. The resulting system means were grouped by state and a probability of system means exceeding the MCL was estimated for each state. Weighted sums of the exceedance probabilities generated for each state were then used to develop regional and national exceedance probabilities. These exceedance probabilities were used to estimate the numbers of systems exceeding specific arsenic concentration thresholds. MLE was used in the EPA’s radon analysis (9). Because of the large numbers of drinking water systems and contaminants to be reviewed under the Six-Year Review, these two methods would be prohibitively expensive. Furthermore, distribution parameters cannot be adequately estimated for systems with a large number of left-censored values. As a result, those systems would be excluded, likely resulting in overestimation of the occurrence. In the applications of the ROS and MLE methods, only the estimated system means are used to assess the regional and national exceedance probabilities. Uncertainty of the estimated system means varies from system to system due largely to the difference in sample size. Treating these estimated means with equal weight may lead to misleading estimates of exceedance probabilities. For the first Six-Year Review, we used a Bayesian hierarchical modeling approach (10). The Bayesian approach was initially motivated for its ease in computation. The national distribution of system mean contaminant concentration was not directly observed. By using a hierarchical structure, individual system means are linked to form the national distribution of system means. This hierarchical structure facilitates our computational strategy, which addresses the problem of uneven sample sizes and provides a Bayesian framework for combining and updating information. Although it is possible to use the conventional MLE approach in this study, we found that a Bayesian approach has advantages over the MLE approach in addition to the computational simplicity. A fundamental difference between the Bayesian approach and the MLE approach is that the Bayesian approach makes inference over the entire support of the posterior distribution, while the MLE approach makes inference based on a point estimate. A Bayesian method does not depend on the assumption of asymptotic normality underlying the MLE method. Using recently developed computer-intensive sampling-based methods, we can estimate the full distributional profile of the posterior distributions of system mean concentrations and convert them into uncertainty in exceedance probabilities without further assuming log-normality. Furthermore, the Bayesian approach provides a framework for accumulating information because the current Six-Year Review posterior distributions can be used as priors for the next Six-Year Review. Although developed independently, our model is similar to the model of Lockwood et al. (11). The main difference between the two models relates to the way the among-system dependency was modeled (see Section 2.2).

2 Methods 2.1. Data Collection and Management. The first Six-Year Review assessed 60 contaminants, including regulated

TABLE 1. Contaminants Analyzed for All Systems in Support of the Six-Year Review (MCL in mg/L) Inorganic Chemicals Inorganic Chemicals antimony (0.006) beryllium (0.004) chromium (0.1) fluoride (4) selenium (0.05)

barium (2) cadmium (0.005) cyanide (0.2) mercury (0.002) thallium (0.002)

Synthetic Organic Chemicals alachlor (0.002) benzo[a]pyrene (0.0002) bis(2-ethylhexyl)phthalate (0.006) chlordane (0.002) dalapon (0.2)

atrazine (0.003) bis(2-ethylhexyl)adipate (0.4) carbofuran (0.04) 2,4-D (0.07) 1,2-dibromo-3-chloropropane (0.0002) diquat (0.02) endrin (0.002) glyphosate (0.7)

dinoseb (0.007) endothall (0.1) ethylene dibromide (0.00005) heptachlor (0.0004) heptachlor epoxide (0.0002) hexachlorobenzene hexachlorocyclopentadiene (0.001) (0.05) lindane (0.0002) methoxychlor (0.04) oxamyl (0.2) PCBs (0.0005) pentachlorophenol (0.001) picloram (0.5) simazine (0.004) toxaphene (0.003) 2,4,5-TP (0.05) Volatile Organic Chemicals benzene (0.005) 1,4-dichlorobenzene (0.075) 1,2-dichloroethane (0.005) cis-dichloroethylene (0.07) dichloromethane (0.005) ethylbenzene (0.7) styrene (0.1) toluene (1) 1,1,1-trichloroethane (0.2) trichloroethylene (0.005) xylenes (total) (10)

carbon tetrachloride (0.005) o-dichlorobenzene (0.6) 1,1-dichloroethylene (0.007) trans-dichloroethylene (0.1) 1,2-dichloropropane (0.005) monochlorobenzene (0.1) tetrachloroethylene (0.005) 1,2,4-trichlorobenzene (0.07) 1,1,2-trichloroethane (0.005) vinyl chloride (0.002)

inorganic chemicals, synthetic organic chemicals, and volatile organic chemicals (Table 1). Arsenic was evaluated separately through the new arsenic rule (also see ref 11), and dioxin, asbestos, and nickel were not included because there were no or little data available for a thorough national assessment. Complete analytical records of contaminants in drinking water from public water systems (PWSs) collected under SDWA are unavailable in a form that can be processed for a comprehensive national overview of occurrence and exposure. The EPA’s Safe Drinking Water Information System/Federal Version (SDWIS/Fed) maintains a variety of water system inventory and operation information as well as compliance program information. However, for most contaminants, the only analytical results filed in SDWIS/Fed are those related to drinking water regulation violations. The detailed analytical results of most contaminants are not fully reported to SDWIS/Fed and are stored in individual state databases. As a result, state data sets (comprising SDWA compliance monitoring data from PWSs) constitute the data sources for the Six-Year Review. It was not feasible to collect and manage the compliance monitoring data for all states and territories due to many factors. A random sample may include states that are unwilling to provide data because of no legal mandate or states that have no data because they received waivers for monitoring certain pollutants based on previous evidence of low or zero occurrence. In a previous study, the EPA reviewed the occurrence of regulated contaminants in PWSs using data sets voluntarily VOL. 38, NO. 4, 2004 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

1177

FIGURE 1. Distribution of state ranking for manufacturing establishments per square mile versus farm agriculture chemical expenses (+: 16 cross-section states, o: all other states). provided by 12 states (12). Data from eight of these 12 states were selected for use in the development of an initial analysis of a national cross-section of contaminant occurrence. The results of this prior work, and the report generated (12), have been widely reviewed by the EPA stakeholders and were also peer-reviewed. For the current study, the EPA used the same method to select an additional eight states, resulting in a 16-state cross-section. The cross-section development was based on researched “pollution potential indicators”, published quantitative measures that serve as indicators of the potential pollution from manufacturing and agriculture in each state. States were then ranked from 1 to 50 for each indicator, and the resulting rankings were divided into quartiles. States to be selected for the national cross-section were to be distributed in a balanced manner across all quartiles. Availability of complete state data sets and a balanced geographic distribution of states for selection were also considered. (Consideration of states that together provide a geographic diversity is one means to include contaminant occurrence data from a wide range of climatic and hydrologic conditions across the United States.) The 50 states were plotted on a two-dimensional scatter plot, with the x- and y-axes representing their respective manufacturing and agricultural ranking (Figure 1). We note that the 16 states are not randomly selected. Consequently, the data we used may be biased. However, results based on the cross-section should be indicative of national occurrence based on the broad range of pollution potential indicators, as well as geographic diversity. The data sets used in this study are SDWA compliance monitoring data from the 16 states. There are more than 13 million individual concentration values from approximately 41 000 PWSs. For the two chemicals we used as examples, there are approximately 93 000 concentration values from 21 000 PWSs for fluoride and 47 000 observations from 18 000 PWSs for thallium. For a complete description of the 16state cross-section development, see ref 13. The PWSs were divided into 10 groups according to their source water type 1178

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 38, NO. 4, 2004

(groundwater or surface water) and population served (e500, 501-3300, 3301-10 000, 10 001-50 000, and >50 000), a common stratification used by the EPA. Data included are those collected from 1993 to 1997. There are multiple MRLs for each pollutant reflecting the difference among states as well as improvements of analytical methods made over time. 2.2. Bayesian Hierarchical Modeling. Bayesian methods are currently experiencing an increasing popularity in the sciences as a means of probabilistic inference (e.g., ref 14). Among their advantages are the ability to incorporate prior information, the ease of incorporating into a formal decision analytic context, the explicit handling of uncertainty, and the straightforward ability to assimilate new information in contexts such as adaptive management. Introduction to Bayesian statistics and data analysis can be found in refs 10, 15, and 16. All drinking water systems in the United States are subject to the same federal regulations. Because it is difficult to collect monitoring data from all PWSs in the United States, the national distribution of system mean contaminant concentrations required by SDWA should be treated as the probability distribution that governs the mean contaminant concentrations of each system. We can imagine that individual log concentration values from a given system are generated from a normal distribution with an unknown system mean and variance. The system means in a given group, in turn, are generated from a higher level distribution with an unknown group mean and variance. The group means are further generated from a super distribution. In other words, the observed log concentration data are modeled conditionally on their respective system mean and variance, which themselves are given a probabilistic specification using stratum mean and variance, and the strata means are modeled as from a common population distribution. An important feature of the hierarchical model is that the hierarchical probability distribution imposes dependence into parameters, which allows a hierarchical model to have enough parameters to form a realistic model without

overfitting the data (10). Details of the Bayesian hierarchical modeling approach can be found in refs 10 and 17. The model is summarized in eq 1

TABLE 2. Summary Statistics of Generated System Means µ and Standard Deviations σ

yijk ∼ N(yijk|µij, σ21)I(, Sijk) µij ∼

µij σ1

N(µij|µi, σ22)

µi ∼ N(µi|µ,σ23)

(1)

where yijk is the kth observed log concentration value from system j, in group i (defined by source water type and population range), µij is the system mean, and σ21 is the variance. µi is the group mean, and σ22 is the group variance of system means, while µ and σ23 are the hypermean and variance of the group means. Note that the hierarchical notation in eq 1 indicates conditional distributions, that is, yijk is normal conditional on µij and σ21, µij is normal conditional on µi and σ22, and µi is normal conditional on µ and σ23. The notation I(, Sijk) indicates that the corresponding concentration value is less than Sijk or left-censored. Numerically, when evaluating the value of a censored data point, a random variate is generated using the current guess of the model parameters (from the MCMC process). The censored value is replaced by the first generated value that is less than Sijk. A new set of model parameter estimates are produced using uncensored data as well as previously generated censored values. The newly estimated model parameters are used to generate a set of new values for the censored data (see discussions about BUGS implementation of interval censorship in ref 18). The prior distributions of µ, σ21, σ22, and σ23 are

µ ∼ N(µ| 0, 10 000) 1/σ2r ∼ gamma(1/σ2r| 0.001, 0.001), r ) 1, 2, 3

(2)

where N(0, 10 000) is the normal distribution with mean 0 and variance 10000 and gamma(0.001, 0.001) is the gamma distribution with shape and scale parameters of 0.001. The prior distributions for 1/σ21, 1/σ22, 1/σ23, and µ are considered “noninformative” or vague. The distribution N(0, 10 000) has a mean 0 and a standard deviation of 100. The width of the 95% prior interval is approximately (200, that is, it is practically “flat” in the region of interest. The standard noninformative prior for a variance parameter is p(σ2)∝1/σ2, which arises from assuming that the log of the variance parameter has a uniform distribution on (- ∞, + ∞). This prior is improper, which may lead to an improper posterior distribution. The vague gamma prior, gamma(0.001, 0.001), is instead commonly used. See refs 17 and 18 for details. We used noninformative priors because we have little prior knowledge about these parameters for the 60 chemicals, and the data have a large sample size which would dominate the posterior for any reasonable priors. For the next Six-Year Review, the posterior distributions of these parameters will be used as priors. A Markov chain Monte Carlo simulation (MCMC) method (19) is used for simultaneously estimating the distribution parameters by sampling the parameters from their joint posterior distribution. The MCMC method is implemented using the freely available WinBUGS software (18). The national distribution as well as the probabilities of exceeding certain thresholds is estimated in the same Monte Carlo simulation. The MCMC method allows us to sample µij and σ21 from their joint posterior distribution. All inferences are made based on these posterior samples. The burn-in period required is less than 100 based on the Raftery and Lewis

min

25%

-6.717 8.26 × 10-4

-1.931 0.096

µij -20.91 σ1 0.752

median

mean

Simulated Fluoride -1.148 0.210

75%

max

1.145 -0.3617 3.175 0.411 0.462 16.58

Simulated Thallium

-11.88 -10.05 -10.06 -8.239 1.074 1.168 1.185 1.277

3.581 2.136

procedure (22, 23). To ensure convergence and independence, we set the burn-in period to be 1000 and took 500 samples for each of the unknown quantities from the next 1000 MCMC iterations (one sample every other iteration). Because we assume the log concentration values follow a normal distribution at the system level, the arithmetic system mean concentration, mnij, is calculated from the generated log mean µij and log variance σ21, i.e., mnij ) exp[µij + 0.5 × σ21]. For each set of system means, an empirical cumulative distribution function (CDF) is generated as an estimate of the national distribution. From this empirical CDF we calculate the proportions of systems with mean concentrations exceeding certain thresholds as estimates of the exceedance probabilities. By repeated sampling of the system means, we have many empirical CDFs of the national distribution and many estimates of the exceedance probabilities. These empirical CDFs are used to summarize the national distribution as well as the uncertainty about the distribution. Because information such as source water type (groundwater or surface water) and the population served are available for each of the PWSs, the national distribution of system mean concentrations can be used to make population-based inferences, such as the proportions of population served by various types of PWSs with elevated levels of contaminants. Note that the estimated means and variances are based on log-transformed concentration values. Thus, it is necessary to obtain the system means in the original metric. The proportion of the estimated system means above a threshold value is an estimate of the probability of exceeding the threshold. If the MCMC process produces 500 pairs of samples of µij and σ21, we have 500 empirical CDFs of the national distribution, 500 samples of the arithmetic system mean, and hence, 500 estimates of the exceedence probability for each threshold. From these 500 probabilities, we can calculate the mean, median, and the 90% credible interval to summarize the uncertainty of the quantity. Because the exceedance probabilities are almost always very small, their variances tend to be very small. In addition, the empirical CDF is evaluated using a large number of system means (usually > > 15 000), resulting in a stable CDF. Based on initial tests on selected chemicals, a sample size of 500 is sufficient for evaluating exceedance probabilities. It is reasonable to assume that system means are correlated. Lockwood et al. (11) modeled the correlation using a distance-based covariance model. Since their subject was source (untreated) water arsenic concentrations, a distancebased approach is sensible. For the finished (treated) drinking water in our case, potential pollution in a PWS is more related to the population served (the size) as well as the source water type (groundwater versus surface water). We divided systems into 10 groups based on source water type and service population size, and the correlations among systems were indirectly modeled through the hierarchical structure. Specifically, system means µijs within a stratum are connected by the population distribution of system means for the same stratum, while the stratum means µi are connected by their population distribution. In other words, the correlations VOL. 38, NO. 4, 2004 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

1179

FIGURE 2. Results of a simulated fluoride data set with known (“true”) CDF.

FIGURE 3. Results of a simulated thallium data set with known (“true”) CDF. among the µij’s were implicitly modeled through the hierarchical model. This approach simplifies computation and allows us to use noninformative priors. We note that we assumed an a´ priori constant withinsystem variance (σ21), which is necessary because many systems in the data set have only one observation. However, this prior constant variance assumption will not result in a constant posterior variance for all systems. This is because the mean µij is modeled as a random variable. In other words, the posterior predictive distribution of a future observation from system j in the ith group Y˜ ij is estimated by the conditional distribution p(Y˜ ij|µij,σ21, Y). The variance of Y˜ ij is a function of both µij and σ21, and hence, is system specific. 2.3. Model Evaluation. The model presented in eqs 1 and 2 was evaluated using a simulation study. We generated seven synthetic data sets, each containing the same number of observations from the same number of PWSs as contained in the actual 16-state fluoride or thallium data sets. We present the results from the two data sets with within-system observations generated from a Weibull distribution. The Weibull distribution is a light tail distribution, while the lognormal is a heavy tail distribution (see ref 20, pp 69-70). Using a Weibull distribution would result in a very different extreme value pattern than the log-normal distribution would, allowing an evaluation of the robustness of the model 1180

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 38, NO. 4, 2004

against the log-normal assumption at the system level. For each system, a system mean is generated from a log-normal distribution: LN(-1.15, 1.33) (log-mean and log-variance) for the simulated fluoride data and LN(-10, 5.76) for simulated thallium data. The system variances are generated from the inverse gamma distribution with parameters (20, 9) for fluoride and from the inverse-gamma (15, 20) for thallium. The generated system means (µij) and standard deviations (σ1) are summarized in Table 2. Observations within a PWS are generated from a Weibull distribution with a shape parameter of 2. The scale parameter of the Weibull distribution is calculated for each system to yield the mean and variance assigned for the system. The generated data were censored at 0.11 mg/L for the simulated fluoride and 0.001 mg/L for the simulated thallium to produce the similar fractions of censored values in the simulated data as in the 16-state data sets (∼25% and ∼96%, respectively). The Bayesian hierarchical model-estimated national distributions are almost identical to the true distributions (Figures 2 and 3). When the log-normal distribution was used at the system level, the results (not shown) were similar to those presented in Figures 2 and 3. The estimated national distribution of system mean fluoride concentrations is presented using a CDF (Figure 4). In the figure, the black solid line is the median, and the black

FIGURE 4. The CDF of the national distribution of system mean fluoride concentrations using the 16-state data set.

TABLE 3. Selected Percentiles (mg/L) of the Estimated National Distribution of System Mean Fluoride Concentrations percentile 1-line 0.5-line 0-line model

5%

10%

25%

50%

75%

90%

95%

0 0.05 0.10 0.05

0 0.05 0.10 0.07

0.08 0.10 0.13 0.13

0.21 0.23 0.24 0.27

0.59 0.60 0.60 0.60

1.10 1.10 1.10 1.11

1.68 1.68 1.68 1.67

dotted lines represent the 95% credible interval (which is very narrow due to the large sample size of system means). The estimated CDF is compared to three empirical CDFs calculated from system means where all censored values have been substituted with (1) 0 (the 0-line), (2) half of the MRL values (the 0.5-line), and (3) the MRL values (the 1-line). This graphical comparison indicates that the proposed model yields desired results in that the four CDFs converge at higher concentration values where left censoring is less likely to occur. In addition, we note that most of the fluoride MRL values are 0.1 mg/L, which is reflected in the figure as the widest gap between the model-estimated CDF and the 1-line. The second most common MRL value is 0.2 mg/L, where the 0.5-line diverted slightly from the model estimated CDF. As expected, the estimated CDF is bounded by the 0- and 1-lines

and is close to the 0.5-line. Because the censoring occurs primarily at the lower end of the distribution, the differences among the four CDFs appear at the left tail of the distribution. Table 3 also presents a comparison of selected percentiles of the estimated fluoride national distribution, based on the predictions from the model and the three empirical CDFs evaluated from the substitutions. As expected, the difference between the model-predicted and the empirical percentiles are principally in the lower half of the distribution.

3. Results and Discussion We present the results for fluoride and thallium as examples. The estimated CDF of national distributions (Figures 4 and 5) represents the fraction of systems with a mean concentration less than the given concentration threshold. By separating the generated system means into the 10 groups described earlier, we are able to estimate the fraction of systems with mean concentrations exceeding selected thresholds as well as the fraction of the population within each group exposed to concentrations higher than the selected thresholds. These results, along with results for all 60 chemicals listed in Table 1 and the WinBUGS programs, are presented elsewhere (13). The total number of community water systems plus nontransient, noncommunity water systems listed in an EPA

FIGURE 5. The CDF of the national distribution of system mean thallium concentrations using the 16-state data set. VOL. 38, NO. 4, 2004 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

1181

handbook is 65 030, serving a population of 213 008 182 (21). (The handbook also presents the system and populationserved numbers stratified by source water type and population-served size categories.) This information can be used to estimate the total number of systems nationwide (and the number of people they serve) having mean concentrations exceeding selected thresholds, assuming the 16-state crosssection provides a representative national sample. For example, the expected number of systems nationally expected to have a mean fluoride concentration exceeding 4 mg/L for groundwater systems serving 500 people or less is 256. This estimate is derived by multiplying the estimate of probability of exceedance (0.590%) by the total number of groundwater systems nationally that serve 500 people or less (43 498). Using the population-based probability in the same group, about 46 000 (from a total of 6 498 708) people are expected to be exposed to a mean fluoride concentration above 4 mg/L. This stratified occurrence information can be used to support some aspects of a regulatory impact assessment (e.g., supporting a cost-benefit analysis that might be used to determine whether a new regulatory MCL should be imposed). For example, all PWSs are currently required to maintain a mean fluoride concentration less than 4 mg/L. A total of about 332 (284-385) systems are estimated to exceed this threshold, and most of them (256 or 212-299) are groundwater systems serving less than 500 people. There are about 191 000 (118 200-301 000) people receiving water with a mean fluoride concentration above 4 mg/L, with most of them (178 000 or 112 000-268 900) receiving their water supply from groundwater PWSs. If PWSs are required to maintain a mean fluoride concentration below 2 mg/L, the affected PWSs are expected to be 1885 (still primarily small groundwater systems), but the new requirement will benefit about 2 million people, and a substantial number of them will be in areas served by large surface water PWSs. The use of the hierarchical modeling approach in this study has at least three advantages. First, the model is robust against departures from the log-normal assumption. This can be explained on two levels: (1) The national distribution is expressed as a mixture of many log-normal distributions, and the resulting national distribution is not subject to the log-normality assumption. This flexibility makes the model prediction more realistic. (2) At the system level, log-normality is conditional on the system mean, and the system mean is modeled as a random variable. As a consequence, the log-normal assumption at the system level is also relaxed. The second advantage of the hierarchical approach is the relatively light computational burden. Less computation is needed because the hierarchical model produces the national distribution in one model run, while for the ROS plotting position method to be used properly, mean concentration distributions for each PWS have to be estimated separately. The third advantage of using the hierarchical modeling approach is its Bayesian feature, which allows combining information from all systems. This feature is desired because not all the water supply systems have the same amount of data, and for some systems all measurements are below the MRL. As a result, if the national distribution were estimated based on individual system distributions as in the ROS plotting position method, uncertainty in estimated individual system distributions would vary significantly. For future SixYear Reviews, results from this study will be used as priors,

1182

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 38, NO. 4, 2004

providing a framework for systematic information accumulation.

Acknowledgments This work was partially funded by the U.S. EPA (Contract # 68-C-99-206). The authors thank Walter Grayman, Elizabeth Margoches, Michael Messner, Clifton Sutton, and Inge Werner for constructive comments and suggestions on an earlier draft of this paper. Comments and suggestions from the associate editor and three reviewers are greatly appreciated. The views expressed in this paper are those of the authors and do not necessarily reflect the views or policies of the United States Environmental Protection Agency.

Literature Cited (1) Scroggin, D. G. Hazard. Waste Hazard. Mater. 1994, 11(1), 1-4. (2) Gleit, A. Environ. Sci. Technol. 1985, 19, 1201-1206. (3) Gilliom, R. J.; Helsel, D. R. Water Resour. Res. 1986, 22, 135146. (4) Helsel, D. R.; Gilliom, R. J. Water Resour. Res. 1986, 22, 147155. (5) Helsel, D. R.; Cohn, T. A. Water Resour. Res. 1988, 24(12), 19972004. (6) Little, R. J. A.; Rubin, D. B. Statistical Analysis with Missing Data; Wiley: New York, 1987. (7) Tanner, M. A. Tools for Statistical Inference: Methods for the Exploration of Posterior Distributions and Likelihood Functions; Springer-Verlag: New York, 1996. (8) USEPA. Arsenic Occurrence in Public Drinking Water Supplies; EPA Office of Ground Water and Drinking Water: 2000a. (9) USEPA. Methods, Occurrence, and Monitoring Document for Radon in Drinking Water; EPA Office of Ground Water and Drinking Water: 2000b. (10) Gelman, A.; Carlin, J. B.; Stern, H. S.; Rubin, D. B. Bayesian Data Analysis; Chapman & Hall: London, 1995. (11) Lockwood, J. R.; Schervish, M. J, Gurian, P.; Small, M. J. J. Am. Stat. Assoc. 2001, 96(456), 1184-1193. (12) USEPA. A Review of Contaminant Occurrence in Public Water Systems; EPA816-R-99-006; U.S. EPA, Office of Ground Water and Drinking Water: 1999. (13) USEPA. Occurrence Estimation Methodology and Occurrence Findings Report for the Six-Year Review of Existing National Primary Drinking Water Regulations; EPA-815-R-03-006; Office of Water: June 2003, 874p. http://www.epa.gov/safewater/ standard/review/pdfs/support•6yr•occurrencemethods• final.pdf (14) Malakoff, D. M. Science 1999, 286, 1460-1464. (15) Box, G. E. P.; Tiao, G. C. Bayesian Inference in Statistical Analysis; Addison-Wesley: Reading, MA, 1973. (16) Bernardo, J. M.; Smith, A. F. M. Bayesian Theory; Wiley: West Sussex, England, 1994. (17) Congdon, P. Bayesian Statistical Modelling; Wiley: West Sussex, England, 2001. (18) Spiegelhalter, D. J.; Thomas, A.; Best, N.; Gilks, W. R. BUGS 0.5: Bayesian Inference Using Gibbs Sampling Manual; Medical Research Council Biostatistics Unit, Institute of Public Health: Cambridge, UK, 1996. (19) Gilks, W. R.; Richardson, S.; Spiegelhalter, D. J. Markov Chain Monte Carlo in Practice; Chapman & Hall: London, 1996. (20) Ott, W. R. Environmental Statistics and Data Analysis; Lewis Publishers: Boca Raton, FL, 1995. (21) USEPA. Water Industry Baseline Handbook, 2nd ed. (draft); 2000c; March 17, 2000. (22) Raftery, A. E.; Lewis, S. M. Stat. Sci. 1992, 7, 493-497. (23) Raftery, A. E.; Lewis, S. M. The number of iterations, convergence diagnostics and generic Metropolis algorithms. In Practical Markov Chain Monte Carlo; Gilks, W. R., Spiegelhalter, D. J., Richardson, S., Eds.; Chapman and Hall: London, U.K., 1995.

Received for review April 10, 2002. Revised manuscript received November 18, 2003. Accepted November 21, 2003. ES020686Q