Exercise_07_TreeRings.Rmd

Fusing Times-Series Data: Tree Rings and Forest Inventory
========================================================

In this exercise we will extend the state-space framework to combine multiple data streams with different observation errors and to separate observation error from process error. We will also demonstrate how to add hierarchical random effects to partition the process error into multiple sources.

Specifically, we will be building upto the model presented by Clark et al. 2007 Ecological Applications in order to combine tree ring data with forest inventory data. Unlike the original model, which was written all in R, we will rewrite this model into JAGS, which makes it easier to see what is going on and to modify the model. In this exercise we will utilize data from a collection of small plots at the Harvard Forest, Petersham, MA.

We will divide this analysis into a number of steps, which we will encapsulate into functions to make them easier to understand and run. Thus we will begin by defining these functions. Specifically, the steps will be:

1. load forest inventory data
2. load tree ring data
3. match the tree core and inventory data for individual trees and merge these data sets into one data frame
4. format this data into a list for input into JAGS
5. run the JAGS model
6. visualize the output

For steps 2-4 we will leverage functions already written to deal with these steps that are part of the PEcAn system we will explore next week. Specifically, they are within PEcAn's land data R package, which can be downloaded and installed off Github using devtools


```{r}
if(!require(PEcAn.data.land)){
  install.packages(c("digest","dplR"))
  devtools::install_github("PecanProject/pecan/utils")
  devtools::install_github("PecanProject/pecan/db")
  devtools::install_github("PecanProject/pecan/modules/data.land")
  require(PEcAn.data.land)
}
require(rjags)
```

```{r, echo=FALSE}
#plots a confidence interval around an x-y plot (e.g. a timeseries)
ciEnvelope <- function(x,ylo,yhi,...){
  polygon(cbind(c(x, rev(x), x[1]), c(ylo, rev(yhi),
                                      ylo[1])), border = NA,...) 
}
```


```{r}
## 1. Read tree data
trees <- read.csv("data/H2012AdultFieldData.csv")

## 2. Read tree ring data
rings <- Read_Tucson("data/TUCSON/")

## 3. merge inventory and tree ring data, extract most recent nyears
combined <- matchInventoryRings(trees,rings,nyears=15)

## take a look at the first few rows of data to see the structure
knitr::kable(combined[1:5,])

## 4. organize data into a list
data <- buildJAGSdata_InventoryRings(combined)

# y = increment (tree x year)
# z = dbh (tree x year)
# make sure to take a look at all the priors!
str(data)

```

Now that we have the data prepped we need to fit the model itself. The bulk of this code is just the same JAGS syntax we've used before, so lets focus on the JAGS code itself. To begin with, lets look back at the JAGS code for the random walk

```
model{
  
  #### Data Model
  for(i in 1:n){
    y[i] ~ dnorm(x[i],tau_obs)
  }
  
  #### Process Model
  for(i in 2:n){
    x[i]~dnorm(x[i-1],tau_add)
  }
  
  #### Priors
  x[1] ~ dnorm(x_ic,tau_ic)
  tau_obs ~ dgamma(a_obs,r_obs)
  tau_add ~ dgamma(a_add,r_add)
}
```

Since we're fusing two data sources, we'll need to add a second data model. We'll also modify our process model to include a mean growth rate term. Finally, we'll need to specify priors on both observation errors, the process error, and the mean.

```
model{

  #### Data Model: DBH
  for(i in 1:n){
    z[i] ~ dnorm(x[i],tau_dbh)
  }

  #### Data Model: growth
  for(i in 2:n){
    inc[i] <- x[i]-x[i-1]
    y[i] ~ dnorm(inc[i],tau_inc)
  }

  #### Process Model
  #### Dnew is the expected new diameter given the previous diameter and the mean growth rate
  for(i in 2:n){
    Dnew[i] <- x[i-1] + mu  
    x[i]~dnorm(Dnew[i],tau_add)
  }

  #### Priors
  x[1] ~ dnorm(x_ic,tau_ic)     ## initial DBH
  tau_dbh ~ dgamma(a_dbh,r_dbh) ## observation error: DBH
  tau_inc ~ dgamma(a_inc,r_inc) ## observation error: tree rings
  tau_add ~ dgamma(a_add,r_add) ## process error: growth
  mu ~ dnorm(0.5,0.5)           ## mean growth
}
```

