Case Study 1: Distribution and Abundance of Tree Species Along Climate Gradients

1. Case Study 1:Distribution and Abundance of Tree Species Along Climate Gradients (Sometimes the PDF is more interesting than the actual model…) Canham, C. D. and R. Q. Thomas. 2010. Frequency, not relative abundance, of temperate tree species varies along climate gradients in eastern North America. Ecology 91:3433-3440

2. The Data • Relative abundance of the 24 most common tree species in the northeastern US • In ~ 20,000 US Forest Service Forest Inventory and Analysis (FIA) plots from 19 northeastern states (Maine to Wisconsin, south to Kentucky and Virginia) • Climate (mean annual temperature and average annual precipitation) for each plot, averaged over the period since the previous census • Using 800-m resolution gridded PRISM climate data, with bi-linear interpolation to true plot locations

3. The Basic Question • How does the abundance of tree species vary along climate gradients? http://www.fs.fed.us/ne/delaware/atlas/s318.html

4. Basic Approach • Develop regression models that predict abundance of a given tree species in a plot as a function of climate at the location of the plot… • Initial decisions: • What to use as a measure of abundance? • I chose relative abundance over absolute abundance • What sorts of functions could describe variation in relative abundance along climate gradients • Compare a Gaussian function with a null model that was flat

5. But is a flat line an appropriate null hypothesis? • What about range limits? • And shouldn’t even the Gaussian model have truncated tails?

6. So, what do the data look like? Acer rubrum (red maple) Quercus prinus ( chestnut oak) It doesn’t look like there’s much hope for nice Gaussian niche curves…

7. A first regression model (in R) # Square Gaussian - predicts mean of 0 when below lo or above hi square.gauss.model <- function(a,m,b,lo,hi,X) { ifelse(X < lo,0,ifelse(X > hi,0,a*exp(-0.5*(((X-m)/b)^2)))) } • 5 Parameters: • “a” determines height of curve at mode • “m” temperature at peak of curve • “b” breadth of the curve • “lo” lower temperature limit • “hi” upper temperature limit Note that a separate “null” model is not really necessary, because if “b” is large, the curve is flat

8. What might be an appropriate PDF? The data are clearly not normally distributed… In fact, there are a whole lot of zeros….

9. Divide the analysis up into two parts • Predict probability of “presence” (i.e. non-zero abundance) – • Separately predict relative abundance when present These two components of abundance have very different ecological meanings….

10. A PDF for relative abundance when present:The Gamma Function One likely choice for plots where a species is present… Gamma PDF with scale parameter = 22

11. First ResultsRelative Abundance (when present) Bottom line: for all 24 species, relative abundance is highly variable, but shows very little trend in the mean across the entire range of a species’ climate niche limits Note: these results are for relative abundance, given that a species is present…

12. But what about all of those zeros?Creating a PDF that can model both “presence” and “relative abundance” • “Zero-Inflated Distributions” # ZERO INFLATED NORMAL PDF zinf_norm_PDF <- function(x,mean,sd,pz) { log(ifelse(x==0,pz + (1-pz)*dnorm(0,mean,sd,log=F), (1-pz)*dnorm(x,mean,sd,log=F))) } # ZERO INFLATED GAMMA PDF zinf_gamma_PDF_climate <- function(x,mean,scale,px) { shape <- mean/scale loglh <- log(ifelse(x==0,pz + (1-pz)*dgamma(0,shape=shape,scale=scale,log=F), (1-pz)*dgamma(x,shape=shape,scale=scale,log=F))) return(loglh) }

13. But why stop there? Is presence constant across climate gradients, or does it vary? # ZERO INFLATED GAMMA PDF with pz a gaussian function of the # independent (climate) variable zinf_gamma_PDF_climate <- function(x,mean,scale,pa,pm,pb,px) { shape <- mean/scale pz <- 1 - pa*(exp(-0.5*(((px-pm)/pb)^2))) # px is temp or precip in the plot loglh <- log(ifelse(x==0,pz + (1-pz)*dgamma(0,shape=shape,scale=scale,log=F), (1-pz)*dgamma(x,shape=shape,scale=scale,log=F))) return(loglh) }

14. One final complication… • How can I factor range limits into the PDF? • If probability of presence is modeled as a Gaussian function of climate, but with truncated tails (climatic limits), what likelihood should be assigned to plots outside the estimated climatic limits? (my answer – 1 in a million…) # CLIMATE DEPENDENT ZERO INFLATED GAMMA, WITH LIMITS PDF - # ARBITRARILY SET Prob(x>0|X=0) = 0.000001 zinf_limits_gamma_PDF_climate <- function(x,mean,scale,pa,pm,pb,px) { shape <- mean/scale pz <- 1 - pa*(exp(-0.5*(((px-pm)/pb)^2))) loglh <- log(ifelse(mean == 0, ifelse(x == 0,0.999999,0.000001), ifelse(x==0,pz + (1-pz) * dgamma(0,shape=shape,scale=scale,log=F), (1-pz)*dgamma(x,shape=shape,scale=scale,log=F)))) return(loglh) }

18. Extending the approach to western US species

19. Is a symmetric gaussian model a good fit? • A simple asymmetric version of the gaussian function: