# Presentation: 804896

## Bayesian Methods for Regional Eutrophication Models

- E. Conrad Lamon III
- Dept. of Environmental Studies
- Louisiana State University
- and
- Craig A. Stow
- Department of Environmental Health Sciences
- University of South Carolina

- STAR Progress Review Workshop
-    Old Town Alexandria, VA
- June 16-18, 2004

**Notes:**
As you can see, the goal of this research is to develop regional models.  

More specifically, the question is:
“What are the regions within which Lake Eutrophication models work”?

Among other benefits, regional models have the advantage of allowing inference regarding unmonitored sites.

Bayesian Methods were chosen due to the nature of the question. 
1st it regards uncertainty in defining regions and ,then, 

given the regions, what are the models?

So, model structure is hierarchical in nature.

## Overview

- Goals and Objectives
- Approach
- Preliminary Findings
- Significance
- Next Steps

**Notes:**
As instructed, I’ll follow this outline....

## Goals and Objectives

- Use modern classification and regression trees and hierarchical Bayesian techniques to link multiple environmental stressors to biological responses and quantify uncertainty in model predictions and parameters.

**Notes:**
Stressor response models, and more generally, model based inference may improve upon the current approach to Nutrient Criteria development, based (simply) on the assumption that the 25th percentile of nutrient levels represent an un-impacted state.  

Instead of y = ybar as our (implied) model, assume some variability in eutrophication response Y is explained by stressor X.  This will allow determination of nutrient criteria based on expected impact on Y.

With this statement of goals, we are compelled to comply with the recommendations of the National Research Council (NRC 2001) , namely:........

## Guidance for TMDL model selection (NRC 2001) 

- report prediction uncertainty
- be consistent with the amount of data available
- flexible enough to permit updates and improvements

**Notes:**
1) report prediction uncertainty, 2) be consistent with the amount of data available and 3) be flexible enough to permit updates and improvements. 

(The Bayesian paradigm fills the bill.)

## Approach

## Tree based methods

- are a flexible approach useful for variable subset selection
- when the analyst suspects global non-linearity
- and cannot (or does not want to) specify the functional form of  possible interactions *a priori*.

- Approach

**Notes:**
We are using several tree based methods in this research.

Generally 1, 2, 3.....

## Example tree

- Approach

**Notes:**
This is a figure from another paper to illustrate regression trees (CART).  Numbers in each node are the mean value of y within that node, while numbers below each node are the deviance of y (analogous to variance).

A split is chosen at each node that maximizes  D, the difference between that nodes’ deviance (variance) and the sum of the deviances of the two child nodes.  Each split chosen in a CART regression tree model maximizes the reduction in deviance.  This algorithm is often referred to as “greedy” in the statistical literature, because D is calculated considering only the current node and the resulting child nodes.  Greedy because there may be a split of the current node that yields a less than maximum D, but provides a larger delta D in subsequent splits.

At each node : regression trees from the CART procedure produce subsets on the design matrix X that result in maximum reduction of deviance (variance) in the two resulting child nodes.  This yields subsets of X that have relatively homogeneous Y values.  If X is restricted to geographic variables such as Lat and Lon, these subsets represent regions with homogeneous Y values (“small” variability about the mean). 

## Methods

- _C_lassification _A_nd _R_egression _T_rees (CART),
- it’s Bayesian analogue, BCART
- a recently developed enhancement to the BCART procedure, which includes BCART as a model subclass, known as Bayesian Treed (BTREED) models, and
- Bayesian Hierarchical Models

- Approach

**Notes:**
Wild for acronyms!!!

## BCART and BTREED Models

- Will be used with the EPA Nutrient Criteria Database to identify and estimate regional eutrophication stressor – response models for EPA STAR funded research.
- Lamon and Stow, 2004, _Water Research_, 38(11): 2764-2774*.*

- Approach

**Notes:**
Of course, I then pursued some Bayesian Analogues to the tree based CART procedures.  These methods also take advantage of MCMC procedures to estimate Bayesian Hierarchical models.  There will be two stages in the initial hierarchical models explored (T and theta)

This talk is based on USEPA NES data, as proof of concept.

Further work will use the EPA Nutrient Criteria Database

## Bayesian Treed models

- Bayesian Hierarchical model to:
- Select subsets on *X*     *X**s*
- Fit linear models to these subsets *X**s*
- Tree structured models
- “ANOVA in Reverse”
- “Leaves” contain linear models, not just a mean (like in CART models)

- Approach

**Notes:**
Today I’m reporting results obtained using BTREED, which are Bayesian Hierarchical models that select subsets of observations on partitions of X, for which linear models are estimated.

Further they identify subsets of X for which LM performance is improved as measured by predictive LL.

It’s helpful to contrast with ANOVA to understand the CART and BTREED processes.

With ANOVA, we subset the data along the axes of X, then seek and answer to the question “Do the subset means differ from each other???”

CART and other tree structured models seek subsets for which the means, (or variance about the means), differ.