This code would work perfectly if we only had only measured a single tree, but we measured a number of trees so next need to modify the code to work with tree-by-year matrices of DBH and growth.

```
model{

  ### Loop over all individuals
  for(i in 1:ni){
  
  #### Data Model: DBH
  for(t in 1:nt){
  z[i,t] ~ dnorm(x[i,t],tau_dbh)
  }
  
  #### Data Model: growth
  for(t in 2:nt){
  inc[i,t] <- x[i,t]-x[i,t-1]
  y[i,t] ~ dnorm(inc[i,t],tau_inc)
  }
  
  #### Process Model
  for(t in 2:nt){
  Dnew[i,t] <- x[i,t-1] + mu
  x[i,t]~dnorm(Dnew[i,t],tau_add)
  }
  
  x[i,1] ~ dnorm(x_ic,tau_ic)
  }  ## end loop over individuals
  
  #### Priors
  tau_dbh ~ dgamma(a_dbh,r_dbh)
  tau_inc ~ dgamma(a_inc,r_inc)
  tau_add ~ dgamma(a_add,r_add)
  mu ~ dnorm(0.5,0.5)
}
```

Finally, since growth is indexed by both tree and year, lets add random effects for both individuals and years. In this case our process model now becomes

Dnew[i,t] <- x[i,t-1] + mu + ind[i] + year[t]

where ind and year are the random effects for individual and year respectively. Next, we'll need to specify the distributions that these random effects are drawn from, as well as the priors on the random effect variances

```
model{
  
  ### Loop over all individuals
  for(i in 1:ni){
  
  #### Data Model: DBH
  for(t in 1:nt){
    z[i,t] ~ dnorm(x[i,t],tau_dbh)
  }
  
  #### Data Model: growth
  for(t in 2:nt){
    inc[i,t] <- x[i,t]-x[i,t-1]
    y[i,t] ~ dnorm(inc[i,t],tau_inc)
  }
  
  #### Process Model
  for(t in 2:nt){
    Dnew[i,t] <- x[i,t-1] + mu + ind[i] + year[t]
    x[i,t]~dnorm(Dnew[i,t],tau_add)
  }
  
  ## individual effects
  ind[i] ~ dnorm(0,tau_ind)
  
  ## initial condition
  x[i,1] ~ dnorm(x_ic,tau_ic)
  }  ## end loop over individuals
  
  ## year effects
  for(t in 1:nt){
    year[t] ~ dnorm(0,tau_yr)
  }
  
  
  #### Priors
  tau_dbh ~ dgamma(a_dbh,r_dbh)
  tau_inc ~ dgamma(a_inc,r_inc)
  tau_add ~ dgamma(a_add,r_add)
  tau_ind ~ dgamma(1,0.1)
  tau_yr  ~ dgamma(1,0.1)
  mu ~ dnorm(0.5,0.5)
  
  }
```

Putting this all together gives the following R code for the base case (no random effects) 

