Mark Delorey F. Jay Breidt Colorado State University

1. This research is funded by U.S.EPA – Science To Achieve Results (STAR) Program Cooperative Agreement # CR - 829095 Distribution Function Estimation in Small Areas for Aquatic ResourcesSpatial Ensemble Estimates of Temporal Trendsin Acid Neutralizing Capacity Mark Delorey F. Jay Breidt Colorado State University

2. Project Funding • The work reported here was developed under the STAR Research Assistance Agreement CR-829095 awarded by the U.S. Environmental Protection Agency (EPA) to Colorado State University. This presentation has not been formally reviewed by EPA.  The views expressed here are solely those of the presenter and the STARMAP, the Program he represents. EPA does not endorse any products or commercial services mentioned in this presentation.

3. Outline • Statement of the problem: How to get a set of estimates that are good for multiple inferences of acid trends in watersheds? • Hierarchical model and Bayesian inference • Constrained Bayes estimators • adjusting the variance of the estimators • Conditional auto-regressive (CAR) model • introducing spatial correlation • Constrained Bayes with CAR • Summary

4. The Problem • Evaluation of the Clean Air Act Amendments of 1990 • examine acid neutralizing capacity (ANC) • surface waters are acidic if ANC < 0 • supply of acids from atmospheric deposition andwatershed processes exceeds buffering capacity • Temporal trends in ANC within watersheds (8-digit HUC’s) • characterize the spatial ensemble of trends • make a map, construct a histogram, plot an empirical distribution function

5. Data Set • 86 HUC’s in Mid-Atlantic Highlands • ANC in at least two years from 1993–1998 • HUC-level covariates: • area • average elevation • average slope, max slope • percents agriculture, urban, and forest • spatial coordinates • dry acid deposition from NADP

6. Region of Study

7. Locations of Sites

8. Small Area Estimation • Probability sample across region • regional-level inferences are model-free • samples are not sufficiently dense in small watersheds (HUC-8) • need to incorporate auxiliary information through model • Two standard types of small area models (Rao, 2003) • area-level: watersheds • unit-level: site within watershed

9. Two Inferential Goals • Interested in estimating individual HUC-specific slopes • Also interested in ensemble:spatially-indexed true values:spatially-indexed estimates: • subgroup analysis: what proportion of HUC’s have ANC increasing over time? • “empirical” distribution function (edf):

10. Deconvolution Approach • Treat this as measurement error problem: • Deconvolve: • parametric: assume F in parametric class • semi-parametric: assume F well-approximated within class (like splines, normal mixtures) • non-parametric: assume EF[ei] is smooth • Not so appropriate for heteroskedastic measurements, explanatory variables, two inferential goals

11. Hierarchical Area-Level Model • Extend model specification by describing parameter uncertainty: • Prior specification:

12. Bayesian Inference • Individual estimates: use posterior meanswhere • Do Bayes estimates yield a good ensemble estimate? • use edf of Bayes estimates to estimate F? • No: Bayes estimates are “over-shrunk” • too little variability to give good representation of edf (Louis 1984, Ghosh 1992)

13. Adjusted Shrinkage • Posterior means not good for both individual and ensemble estimates • Improve by reducing shrinkage • sample mean of Bayes estimates already matches posterior mean of • adjust shrinkage so that sample variance of estimates matches posterior variance of true values • Louis (1984), Ghosh (1992) • Cressie and Stern (1991)

14. Constrained Bayes Estimates • Compute the scalars • Form the constrained Bayes (CB) estimates aswhere

15. Shrinkage Comparisons for the Slope Ensemble

16. Numerical Illustration • Compare edf’s of estimates to posterior mean of F: • Comparison of ensemble estimates at selected quantiles:

17. Estimated EDF’s of the Slope Ensemble CB Posterior Mean Bayes

18. Spatial Model • Letwhere  is an unknown coefficient vector, C = (cij) represents the adjacency matrix,  is a parameter measuring spatial dependence,  is a known diagonal matrix of scaling factors for the variance in each HUC, and  is an unknown parameter. • Adjacency matrix C can reflect watershed structure

19. Conditional Auto Regressive (CAR) Model • Let Ahdenote a set of neighboring HUCs for HUC h • The previous formulation is equivalent to: • Cressie and Stern (1991)

20. HUC Structure • First level (2-digit) divides U.S. into 21 major geographic regions • Second level (4-digit) identifies area drained by a river system, closed basin, or coastal drainage area • Third level (6-digit) creates accounting units of surface drainage basins or combination of basins • Fourth level (8-digit) distinguishes parts of drainage basins and unique hydrologic features

21. Neighborhood Structure • All watersheds within the same HUC-6 region were considered part of same neighborhood • No spatial relationship among HUC-4 regions or HUC-2 regions considered at this point

22. Model Specifications • Adjacency matrix: •  = diaghh = , h = 1,…,m; nh = # neighbors of HUC hhk = 0, h ≠ k •  can be fixed or random

23. Constrained Bayes with CAR ( fixed) • If  is known, Stern and Cressie (1999) show how to solve for H1(Y) and H2(Y) under the mean and variance constraints:andrespectively, where

24. When  is Unknown or Random • We place a uniform prior on  and minimize the Lagrangian:whereto get a system of equations that can be used to solve for • Posterior quantities can be estimated using BUGS or other software

25. Spatial Structure

26. Summary • In Bayesian context, posterior means are overshrunk; in order to obtain estimates appropriate for ensemble, need to adjust • In CAR, if  is known, can find CB estimators following Stern and Cressie (1999); if  is unknown, can still find CB estimators numerically • Contour plot indicates that trend slopes of ANC are smoothed and somewhat homogenized within HUC

27. Ongoing Work • Replace spatial CAR with geostatistical model; model site responseswhere • Is CB estimate of rate the same as rate from CB estimates?

28. Other Issues • Restrict to acid-sensitive waters • Combine probability and convenience samples • Modify spatial structure