With BTM we take this a step further by seeking an answer to the question 
“What are the subsets for which LM’s differ?”

## Bayesian Treed model specification

* **y*|*x*, with *x* = (*x*1, *x*2, ..., *x*p),
- where *p* = number of predictor variables.
**two components of model**
- 1. tree *T* with *b* bottom nodes,
- 2. parameter vector ** = (**1, **2,..., **b),
- where **i is associated with the *i*th bottom node.  If *x* is in the *i*th node, then *y*|*x* = *f*(*y*|**i), where *f* is a parametric family indexed by **i.

- Approach

**Notes:**

1. bottom nodes = (subsets),
2. theta specifies the parameters of a LM
This structure consists of two components, a tree T with b bottom nodes, and a parameter vector  = (1, 2,..., b), where i is associated with the ith bottom node.  If x is in the ith node, then y|x = f(y|i), where f is a parametric family indexed by i.  The tree is fully specified by (, T), and the Bayesian approach requires that we specify a prior, p(, T).  

A spatial level of the hierarchy....in which regions will be randomly chosen from relatively non-informative prior distributions.
	The second level of the hierarchy.....will consist of linear models fit to data in the previously selected regions.

## Bayesian Treed model specification (cont.)

- Tree is fully specified by **(**********, *****T*****)**
- need a prior,
- 			***p*****(**********, *****T*****).**
- Because ****** indexes a parametric model for each ***T***, we can use Bayes theorem such that
*	**	****p*****(**********, *****T*****) = *****p*****(********** | *****T*****)*****p*****(*****T*****)**.
- So, specify _prior_ in two stages:
- 1 – on the tree space, ***p(T)****,* and
- 2 – on the distribution of ***Y*** at the bottom nodes, conditional on ***T*****, *****p(****************T)*****.**

- Approach

**Notes:**
Because  indexes a parametric model for each T, we can use Bayes theorem such that p(, T) = p( | T)p(T).  Consequently, prior specification is done in two stages: first by specifying a prior on the tree space, p(T), and a prior on the distribution at the bottom nodes, conditional on T, p(T).  The same approach may be used for p(T) in both classification and regression trees, while the parametric specification of p(, T) will differ.

## Bayesian Treed model search

- MCMC used to stochastically search for high posterior probability trees *T*.
- Metropolis –Hastings algorithm simulates a Markov chain with limiting distribution *p(T|Y,X)*
- Chipman, George and McColloch, 2000, JASA.
- http://gsbwww.uchicago.edu/fac/robert.mcculloch/research/papers/index.html

- Approach

**Notes:**
Bayes theorem is then used with the data to update the probability of selecting that region in subsequent random draws, and 

The process is repeated many, many times,

“Bayes theorem is a data-driven probability transformation machine which maps prior distributions to posterior distributions”.  J.M. Bernardo, Bayesian Statistics, 2003.

## Data

- Response variables may be
- either continuous (such as biological indices of abundance) or
- discrete (such as designated use attainment classes).
- EPA NES example: response variable is lake-wide, summer average log10 Chlorophyll *a* concentration.

- Approach

**Notes:**
HERE, response variable is lake-wide, summer average log10 Chlorophyll a concentration from the USEPA NES dataset.

Continuous bio. Index on algal abundance.

## Data

- Predictor variables in tree based methods may also be continuous or discrete, and may include :
- source agency, basin, sub-watersheds, states, EPA regions, latitude and longitude, and many continuous predictors related to water chemistry, water use, discharges or pollutant loading.

- Approach

## Data

- For the EPA NES example, Latitude and Longitude were used in the tree portion,
- and
- 		log10Qin,  	log10 Z       Log10w
- 		In-lake log10 TP  	In-lake log10 TN
- For the linear model within each bottom node (leaf)

- Approach

**Notes:**
HERE, predictors are limited to Latitude and longitude for the tree portion, and 
		log10Qin
		log10 Z
		Log10 Tw
		Inlake log10 TP 
		in lake log10 TN
For the linear model within each bottom node, or leaf of the tree portion of the hierarchy.

## Preliminary Findings

- Lamon, E.C., and C.A. Stow, 2004. Bayesian Methods for regional-scale eutrophication models,
- _Water Research_, 38(11): 2764-2774*.*

- RD-83088701-0

**Notes:**
Partial funding of this work obtained from US EPA.

## log10(Chla) ~ latitude + longitude

**0.798**

** **

**0.432**

** **

**0.925**

** **

**1.420**

** **

**0.998**

** **

**0.796**

** **

**1.240**

** **

**1.600**

** **

**1.210**

** **

**0.900**

** **

**0.870**

** **

**1.480**

** **

**1.140**

** **

**0.734**

** **

**MSE=0.1092**
**14 leaves**

- Preliminary Findings

- log10(Chl*a*) ~ latitude + longitude

**Notes:**
Figure 1 – Partitions of a CART model of the form: log10(Chlorophyll a) ~ latitude + longitude.  Partitions represent regions of homogeneous mean log10 Chlorophyll a concentrations, represented by the number in each partition.

