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 up to the model presented by Clark et al. 2007 Ecological Applications that combines 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

```{r}
library(rjags)
library(ecoforecastR)
```

Steps 1-4 have already been done for you and leverages functions that are part of the [PEcAn system](https://pecanproject.org). Specifically, they are within PEcAn's land data R package, which can be downloaded and installed off Github using devtools

```
if(!require(PEcAn.data.land)){
  library(devtools)
  install.packages(c("digest","dplR","PeriodicTable"),repos = "https://cloud.r-project.org")
  devtools::install_github("PecanProject/pecan/base/logger")
  devtools::install_github("PecanProject/pecan/base/remote")
  devtools::install_github("PecanProject/pecan/base/utils")
  devtools::install_github("PecanProject/pecan/base/db")
  devtools::install_github("PecanProject/pecan/modules/data.land")
  library(PEcAn.data.land)
}

## 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)

```
Instead we'll just load and investigate the final, prepared data object

```{r}
load("data/Activity07.RData")
# 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, x[i-1], and the mean growth rate, mu
  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, fig.asp=0.8}

n.iter = 500                           ## INCREASE THIS NUMBER FOR ACTUAL ANALYSES **************************************

## 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 (subtract observed diameter increments off from the observed diameter)
  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}
#### Helper function
##' @param w mcmc object containing matrix outputs
##' @param pre prefix (variable name) for the matrix variable to be extracted
##' @param numeric boolean, whether to coerce class to numeric
parse.MatrixNames <- function(w, pre = "x", numeric = FALSE) {
  w <- sub(pre, "", w)
  w <- sub("[", "", w, fixed = TRUE)
  w <- sub("]", "", w, fixed = TRUE)
  w <- matrix(unlist(strsplit(w, ",")), nrow = length(w), byrow = TRUE)
  if (numeric) {
    class(w) <- "numeric"
  }
  colnames(w) <- c("row", "col")
  return(as.data.frame(w))
} # parse.MatrixNames

  #### 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
  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 track down additional covariate data or 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 **

## Applying this to your own forecast

* If you already have a state-space model running for your own model (Exercise 06) you can add additional data constraints by adding additional likelihoods to the JAGS code. For simple cases this is relatively straightforward, but you'll want to think very carefully about many of the issues raised in Chapters 6 and 9 when constructing your likelihood (model and data biases, autocorrelation, spatial and temporal aggregation, proxy data, non-Gaussian error, etc.)

* If you are forecasting multiple variables (e.g. different sites, individuals, species, etc.), the code above provides a useful template for how to organize loops over multiple timeseries.

* The code above also provides a useful template for adding different random effects 

* While the `fit_dlm` function discussed in Exercise 06 doesn't yet support data fusion models or multivariate models, the `ParseFixed` function can be used to help expand a basic random walk model to a more complex model.

* When developing your own more sophisticated analyses, make sure to start simple and add complexity incrementally. By only adding one new thing at a time it becomes much easier to debug when something goes wrong and to stop when you push the model beyond the constraints the data are able to provide.