```{r}

n.iter = 500

## this code fuses forest inventory data with tree growth data (tree ring or dendrometer band)
## for the same plots. Code is a rewrite of Clark et al 2007 Ecol Appl into JAGS
TreeDataFusionMV = "
model{

  ### Loop over all individuals
  for(i in 1:ni){
  
  #### Data Model: DBH
  for(t in 1:nt){
  z[i,t] ~ dnorm(x[i,t],tau_dbh)
  }
  
  #### Data Model: growth
  for(t in 2:nt){
  inc[i,t] <- x[i,t]-x[i,t-1]
  y[i,t] ~ dnorm(inc[i,t],tau_inc)
  }
  
  #### Process Model
  for(t in 2:nt){
  Dnew[i,t] <- x[i,t-1] + mu
  x[i,t]~dnorm(Dnew[i,t],tau_add)
  }
  
  x[i,1] ~ dnorm(x_ic,tau_ic)
  }  ## end loop over individuals
  
  #### Priors
  tau_dbh ~ dgamma(a_dbh,r_dbh)
  tau_inc ~ dgamma(a_inc,r_inc)
  tau_add ~ dgamma(a_add,r_add)
  mu ~ dnorm(0.5,0.5)
}"

  ## state variable initial condition
  z0 = t(apply(data$y,1,function(y){-rev(cumsum(rev(y)))})) + data$z[,ncol(data$z)] 
  
  ## JAGS initial conditions
  nchain = 3
  init <- list()
  for(i in 1:nchain){
    y.samp = sample(data$y,length(data$y),replace=TRUE)
    init[[i]] <- list(x = z0,tau_add=runif(1,1,5)/var(diff(y.samp),na.rm=TRUE),
                      tau_dbh=1,tau_inc=500,tau_ind=50,tau_yr=100,ind=rep(0,data$ni),year=rep(0,data$nt))
  }
  
  ## compile JAGS model
  j.model   <- jags.model (file = textConnection(TreeDataFusionMV),
                           data = data,
                           inits = init,
                           n.chains = 3)
  ## burn-in
  jags.out   <- coda.samples (model = j.model,
                              variable.names = c("tau_add","tau_dbh","tau_inc","mu","tau_ind","tau_yr"),
                              n.iter = min(n.iter,2000))
  plot(jags.out)
  
  ## run MCMC
  jags.out   <- coda.samples (model = j.model,
                              variable.names = c("x","tau_add","tau_dbh","tau_inc","mu",
                                                 "tau_ind","tau_yr","ind","year"),
                              n.iter = n.iter)
```

Next, lets generate some diagnostic plots to look at the model. First, lets plot the posterior CI for growth and DBH and compare these to observations. Since we have scores of cores and trees, we'll pick a random subset of trees to check. One thing that's critical to note is that for the confidence intervals on growth that these are calculated pathwise -- we're looking at the growth from a whole MCMC iteration -- rather than pairwise (i.e. subtracting the posterior distribution for DBH at one point from the posterior distribution of DBH at the next). Because there's high correlations between successive time points, the pathwise uncertainty estimates are considerably lower in uncertainty -- essentially saying that we know can know the growth rate of the tree better than we can know the actual size of the tree

```{r, fig.height=8}
  #### Diagnostic plots
  
  ### DBH
  layout(matrix(1:8,4,2))
  out <- as.matrix(jags.out)
  x.cols = which(substr(colnames(out),1,1)=="x")   ## which columns are the state variable, x
  ci <- apply(out[,x.cols],2,quantile,c(0.025,0.5,0.975))
  ci.names = parse.MatrixNames(colnames(ci),numeric=TRUE)
  
  smp = c(sample.int(data$ni,3),49)  ## I've rigged the sampling to make sure you see tree 49!
  for(i in smp){
    sel = which(ci.names$row == i)
    plot(data$time,ci[2,sel],type='n',ylim=range(ci[,sel],na.rm=TRUE),ylab="DBH (cm)",main=i)
    ciEnvelope(data$time,ci[1,sel],ci[3,sel],col="lightBlue")
    points(data$time,data$z[i,],pch="+",cex=1.5)
  }
  
  ## growth
  for(i in smp){
    sel = which(ci.names$row == i)
    inc.mcmc = apply(out[,x.cols[sel]],1,diff)
    inc.ci = apply(inc.mcmc,1,quantile,c(0.025,0.5,0.975))*5
    
    plot(data$time[-1],inc.ci[2,],type='n',ylim=range(inc.ci,na.rm=TRUE),ylab="Ring Increment (mm)")
    ciEnvelope(data$time[-1],inc.ci[1,],inc.ci[3,],col="lightBlue")
    points(data$time,data$y[i,]*5,pch="+",cex=1.5,type='b',lty=2)
  }
```