Here we se a conventional CART model fit to lake-wide, summer average log10 Chla using only latitude and longitude to subset the data.
  Numbers in each partition are mean log10 Chla in each
 subset, or region.

MSE and complexity of model.

If we loosen the spatial restriction on predictors by including all of them, MSE is improved.

## log10(Chla)~TP+TN+Lat+Lon+Qin+Z+w

- |

**TP < 0.0395**

** **

**TP < 0.0145**

** **

**Lon <**

**-**

**110.88**

** **

**TN < 0.832**

** **

**TN <  2.175**

** **

**TN <  0.815**

** **

**Z < 7.997**

** **

**Z <  0.0495**

** **

**Qin < 4.018**

** **

**0.2547**

** **

**0.5869**

** **

**0.7827**

** **

**1.0290**

** **

**0.9766**

** **

**1.2520**

** **

**0.9045**

** **

**1.8810**

** **

**1.3360**

** **

**1.2480**

** **

**MSE=0.0802**
**10 leaves**

- Preliminary Findings

- log10(Chla)~TP+TN+Lat+Lon+Qin+Z+w

**Notes:**
 Figure 2  – Dendrogram of CART model: log10(Chla)~TP+TN+Lat+Lon+Qin+Z+w.  Because 7 predictors are used, a two-dimensional plot of resulting partitions such as in Figure 1 is not possible.  Numbers below each of the 10 terminal nodes (leaves) of this tree model are the predicted log10 chlorophyll a concentration.  

CART Model using these predictors, not just lat and lon. 
 Note the improved MSE.  Still, numbers at bottom of 
branches are mean log10 Chla for subsets defined by these 
rules.

If we loosen the spatial restriction on predictors by including all of them, MSE is improved.

## Results

- Bayesian Treed Model

- Preliminary Findings

## BTREED Model – partitions (T)

**MSE=0.074**
**4 leaves**

- Preliminary Findings

- BTREED Model – partitions (*T*)

**Notes:**
Figure 3 – Partitions on Latitude and Longitude of the Bayesian Treed model.  This prior resulted in a four terminal node model.  A larger prior mean may result in a model with more terminal nodes.

Note improvement in MSE, with fewer partitions, because
 some variation in log10 Chla is now explained by the full
 regression model in the end nodes, not just means as
 before.

Now this is where visualization of results becomes a bit problematic.

## Preliminary Findings

**Notes:**
X’s are standardized before model fitting, so magnitude of ’s is proportinal to “importance”.

Table entries are:
	Bayesian parameter estimate
	OLS parameter estimate
	OLS parameter SE
And, bold OLS entries are significantly different from zero at  = 0.05

TP is 2x more important in NE than anywhere else.

TN n.s. (actually negative ) in NE
	and most important in the SE!!!!
Z in SE is unimportant.

## N E

- partres

- -0.6

- -0.2

- 0.2

- 0.4

- 0.6

- -0.6

- -0.2

- 0.0

- 0.2

- 0.4

- 0.6

- n=164

- slope= 0.356836

- slope= 0.343962

- slope= 0.799635

- slope= -0.158648

- slope= 0.331859

- slope= 0.174524

- slope= -0.404363

- slope= -0.085098

- slope= 0.352297

- slope= 0.107354

- slope= 0.378674

- slope= 0.331706

- slope= 0.733359

- slope= 0.4172

- slope= 0.022849

- slope= -0.275202

- slope= 0.031206

- slope= -0.086991

- slope= -0.322522

- slope= 0.016941

- Preliminary Findings

**Notes:**
Figure 4 – Partial residuals plots for the four terminal node regression models of the Bayesian Treed Model.   Solid lines are the (shrunken) Bayesian regression coefficients, while dotted lines are the OLS estimates of these coefficients.

The models in each region differ.

TP is 2x more important in NE than anywhere else.

TN n.s. (actually negative ) in NE
	and most important in the SE!!!!
Z in SE is unimportant.

## Next Steps

## Next Steps

- More predictor variables
- Apply these methods to the Nutrient Criteria Database
- Use resultant tree structures to identify important hierarchical structure
- Explore these structures with other Hierarchical Bayesian methods
- Non-linear specification? Spline basis functions in leaf model or inclusion of all predictors in tree
- Tools

**Notes:**
1. Watershed char., inflow characteristics, waterbody morphometry, etc.
2. Training and test data. To evaluate out-of-sample pred. uncertainty.
3. Also explore geographical surrogates (ecoregions, etc.)
4. WinBUGS for DAGs
5. As in intermediate tree of WR paper.
Visualization (dendrograms and partial residuals plots) and 
Implementation (in Excel?).

Prediction: BMA

## Thanks!

- EPA STAR program for funding.
- Hugh Chipman, Univ. of Waterloo and
- 	Robert McColloch, University of Chicago for BCART/BTREED computer code.

- RD-83088701-0