Stanimirova_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=20)

## 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) #prior on the uncertainty in the dbh - prior on precision 
  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.height=8}

n.iter = 20000

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

# #check for MCMC convergence
# gelman.diag(jags.out)
# GBR <- gelman.plot(jags.out)
# 
# ## determine the first iteration after convergence
# burnin <- GBR$last.iter[tail(which(GBR$shrink[,1,2] > 1.1),1)+1]
# ## check for no burn-in case
# if(length(burnin) == 0) burnin = 1
# ## remove burn-in
# jags.burn <- window(jags.out,start=burnin)
# ## check diagnostics post burn-in
# gelman.diag(jags.burn)
# plot(jags.burn)
  
## we need more samples because the MCMC hasn't converged
  ## 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 now 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.


```{r, fig.height=8}
TreeDataFusionMVRE = "
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)
  
  }"
  
## compile JAGS model
j.model.RE   <- jags.model (file = textConnection(TreeDataFusionMVRE),
                          data = data,
                          inits = init,
                          n.chains = 3)

jags.out.RE   <- coda.samples (model = j.model.RE,
                              variable.names = c("tau_add","tau_dbh","tau_inc","mu","tau_ind","tau_yr"),
                              n.iter = min(n.iter,2000))

plot(jags.out.RE)

## run MCMC
jags.out.RE   <- coda.samples (model = j.model.RE,
                            variable.names = c("x","tau_add","tau_dbh","tau_inc","mu",
                                                "tau_ind","tau_yr","ind","year"),
                            n.iter = n.iter)

```


Diagnostic plots for the posterior CI for growth and DBH and compare these to observations.

Randomly selected samples have relatively tight confidence intervals compared with number 49. The CI for 49 is much wider and some of the points are even outside the CI.
The CI of number 49 didn't improve even thou we added random effects and explained more of the varibaility. The CI intervals for DBH decreased overall and the CI intervals for ring increment remained similar.  


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


Histogram of fixed effect, mu, and the standard deviations to make them easier to interpret


In the NRE (no random effects) model there is a weak negative correlation between tau_add and tau_inc. The rest of the parameters are not correlated. We want to know if the parameters are correlated with each other in order to know if they are separable and if they provide independent information about tree growth. If the two parameters are correlated with each other we cannot pull them apart and we should use their joint probability. If the parameters are independent of each other then we can use their marginal probability. The correlation between variables decreased significantly when we included random effects. 

When comparing the tau histograms between the two models, we can see that the variance of tau_add and the variance in tau_inc have decreased. This makes sense because we have added a random effect for year and a random effect for increment thereby accounting for some of the variance we previously saw in the tau_add and tau_inc. We have successfully partioned the variance further. tau_dbh did not decrease after we added the random effects, in fact it stayed almost exactly the same. This is because the random effects did not help us reduce the dbh variance, which is largely due to sampling error. Based on the mu histogram, it appears that we had constrained mu better with the model when we did not include the random effects (pink). We were able to chip away at the uncertainty and attribute more and more of the variability to other components (for now to random effects but eventually to fixed effects). Although we have tighter variance, we are less certain about where mu is exactly because mu gets more uncertain as we add more parameters. 


```{r}
 
## process model
  vars.RE = (1:ncol(out.RE))[-c(which(substr(colnames(out.RE),1,1)=="x"),grep("tau",colnames(out.RE)),
                          grep("year",colnames(out.RE)),grep("ind",colnames(out.RE)))]


## plot histogram for the random effects model -- mu

par(mfrow=c(1,1))
  for(i in vars.RE){
    hist(out.RE[,i],main=colnames(out.RE)[i], col="blue")
  }
  if(length(vars.RE)>1) pairs(out.RE[,vars.RE])

## plot histogram for the no random effects model -- mu 

  par(mfrow=c(1,1))
  for(i in vars){
    hist(out[,i],main=colnames(out)[i], add=TRUE, col="pink")
    legend("topright",c("RE", "NRE"), fill=c("blue", "pink"))
  }
  if(length(vars)>1) pairs(out[,vars])


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

## histograms of the taus with NRE and RE models 

par(mfrow=c(1,3))
hist(1/sqrt(out[,"tau_add"]),main=colnames(out)[2], col="blue", xlim=c(0.05,0.15))
hist(1/sqrt(out.RE[,"tau_add"]),main=colnames(out.RE)[2], col="yellow", add=TRUE)
hist(1/sqrt(out[,"tau_dbh"]),main=colnames(out)[3], col="blue", xlim=c(0.3,0.8))
hist(1/sqrt(out.RE[,"tau_dbh"]),main=colnames(out.RE)[3], col="yellow", add=TRUE)
hist(1/sqrt(out[,"tau_inc"]),main=colnames(out)[4], col="blue", xlim=c(0.05,0.2))
hist(1/sqrt(out.RE[,"tau_inc"]),main=colnames(out.RE)[4], col="yellow", add=TRUE)
legend("topright",c("NRE", "RE"), fill=c("blue", "yellow"))

cor(prec.RE) ## correlation between taus 
pairs(prec.RE)


```