Second, let's look at the histogram of our fixed effect, mu, and the precisions. Let's also convert the precisions to standard deviations to make them easier to interpret

```{r}
  ## process model
  vars = (1:ncol(out))[-c(which(substr(colnames(out),1,1)=="x"),grep("tau",colnames(out)),
                          grep("year",colnames(out)),grep("ind",colnames(out)))]
  par(mfrow=c(1,1))
  for(i in vars){
    hist(out[,i],main=colnames(out)[i])
  }
  if(length(vars)>1) pairs(out[,vars])

  ## Standard Deviations
  #layout(matrix(c(1,2,3,3),2,2,byrow=TRUE))
  par(mfrow=c(2,3))
  prec = out[,grep("tau",colnames(out))]
  for(i in 1:ncol(prec)){
    hist(1/sqrt(prec[,i]),main=colnames(prec)[i])
  }
  cor(prec)
  pairs(prec)
```


Third, let's look at the random effects. It is easy enough to plot the year effects by year. For the individual effects we'll plot these twice, first ordering the effects by plot and the second ordering them by species.

```{r}
  par(mfrow=c(1,1))
  ### YEAR
  year.cols = grep("year",colnames(out))
  if(length(year.cols>0)){
    ci.yr <- apply(out[,year.cols],2,quantile,c(0.025,0.5,0.975))
    plot(data$time,ci.yr[2,],type='n',ylim=range(ci.yr,na.rm=TRUE),main="Year Effect",ylab="cm")
    ciEnvelope(data$time,ci.yr[1,],ci.yr[3,],col="lightBlue")
    lines(data$time,ci.yr[2,],lty=1,lwd=2)
    abline(h=0,lty=2)
  }
  
  ### INDIV
  ind.cols= which(substr(colnames(out),1,3)=="ind")
  if(length(ind.cols)>0){
    boxplot(out[,ind.cols],horizontal=TRUE,outline=FALSE,col=combined$PLOT,main="Individual Effects By Plot",xlab="cm")
    abline(v=0,lty=2)
    ## calculate plot-level means for random effects
    tapply(apply(out[,ind.cols],2,mean),combined$PLOT,mean)
    table(combined$PLOT)
    
    spp = combined$SPP
    boxplot(out[order(spp),ind.cols],horizontal=TRUE,outline=FALSE,col=spp[order(spp)],main="Individual Effects By Species",xlab="cm")
    abline(v=0,lty=2)
    spp.code = levels(spp)[table(spp)>0]
    legend("bottomright",legend=rev(spp.code),col=rev(which(table(spp)>0)),lwd=4)
    ## calculate species-level means for random effects
    tapply(apply(out[,ind.cols],2,mean),combined$SPP,mean)
  }
```



By default this code is set to run with a small number of years (15), and a much too low number of MCMC iterations (500), just so that the code with "knit" quickly initially. For your analyses you should obviously increase these -- I found that convergence was adequate with around 20,000 samples, though I probably would run 10x longer than that for a publishable analysis. However, such an analysis would take hours to run.

Assignment:

1. Run the model initially with random effects off

2. Rerun the model with random effects on. Compare this to the previous run. What is the relative partitioning of uncertainties in the different versions of the model among observation error, process error, and the different random effects? What does the size of these effects suggest about the drivers of uncertainty in tree growth?

3. Based on the diagnostics, propose an additional effect (fixed or random) to add to the model. Such an effect should plausibly chip away at a sizable fraction of the unexplained variability -- you wouldn't want to propose an effect that isn't associated with systematic variability. 

4. Explain any additional exploratory analyses you would perform (e.g. plotting your proposed covariate against one of the random effects). 

5. Write the JAGS code that would fit the proposed model (note: you don't have to run this model, just propose the code)

** BECAUSE THE PRODUCTION VERSION OF THIS CODE TAKES A LONG TIME TO RUN, PLEASE SUBMIT THE KNIT HTML NOT THE Rmd **