Ecological Archives E096-096-A5

Chris H. Wilson, T. Trevor Caughlin, David J. Civitello, and S. Luke Flory. 2015. Combining mesocosm and filed experiments to predict invasive plant performance: a hierarchical Bayesian approach. Ecology 96:10841092. http://dx.doi.org/10.1890/14-0797.1

Appendix E. Details on model development, particularly working with the gamma and negative binomial parameterizations in JAGS, plus pseudo-code for generating cross-validated posterior predictions.

AppE_text

Pseudo-code for running cross-validation checks:

        

sink("ranmodelgroE.txt") #the sink function drops the code as a text file into your working directory
cat('
    
# Repeat model specification from above 

    }
    
    ',fill=TRUE)
sink()

# data block: define your data inputs here

# The following code sets observations in the dataset to NA
# in our case, this is based on the held-out site 

um<-unique(massdata$sites)
loocv<-vector("list",length(um))

for(i in 1:length(um)) {
  
  mdat<-massdata
  mdat$totmass[which(mdat$Site==um[i])]<-NA
  loocv[[i]]<-mdat
}


# This code block creates a list-vector to save the results of the  #multiple model fits in
# and then  executes a loop to fit the model N times, where N= number #of hold-out sites, or
# whatever

results<-vector("list",length(um))
set.seed(101)
for(i in 1:length(um)) {
  mdat<-loocv[[i]]
  mass_data = list(mass = c((mdat$totmass), rep(NA, times = 100)), light = c(light.dat$massdata.Light, seq(1,100,1)), 
                   site = c(mdat$sites, seq(22,121,1)), seeds = seed_data$seeds[1:24], 
                   biomass = seed_data$biomass[1:24])
  params = list("a", "b", "alpha2", "alpha1", "alpha", "p","rr","lambda","mass", "mu")
  ranM = jags(mass_data ,inits=NULL,unlist(params),model.file="ranmodelgroE.txt",n.chains=5,n.iter=12000,n.burnin=2000,n.thin= 50)
  
  results[[i]]<-gipar(ranM)
  
}
save(results,file="results.Rdata")
 
# The following is a function for extracting out only the predictions for the held-out sites
 MINit<-rep(NA,times=21)
MAXit<-rep(NA,times=21)

for(i in 1:21) {
  MINit[i]<-min(which(massdata$sites==um[i]))
  MAXit[i]<-max(which(massdata$sites==um[i]))
}


nameo<-function(name){
  listo<-vector("list",21)
  for(i in 1:21) {
    
    give<-results[[i]][which(names(results[[i]])==name)]
    listo[[i]]<-give[[1]][,c(MINit[i]:MAXit[i])]
    
  }
  
  return(do.call("cbind",listo))
}

      
        

[Back to E096-096]