For the individual effects: first ordering the effects by plot and the second ordering them by species.

Based on the diagnostic plots there doesn't seem to be a big plot effect. Plot effect would be a random effect. I could make sure that there is no plot effect by calculating the plot means to see if they are significantly different. If I determine that they are significantly different then I  can include plot effect and see if it explains any of the residual variability. However, there is some possible species effects. For example, red maples appear to be growing slower than the average. So it's probably worthwhile to look at species effects. Species would be a fixed effect in terms of how fast red maples are growing. At the biome level the species effect could be a random one.   


```{r}
  par(mfrow=c(1,1))
  ### YEAR
  year.cols.RE = grep("year",colnames(out.RE))
  if(length(year.cols.RE>0)){
    ci.yr.RE <- apply(out.RE[,year.cols.RE],2,quantile,c(0.025,0.5,0.975))
    plot(data$time,ci.yr.RE[2,],type='n',ylim=range(ci.yr.RE,na.rm=TRUE),main="Year Effect",ylab="cm")
    ciEnvelope(data$time,ci.yr.RE[1,],ci.yr.RE[3,],col="lightBlue")
    lines(data$time,ci.yr.RE[2,],lty=1,lwd=2)
    abline(h=0,lty=2)
  }
  
  ### INDIV
  ind.cols.RE= which(substr(colnames(out.RE),1,3)=="ind")
  if(length(ind.cols.RE)>0){
    boxplot(out.RE[,ind.cols.RE],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.RE[,ind.cols.RE],2,mean),combined$PLOT,mean)
    table(combined$PLOT)
    
    spp = combined$SPP
    boxplot(out.RE[order(spp),ind.cols.RE],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.RE[,ind.cols.RE],2,mean),combined$SPP,mean)
  }
```


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. 

I propose a fixed effect - precipitation - which can further the decompose the variability of the year random effect and chip away at the unexplained variability. That way we can account for variability and try to explain its source. A fixed effect describes the slope of the relationship. Precipitation is a fixed effect because the relationship between precipitation and growth is fixed although the amount of precipitation might vary from experiment to experiment. 


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

I would plot the precipitation on the x axis and the random year effect on the y axis. This will tell me how much of year-to-year variability is explained by precipitation. If precipitation explains all of the year to year variability than the random year effect is no longer useful. However, if precipitation explains only half of the year-to-year variability then I will keep the random effect in order to account for the unexplained variability still present. 

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

Data models describe alternative methods of measuring the same thing so in this case it is dbh and tree rings (inc), which can both be used to measure tree growth. The process model includes both random and fixed effects. 

In order to determine the prior for the relationship between precipitation and growth, I would conduct expert elicitation to determine the prior on the slope of the relationship. I would select the functional form based on the natural range of values, and based on a literature review. Since precipitation is a continuous variable I will choose a continuous probability distribution.

beta is the fixed effect for the slope of the relationship between rainfall and tree growth. The slope describes how much growth changes for each change in precipitation. The prior on beta is a normal distribution. This is because the slope between the two can be positive when precipitation is high and negative when the precipitation is low. Therefore, a normal distiribution fits the natural range of values. 

raint[t] - observations of rainfall over the time period 


```
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] + beta*rain[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)
  tau_beta ~ dnorm(mean,sd)
  mu       ~ dnorm(0.5,0.5)
  
  }
```

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