#' # Wolf Population Dynamics:  comparing 3 populations
#'
#'
#' This file contains R code (and output) reproducing all results contained in:
#' Mech and Fieberg, Growth Rates and Variances of Unexploited Wolf Populations in Dynamic Equilibria. JWM.  
#' 
#' Author:  John Fieberg, jfieberg@umn.edu
#'
#' Note:  Results will differ slightly from those in the paper due to MCMC error (saving the seed does not currently 
#' work with the glm module of JAGS - see: http://martynplummer.wordpress.com/2012/08/28/lack-of-reproducibility-in-the-glm-module/)
#' 
#'  Set working directory
  wd <- ifelse(basename(getwd())=="Scripts", 
               gsub("/Scripts", "", getwd()),
               getwd())
  opts_knit$set(root.dir=wd, cache=T)
  rm(list = ls())
  
#' Tried to get a random seed for random number generation & save it so results could be reproduced (but, this does not currently work).
 load("./data/rseed.Rdata")
 if(is.null(seed.save)==T){
   seed.save<-floor(runif(1,0,100000))
   save(seed.save, file="./data/rseed.Rdata")
  }
  set.seed(seed.save)
  
#' Load R libaries of use

#+ warning=F, message=F
  library(ggplot2)
  library(PVAClone)
  library(doBy)
  library(reshape)
  library(R2jags)
  library(R2WinBUGS)
  library(mcmcplots)
 
#' Read in the data
  wdat<-read.csv("./data/Wolfdat.csv")
  
#' Calculate yearly rates of change
  nyr<-nrow(wdat)
  wdat$r1<-c(NA, log(wdat$Denali[2:nyr]/wdat$Denali[1:(nyr-1)]))
  wdat$r2<-c(NA, log(wdat$IsleRoyale[2:nyr]/wdat$IsleRoyale[1:(nyr-1)]))
  wdat$r3<-c(NA, log(wdat$SNF[2:nyr]/wdat$SNF[1:(nyr-1)]))
  
#' Functions to calculate mean, SE(of mean), and variance of lambdas when there are missing data
  mean.narm<-function(x){mean(x, na.rm=T)}
  n.narm<-function(x){length(x[is.na(x)!=T])}
  var.narm<-function(x){var(x, na.rm=T)}
  SE.mean<-function(x){sqrt(var.narm(x)/n.narm(x))}

#' Print matrix containing:  average growth rate (log-scale), SE(average growth rate), Var(growth rate).  Assumptions:  
#' - Growth rates = log(N[t]/N[t-1]) are independent over time 
#' - The distribution is stationary (i.e., not changing over time)
  dsums<-summaryBy(r1+r2+r3~1, FUN=c(mean.narm, SE.mean, var.narm, n.narm), data=wdat)
  sum.mat<- matrix(dsums, 3,4, byrow=F)
  rownames(sum.mat)<-c("Denali", "IsleRoyale","SNF")
  colnames(sum.mat)<-c("r.bar", "SE.r.bar", "Sig2.r", "n")
  sum.mat
  
#' Create Long data set for plotting....  
  wdat2a <- melt(wdat[,1:4], id=c("YR"))
  wdat2b <- melt(wdat[,c(1,5:7)], id=c("YR"))
  names(wdat2a)<-c("yr", "Population", "Count")
  names(wdat2b)<-c("yr", "Population", "r")
  wdat2b$Population<-factor(wdat2b$Population, labels=c("Denali", "IsleRoyale","SNF"))

#' Population size over time  
#+fig.height=6, fig.width=6
  qplot(yr,Count, data=wdat2a, colour=Population, group=Population, geom=c("path", "point")) + theme_bw()+xlab("Year") + ylab("") +
    theme(axis.title = element_text(size = rel(1.5)), axis.text= element_text(size = rel(1.5)),legend.title = element_text(size = rel(1.5)),
          legend.text = element_text(size = rel(1.5)))

#' Yearly growth rates over time  
#+fig.height=6, fig.width=6
  qplot(yr, r,data=wdat2b, colour=Population, group=Population, geom=c("path", "point")) + theme_bw()+
     theme(axis.title = element_text(size = rel(1.5)), axis.text= element_text(size = rel(1.5)),legend.title = element_text(size = rel(1.5)),
          legend.text = element_text(size = rel(1.5))) + xlab("Year") + ylab(expression(lambda[t]))


#' Distribution of yearly growth rates  
#+ fig.height=6, fig.width=6
  qplot(r, data=wdat2b, colour=Population, group=Population, geom=c("density"))+ theme_bw()+
     theme(axis.title = element_text(size = rel(1.5)), axis.text= element_text(size = rel(1.5)),legend.title = element_text(size = rel(1.5)),
          legend.text = element_text(size = rel(1.5)))

#' Figure 1 in paper
par(cex.lab=1.4, cex.axis=1.4, bty="O",mar=c(5.1, 5.1, 3.1, 4.1), xpd=NA, oma=c(0,0,0,4))
  with(wdat2a, plot(yr, Count, xlab="", ylab="",type="n", pch=16, las=1)) 
  with(wdat, lines(YR, Denali, type="o", pch=16, lwd=2))
  with(wdat, lines(YR, IsleRoyale, type="o", pch=16, col="DarkGray", lwd=1))
  with(wdat, lines(YR, SNF, type="o", pch=1, lwd=2, lty=2))
   legend("right", inset=c(-.28,0),cex=1.3, bty="n", legend=c("Denali","IR", "SNF"), pch=c(16,16,1),lty=c(1,1,2),lwd=c(2), ncol=1,col=c("Black", "DarkGray","Black"))
   mtext(outer=F, line=2.7, side=1, "Year", cex=1.6)
  mtext(outer=F, line=3.5, side=2, "Number of Wolves", cex=1.6)
  
#-----------------------------------------------------------------------------------------
#' source R script file that contains all of the models
#' 
  source('./Scripts/PopulationModels.R')
  
#' # Density-Independent Models 
#' 
#' ## No observation error:
#' 
#'  $$N_{t+1}=N_t\exp(a+\epsilon_t)$$
#'    
#'  $$\epsilon_t \sim N(0, \sigma_{proc}^2)$$  
#'   

#' Fit model in JAGS to Denali, Isle Royale, and SNF populations
  N<-(wdat[-c(1:28, 56), 2])
  x<-log(N)
  T<-length(N)
  di.noerror.D <- jags(data=c("x","T"),   parameters.to.save=c("a","sigma"),
    n.iter=25000, n.burnin=5000,model.file=di.noerror, n.thin=1, progress.bar="none" )

  N<-(wdat[, 3])
  x<-log(N)
  T<-length(N)
  di.noerror.IR <- jags(data=c("x","T"),   parameters.to.save=c("a","sigma"),
    n.iter=25000, n.burnin=5000,model.file=di.noerror, n.thin=1, progress.bar="none" )

  N<-(wdat[-c(1:8, 55:56), 4])
  x<-log(N)
  T<-length(N)
  di.noerror.SNF <- jags(data=c("x","T"),   parameters.to.save=c("a","sigma"),
    n.iter=25000, n.burnin=5000,model.file=di.noerror, n.thin=1, progress.bar="none" )
  
  
  
#' ## Poisson observation error:
#' 
#'  $$N_{t+1}=N_t\exp(a+\epsilon_t)$$
#'    
#'  $$\epsilon_t \sim N(0, \sigma_{proc}^2)$$  
#'  
#'  $$O_t \sim Poisson(N_t)$$
#'  
 

#' Fit model in JAGS to Denali, Isle Royale, and SNF populations
  O<-(wdat[-c(1:28, 56), 2])
  T<-length(O)
  di.poiss.D <- jags(data=c("O","T"),   parameters.to.save=c("a","sigma"),
    n.iter=25000, n.burnin=5000,model.file=di.Poisson, n.thin=1, progress.bar="none" )

  O<-(wdat[, 3])
  T<-length(O)
  di.poiss.IR <- jags(data=c("O","T"),   parameters.to.save=c("a","sigma"),
    n.iter=25000, n.burnin=5000,model.file=di.Poisson, n.thin=1, progress.bar="none" )

  O<-(wdat[-c(1:8, 55:56), 4])
  T<-length(O)
  di.poiss.SNF <- jags(data=c("O","T"),   parameters.to.save=c("a","sigma"),
    n.iter=25000, n.burnin=5000,model.file=di.Poisson, n.thin=1, progress.bar="none" )
  
   
#' ## Log-normal observation error:
#' 
#'  $$N_{t+1}=N_t\exp(a+\epsilon_t)$$ or $$X_{t+1} = X_t + a +bN_t + \epsilon_t$$
#'    
#'  $$\epsilon_t \sim N(0, \sigma_{proc}^2)$$  
#'  
#'  $$Y_t \sim N(X_t, \sigma_{obs})$$
#'  
 

#' Fit model in JAGS to Denali, Isle Royale, and SNF populations
  O<-(wdat[-c(1:28, 56), 2])
  y<-log(O)
  T<-length(O)
  di.logN.D <- jags(data=c("y","T"),   parameters.to.save=c("a","sigma", "tau"),
    n.iter=25000, n.burnin=5000,model.file=di.LogN, n.thin=1, progress.bar="none" )

  O<-(wdat[, 3])
  y<-log(O)
  T<-length(O)
  di.logN.IR <- jags(data=c("y","T"),   parameters.to.save=c("a","sigma", "tau"),
    n.iter=25000, n.burnin=5000,model.file=di.LogN, n.thin=1, progress.bar="none" )

  O<-(wdat[-c(1:8, 55:56), 4])
  y<-log(O)
  T<-length(O)
  di.logN.SNF <- jags(data=c("y","T"),   parameters.to.save=c("a","sigma", "tau"),
    n.iter=25000, n.burnin=5000,model.file=di.LogN, n.thin=1, progress.bar="none" )
  
  
  
#' ### Plot posterior densities of mean.r and sig.r
#+ fig.height=8, fig.width=8
  layout(matrix(c(1,1,2,2,3,3,4,4,5,5,0,0,6,6,7,7,8,8), ncol=6, byrow=TRUE))
  par(bty="L", cex.main=1.8, cex.lab=1.6, cex.axis=1.4)

# No observation error: a and sig_p
  plot(density(di.noerror.D$BUGSoutput$sims.matrix[,1]),  lty=1, lwd=2,xlab="", ylab="Density", main=expression(italic(a)),  ylim=c(0,26))
   lines(density(di.noerror.IR$BUGSoutput$sims.matrix[,1]), lty=1, lwd=2,col="DarkGray")
   lines(density(di.noerror.SNF$BUGSoutput$sims.matrix[,1]),  lty=2, lwd=2)
  
  
   plot(density(di.noerror.D$BUGSoutput$sims.matrix[,3]),  lty=1, lwd=2,xlab="", ylab="Density", main=expression(sigma[proc]), xlim=c(0, 0.5),ylim=c(0,22))
   lines(density(di.noerror.IR$BUGSoutput$sims.matrix[,3]), lty=1, lwd=2,col="DarkGray")
   lines(density(di.noerror.SNF$BUGSoutput$sims.matrix[,3]), lty=2, lwd=2)
   
  plot(1:5, 1:5, type="n", xaxt="n", yaxt="n", bty="n", xlab="", ylab="") 
   legend("center", xjust=0.5,  lty=c(1,1,2), lwd=c(2,2,2), col=c("black", "DarkGray", "Black"),
          legend=c("Denali National Park", "Isle Royale","Superior National Forest" ),
          y.intersp=0.8, x.intersp=1.2,
          bty="n", cex=1.4, inset=4, xpd=TRUE) 

# Poisson observation error: a and sig_p  
   plot(density(di.poiss.D$BUGSoutput$sims.matrix[,1]),  lty=1, lwd=2,xlab="", ylab="Density", main=expression(italic(a)), ylim=c(0,26))
   lines(density(di.poiss.IR$BUGSoutput$sims.matrix[,1]), lty=1, lwd=2,col="DarkGray")
   lines(density(di.poiss.SNF$BUGSoutput$sims.matrix[,1]), lty=2, lwd=2)
  
  
   plot(density(di.poiss.D$BUGSoutput$sims.matrix[,3]),  lty=1, lwd=2,xlab="", ylab="Density", main=expression(sigma[proc]), xlim=c(0, 0.5), ylim=c(0,22))
   lines(density(di.poiss.IR$BUGSoutput$sims.matrix[,3]), lty=1, lwd=2,col="DarkGray")
   lines(density(di.poiss.SNF$BUGSoutput$sims.matrix[,3]), lty=2, lwd=2)
    
# Log-normal observation error: a, sig_p, sig_o
   plot(density(di.logN.D$BUGSoutput$sims.matrix[,1]),  lty=1, lwd=2,xlab="", ylab="Density", main=expression(italic(a)), ylim=c(0,26))
   lines(density(di.logN.IR$BUGSoutput$sims.matrix[,1]), lty=1, lwd=2,col="DarkGray")
   lines(density(di.logN.SNF$BUGSoutput$sims.matrix[,1]), lty=2, lwd=2)
  
  
   plot(density(di.logN.D$BUGSoutput$sims.matrix[,3]),  lty=1, lwd=2,xlab="", ylab="Density", main=expression(sigma[proc]), xlim=c(0, 0.5), ylim=c(0,22))
   lines(density(di.logN.IR$BUGSoutput$sims.matrix[,3]), lty=1, lwd=2,col="DarkGray")
   lines(density(di.logN.SNF$BUGSoutput$sims.matrix[,3]), lty=2, lwd=2)
   
   plot(density(di.logN.D$BUGSoutput$sims.matrix[,4]),  lty=1, lwd=2,xlab="", ylab="Density", main=expression(sigma[obs]), xlim=c(0, 0.5), ylim=c(0,22))
   lines(density(di.logN.IR$BUGSoutput$sims.matrix[,4]), lty=1, lwd=2,col="DarkGray")
   lines(density(di.logN.SNF$BUGSoutput$sims.matrix[,4]),lty=2, lwd=2)
    
  
  
  
#' # Ricker models
#' 
#' ## No observation error:
#' 
#'  $$N_{t+1}=N_t\exp(a+bN_t+\epsilon_t)$$
#'    
#'  $$\epsilon_t \sim N(0, \sigma_{proc}^2)$$  
#'   

#' Fit model in JAGS to Denali, Isle Royale, and SNF populations
  N<-(wdat[-c(1:28, 56), 2])
  x<-log(N)
  T<-length(N)
  Ricker.noerror.D <- jags(data=c("x","T"),   parameters.to.save=c("a","b", "sigma"),
    n.iter=25000, n.burnin=5000,model.file=Ricker.noerror, n.thin=1, progress.bar="none" )

  N<-(wdat[, 3])
  x<-log(N)
  T<-length(N)
  Ricker.noerror.IR <- jags(data=c("x","T"),   parameters.to.save=c("a","b", "sigma"),
    n.iter=25000, n.burnin=5000,model.file=Ricker.noerror, n.thin=1, progress.bar="none" )

  N<-(wdat[-c(1:8, 55:56), 4])
  x<-log(N)
  T<-length(N)
  Ricker.noerror.SNF <- jags(data=c("x","T"),   parameters.to.save=c("a","b", "sigma"),
    n.iter=25000, n.burnin=5000,model.file=Ricker.noerror, n.thin=1, progress.bar="none" )
  
  
  
#' ## Poisson observation error:
#' 
#'  $$N_{t+1}=N_t\exp(a+bN_t+\epsilon_t)$$
#'    
#'  $$\epsilon_t \sim N(0, \sigma_{proc}^2)$$  
#'  
#'  $$O_t \sim Poisson(N_t)$$
#'  
 

#' Fit model in JAGS to Denali, Isle Royale, and SNF populations
  O<-(wdat[-c(1:28, 56), 2])
  T<-length(O)
  Ricker.poiss.D <- jags(data=c("O","T"),   parameters.to.save=c("a","b", "sigma"),
    n.iter=25000, n.burnin=5000,model.file=Ricker.Poisson, n.thin=1, progress.bar="none" )

  O<-(wdat[, 3])
  T<-length(O)
  Ricker.poiss.IR <- jags(data=c("O","T"),   parameters.to.save=c("a", "b", "sigma"),
    n.iter=25000, n.burnin=5000,model.file=Ricker.Poisson, n.thin=1, progress.bar="none" )

  O<-(wdat[-c(1:8, 55:56), 4])
  T<-length(O)
  Ricker.poiss.SNF <- jags(data=c("O","T"),   parameters.to.save=c("a","b","sigma"),
    n.iter=25000, n.burnin=5000,model.file=Ricker.Poisson, n.thin=1, progress.bar="none" )
  
   
#' ## Log-normal observation error:
#' 
#'  $$N_{t+1}=N_t\exp(a+bN_t+\epsilon_t)$$ or $$X_{t+1} = X_t + a +bN_t + \epsilon_t$$
#'    
#'  $$\epsilon_t \sim N(0, \sigma_{proc}^2)$$  
#'  
#'  $$Y_t \sim N(X_t, \sigma_{obs})$$
#'  
 

#' Fit model in JAGS to Denali, Isle Royale, and SNF populations
  O<-(wdat[-c(1:28, 56), 2])
  y<-log(O)
  T<-length(O)
  Ricker.logN.D <- jags(data=c("y","T"),   parameters.to.save=c("a","b", "sigma", "tau"),
    n.iter=25000, n.burnin=5000,model.file=Ricker.LogN, n.thin=1, progress.bar="none" )

  O<-(wdat[, 3])
  y<-log(O)
  T<-length(O)
  Ricker.logN.IR <- jags(data=c("y","T"),   parameters.to.save=c("a","b", "sigma", "tau"),
    n.iter=25000, n.burnin=5000,model.file=Ricker.LogN, n.thin=1, progress.bar="none" )

  O<-(wdat[-c(1:8, 55:56), 4])
  y<-log(O)
  T<-length(O)
  Ricker.logN.SNF <- jags(data=c("y","T"),   parameters.to.save=c("a","b", "sigma", "tau"),
    n.iter=25000, n.burnin=5000,model.file=Ricker.LogN, n.thin=1, progress.bar="none" )
  
  
  
#' ### Plot posterior densities of mean.r and sig.r
#+ fig.height=8, fig.width=8
  layout(matrix(c(1,1,2,2,3,3,4,4,5,5,6,6,7,7,0,0,8,8,9,9,10,10,11,11), ncol=8, byrow=TRUE))
  par(bty="L", cex.main=1.8, cex.lab=1.6, cex.axis=1.4, oma=c(0,1,0,0))

# No observation error: a and sig_p
  plot(density(Ricker.noerror.D$BUGSoutput$sims.matrix[,1]), lty=1, lwd=2, xlab="", ylab="Density", main=expression(italic(a)),  ylim=c(0,6),xlim=c(-0.5, 1.5))
   lines(density(Ricker.noerror.IR$BUGSoutput$sims.matrix[,1]),lty=1, lwd=2,col="DarkGray")
   lines(density(Ricker.noerror.SNF$BUGSoutput$sims.matrix[,1]),lty=2, lwd=2)
  
   plot(density(Ricker.noerror.D$BUGSoutput$sims.matrix[,2]), lty=1, lwd=2,xlab="", ylab="Density", main=expression(italic(b)), xlim=c(-0.020, 0.005), ylim=c(0,290))
   lines(density(Ricker.noerror.IR$BUGSoutput$sims.matrix[,2]),lty=1, lwd=2,col="DarkGray")
   lines(density(Ricker.noerror.SNF$BUGSoutput$sims.matrix[,2]), lty=2, lwd=2)
  
  
   plot(density(Ricker.noerror.D$BUGSoutput$sims.matrix[,4]),lty=1, lwd=2,xlab="", ylab="Density", main=expression(sigma[proc]), xlim=c(0, 0.4),ylim=c(0,22))
   lines(density(Ricker.noerror.IR$BUGSoutput$sims.matrix[,4]), lty=1, lwd=2,col="DarkGray")
   lines(density(Ricker.noerror.SNF$BUGSoutput$sims.matrix[,4]), lty=2, lwd=2)

   plot(1:5, 1:5, type="n", xaxt="n", yaxt="n", bty="n", xlab="", ylab="") 
   legend("center", xjust=0.5,  lty=c(1,1,2), lwd=c(2,2,2), col=c("black", "DarkGray", "Black"),
          legend=c("Denali National Park", "Isle Royale","Superior National Forest" ),
          y.intersp=0.8, x.intersp=1.2,
          bty="n", cex=1.4, inset=4, xpd=TRUE) 
  
# Poisson observation error: a and sig_p  
   plot(density(Ricker.poiss.D$BUGSoutput$sims.matrix[,1]), lty=1, lwd=2,xlab="", ylab="Density", main=expression(italic(a)),  ylim=c(0,6),xlim=c(-0.5, 1.5))
   lines(density(Ricker.poiss.IR$BUGSoutput$sims.matrix[,1]), lty=1, lwd=2,col="DarkGray")
   lines(density(Ricker.poiss.SNF$BUGSoutput$sims.matrix[,1]), lty=2, lwd=2)
  
   plot(density(Ricker.poiss.D$BUGSoutput$sims.matrix[,2]), lty=1, lwd=2,xlab="", ylab="Density", main=expression(italic(b)),   xlim=c(-0.020, 0.005),ylim=c(0,290))
   lines(density(Ricker.poiss.IR$BUGSoutput$sims.matrix[,2]), lty=1, lwd=2,col="DarkGray")
   lines(density(Ricker.poiss.SNF$BUGSoutput$sims.matrix[,2]), lty=2, lwd=2)
  
  
   plot(density(Ricker.poiss.D$BUGSoutput$sims.matrix[,4]), lty=1, lwd=2,xlab="", ylab="Density", main=expression(sigma[proc]), xlim=c(0, 0.4),ylim=c(0,22))
   lines(density(Ricker.poiss.IR$BUGSoutput$sims.matrix[,4]), lty=1, lwd=2,col="DarkGray")
   lines(density(Ricker.poiss.SNF$BUGSoutput$sims.matrix[,4]), lty=2, lwd=2)
  
# Log-normal observation error: a, sig_p, sig_o
   plot(density(Ricker.logN.D$BUGSoutput$sims.matrix[,1]), lty=1, lwd=2,xlab="", ylab="Density", main=expression(italic(a)),  ylim=c(0,6), xlim=c(-0.5, 1.5))
   lines(density(Ricker.logN.IR$BUGSoutput$sims.matrix[,1]), lty=1, lwd=2,col="DarkGray")
   lines(density(Ricker.logN.SNF$BUGSoutput$sims.matrix[,1]),lty=2, lwd=2)
  
   plot(density(Ricker.logN.D$BUGSoutput$sims.matrix[,2]), lty=1, lwd=2,xlab="", ylab="Density", main=expression(italic(b)),   xlim=c(-0.020, 0.005),ylim=c(0,290))
   lines(density(Ricker.logN.IR$BUGSoutput$sims.matrix[,2]), lty=1, lwd=2,col="DarkGray")
   lines(density(Ricker.logN.SNF$BUGSoutput$sims.matrix[,2]), lty=2, lwd=2)
  
  
   plot(density(Ricker.logN.D$BUGSoutput$sims.matrix[,4]), lty=1, lwd=2,xlab="", ylab="Density", main=expression(sigma[proc]), xlim=c(0, 0.4),ylim=c(0,22))
   lines(density(Ricker.logN.IR$BUGSoutput$sims.matrix[,4]), lty=1, lwd=2,col="DarkGray")
   lines(density(Ricker.logN.SNF$BUGSoutput$sims.matrix[,4]), lty=2, lwd=2)
     
   plot(density(Ricker.logN.D$BUGSoutput$sims.matrix[,5]),lty=1, lwd=2,xlab="", ylab="Density", main=expression(sigma[obs]), xlim=c(0, 0.3),ylim=c(0,22))
   lines(density(Ricker.logN.IR$BUGSoutput$sims.matrix[,5]), lty=1, lwd=2,col="DarkGray")
   lines(density(Ricker.logN.SNF$BUGSoutput$sims.matrix[,5]), lty=2, lwd=2)
  
  
  
#' # Comparison of model fits by DIC
#' 
  dics<-matrix(NA, 6,3)
  dics[1,]<-c(di.noerror.D$BUGSoutput$DIC, di.noerror.IR$BUGSoutput$DIC, di.noerror.SNF$BUGSoutput$DIC)
  dics[2,]<-c(di.poiss.D$BUGSoutput$DIC, di.poiss.IR$BUGSoutput$DIC, di.poiss.SNF$BUGSoutput$DIC)
  dics[3,]<-c(di.logN.D$BUGSoutput$DIC, di.logN.IR$BUGSoutput$DIC, di.logN.SNF$BUGSoutput$DIC)
  dics[4,]<-c(Ricker.noerror.D$BUGSoutput$DIC, Ricker.noerror.IR$BUGSoutput$DIC, Ricker.noerror.SNF$BUGSoutput$DIC)
  dics[5,]<-c(Ricker.poiss.D$BUGSoutput$DIC, Ricker.poiss.IR$BUGSoutput$DIC, Ricker.poiss.SNF$BUGSoutput$DIC)
  dics[6,]<-c(Ricker.logN.D$BUGSoutput$DIC, Ricker.logN.IR$BUGSoutput$DIC, Ricker.logN.SNF$BUGSoutput$DIC)
  colnames(dics)<-c("Denali", "IR", "SNF")
  rownames(dics)<-c("DI.None", "DI.Poisson", "DI.LogN", "Ricker.None" , "Ricker.Poisson", "Ricker.LogN")
  dics
  
#' Conclusions:  Ricker model fits slightly better if no observation error is assumed to be present, the two
#' models are on equal footing under Poisson observation error, and the DI wins if log-nornmal observation error
#' is present.
#' 
  
#' # Plot data and fit of the models (use models without observation error)
  rickerfit<-function(fit, N, r, type, main){
    T<-length(N)
    plot(N[1:(T-1)], r[2:T], xlab="", ylab="", pch=16, main=main, adj=0)
    ns<-seq(min(N), max(N), length=100)
    lns<-length(ns)
    pns<-matrix(NA, lns, 3)
    for(i in 1:lns){
     pns.temp<-fit$BUGSoutput$sims.matrix[,1] + fit$BUGSoutput$sims.matrix[,2]*ns[i]  
     pns[i,1]<-mean(pns.temp)
     pns[i,2]<-quantile(pns.temp, probs=0.025)
     pns[i,3]<-quantile(pns.temp, probs=0.975)
    } 
    if(type=="ricker"){
    lines(ns, pns[,1], type="l")
    lines(ns, pns[,2], type="l", lty=2)
    lines(ns, pns[,3], type="l", lty=2)}
    if(type=="di"){
    abline(h=fit$BUGSoutput$mean[1])
    abline(h=fit$BUGSoutput$summary[1,3], lty=2)
    abline(h=fit$BUGSoutput$summary[1,7], lty=2)
     }
  }
  
par(mfcol=c(2,3), oma=c(1,2,0,0), mar=c(3.1, 3.1,3.1,2.1), cex.axis=1.6, cex.lab=1.6)  
  N<-(wdat[-c(1:28, 56), 2])
  r<-(wdat[-c(1:28, 56), 5])
  rickerfit(di.noerror.D, N,r, "di", main="A) Denali National Park")
  rickerfit(Ricker.noerror.D, N,r, "ricker", main="D)")

  N<-(wdat[, 3])
  r<-(wdat[, 6])
  rickerfit(di.noerror.IR, N, r, "di", main="B) Isle Royale")
  rickerfit(Ricker.noerror.IR, N,r, "ricker", main="E)") 
  
  r<-(wdat[-c(1:8, 55:56), 7])
  N<-(wdat[-c(1:8, 55:56), 4])
  rickerfit(di.noerror.SNF, N,r, "di", main="C) Superior National Forest")
  rickerfit(Ricker.noerror.SNF, N,r, "ricker", main="F)") 
  mtext(outer=T, side=2, line=-1, cex=1.4, expression(log(N[t+1]/N[t])))
  mtext(outer=T, side=1, line=0,cex=1.4, expression(N[t]))
  


  
#'  Another way to explore fit of Ricker Model  
#+fig.height=4, fig.width=8
  rickerfit2<-function(fit1, fit2, N, main){
    T<-length(N)
    plot(N[1:(T-1)], N[2:T], xlab="", ylab="", pch=16, main=main, adj=0)
    ns<-seq(min(N), max(N), length=100)
    lns<-length(ns)
    lsims1<-dim(fit1$BUGSoutput$sims.matrix)[1]
    lsims2<-dim(fit2$BUGSoutput$sims.matrix)[1]
    pns1<-pns2<-matrix(NA, lns, 3)
    for(i in 1:lns){
     pns.temp1<-ns[i]*exp(fit1$BUGSoutput$sims.matrix[,1] + rnorm(lsims1, 0, fit1$BUGSoutput$sims.matrix[,3]))  
     pns.temp2<-ns[i]*exp(fit2$BUGSoutput$sims.matrix[,1] + fit2$BUGSoutput$sims.matrix[,2]*ns[i] + rnorm(lsims2, 0, fit2$BUGSoutput$sims.matrix[,4]))
     pns1[i,1]<-mean(pns.temp1)
     pns1[i,2]<-quantile(pns.temp1, probs=0.025)
     pns1[i,3]<-quantile(pns.temp1, probs=0.975)
     pns2[i,1]<-mean(pns.temp2)
     pns2[i,2]<-quantile(pns.temp2, probs=0.025)
     pns2[i,3]<-quantile(pns.temp2, probs=0.975)
    } 
    lines(ns, pns1[,1], type="l")
    lines(ns, pns1[,2], type="l", lty=2)
    lines(ns, pns1[,3], type="l", lty=2)
    lines(ns, pns2[,1], type="l", col="gray", lwd=3)
    lines(ns, pns2[,2], type="l", col="gray", lty=2, lwd=3)
    lines(ns, pns2[,3], type="l", col="gray", lty=2, lwd=3)
  }
  
  par(mfrow=c(1,3),oma=c(2,3,0,0))
  N<-(wdat[-c(1:28, 56), 2])
  rickerfit2(di.noerror.D,Ricker.noerror.D, N, main="A) Denali National Park")
  
  N<-(wdat[, 3])
  rickerfit2(di.noerror.IR,Ricker.noerror.IR, N, main="B) Isle Royale")
  
  N<-(wdat[-c(1:8, 55:56), 4])
  rickerfit2(di.noerror.SNF,Ricker.noerror.SNF, N,main="C) Superior National Forest")
  
  mtext(outer=T, side=2, line=0, cex=1.8, expression(N[t+1]))
  mtext(outer=T, side=1, line=1,cex=1.8, expression(N[t]))
  
  
  
#' # Summarize results in 2 Tables for manuscript
#' ## Table 1. Fit of Density Independent Model (Bayesian implementation only)
#' 
#' 
  tab1<-matrix(NA, 7,9)
  colnames(tab1)<-rep(c("PE", "LCL", "UCL"), 3)
  rownames(tab1)<-c(rep("a",3), rep("sigmap",3), "sigobs")
  
# Pull off correct values from each model
# Point Estiamtes
    tab1[c(1,4),1]<-unlist(di.noerror.D$BUGSoutput$mean[c(1,3)])
    tab1[c(1,4),4]<- unlist(di.noerror.IR$BUGSoutput$mean[c(1,3)])
    tab1[c(1,4),7]<- unlist(di.noerror.SNF$BUGSoutput$mean[c(1,3)])
  
    tab1[c(2,5),1]<-unlist(di.poiss.D$BUGSoutput$mean[c(1,3)])
    tab1[c(2,5),4]<- unlist(di.poiss.IR$BUGSoutput$mean[c(1,3)])
    tab1[c(2,5),7]<- unlist(di.poiss.SNF$BUGSoutput$mean[c(1,3)])
  
    tab1[c(3,6),1]<-unlist(di.logN.D$BUGSoutput$mean[c(1,3)])
    tab1[c(3,6),4]<- unlist(di.logN.IR$BUGSoutput$mean[c(1,3)])
    tab1[c(3,6),7]<- unlist(di.logN.SNF$BUGSoutput$mean[c(1,3)])
  
# Conflimits  
  #no error  
    tab1[1,2:3]<- unlist(di.noerror.D$BUGSoutput$summary[1,c(3,7)])
    tab1[4,2:3]<- unlist(di.noerror.D$BUGSoutput$summary[3,c(3,7)])
  
    tab1[1,5:6]<- unlist(di.noerror.IR$BUGSoutput$summary[1,c(3,7)])
    tab1[4,5:6]<- unlist(di.noerror.IR$BUGSoutput$summary[3,c(3,7)])
  
    tab1[1,8:9]<- unlist(di.noerror.SNF$BUGSoutput$summary[1,c(3,7)])
    tab1[4,8:9]<- unlist(di.noerror.SNF$BUGSoutput$summary[3,c(3,7)])
  
  #Poisson error
    tab1[2,2:3]<- unlist(di.poiss.D$BUGSoutput$summary[1,c(3,7)])
    tab1[5,2:3]<- unlist(di.poiss.D$BUGSoutput$summary[3,c(3,7)])
  
    tab1[2,5:6]<- unlist(di.poiss.IR$BUGSoutput$summary[1,c(3,7)])
    tab1[5,5:6]<- unlist(di.poiss.IR$BUGSoutput$summary[3,c(3,7)])
  
    tab1[2,8:9]<- unlist(di.poiss.SNF$BUGSoutput$summary[1,c(3,7)])
    tab1[5,8:9]<- unlist(di.poiss.SNF$BUGSoutput$summary[3,c(3,7)])

  #lognormal error
    tab1[3,2:3]<- unlist(di.logN.D$BUGSoutput$summary[1,c(3,7)])
    tab1[6,2:3]<- unlist(di.logN.D$BUGSoutput$summary[3,c(3,7)])
  
    tab1[3,5:6]<- unlist(di.logN.IR$BUGSoutput$summary[1,c(3,7)])
    tab1[6,5:6]<- unlist(di.logN.IR$BUGSoutput$summary[3,c(3,7)])
  
    tab1[3,8:9]<- unlist(di.logN.SNF$BUGSoutput$summary[1,c(3,7)])
    tab1[6,8:9]<- unlist(di.logN.SNF$BUGSoutput$summary[3,c(3,7)])
  
  #Obs error
    tab1[7,1:3]<-unlist(di.logN.D$BUGSoutput$summary[4,c(1,3,7)])
    tab1[7,4:6]<-unlist(di.logN.IR$BUGSoutput$summary[4,c(1,3,7)])
    tab1[7,7:9]<-unlist(di.logN.SNF$BUGSoutput$summary[4,c(1,3,7)])

  
# Format the table
  tab1c<-matrix(NA, 7,6)
  tab1c[,1]<-as.character(round(tab1[,1],2))
  tab1c[,3]<-as.character(round(tab1[,4],2))
  tab1c[,5]<-as.character(round(tab1[,7],2))
  tab1c[,2]<-paste("(", round(tab1[,2],2), ", ", round(tab1[,3],2), ")", sep="")
  tab1c[,4]<-paste("(", round(tab1[,5],2), ", ", round(tab1[,6],2), ")", sep="")
  tab1c[,6]<-paste("(", round(tab1[,8],2), ", ", round(tab1[,9],2), ")", sep="")
  tab1c<-as.data.frame(tab1c)

  rownames(tab1c)<-c(paste(rep("a",3), c("none", "Poisson", "logN"), sep="."),  
                     paste(rep("sig_p",3), c("none", "Poisson", "logN"), sep="."), "sig_o.logN")
  colnames(tab1c)<-c("Denali.PE", "Denali.95%CI","IR.PE", "IR.95%CI","SNF.PE", "SNF.95%CI")
  tab1c
  
#' ## Table 2. Fit of Ricker Model (Bayesian implementation only)
#' 
#' 
  tab2<-matrix(NA, 10,9)
  colnames(tab2)<-rep(c("PE", "LCL", "UCL"), 3)
  rownames(tab2)<-c(rep("a",3),rep("b",3), rep("sigmap",3), "sigobs")
  
# Pull off correct values from each model
# Point Estiamtes
    tab2[c(1,4,7),1]<- unlist(Ricker.noerror.D$BUGSoutput$mean[c(1:2,4)])
    tab2[c(1,4,7),4]<- unlist(Ricker.noerror.IR$BUGSoutput$mean[c(1:2,4)])
    tab2[c(1,4, 7),7]<- unlist(Ricker.noerror.SNF$BUGSoutput$mean[c(1:2,4)])
  
    tab2[c(2,5,8),1]<-unlist(Ricker.poiss.D$BUGSoutput$mean[c(1:2,4)])
    tab2[c(2,5, 8),4]<- unlist(Ricker.poiss.IR$BUGSoutput$mean[c(1:2,4)])
    tab2[c(2,5, 8),7]<- unlist(Ricker.poiss.SNF$BUGSoutput$mean[c(1:2,4)])
  
    tab2[c(3,6,9),1]<-unlist(Ricker.logN.D$BUGSoutput$mean[c(1:2,4)])
    tab2[c(3,6, 9),4]<- unlist(Ricker.logN.IR$BUGSoutput$mean[c(1:2,4)])
    tab2[c(3,6, 9),7]<- unlist(Ricker.logN.SNF$BUGSoutput$mean[c(1:2,4)])
  
# Conflimits  
  #no error  
    tab2[1,2:3]<- unlist(Ricker.noerror.D$BUGSoutput$summary[1,c(3,7)])
    tab2[4,2:3]<- unlist(Ricker.noerror.D$BUGSoutput$summary[2,c(3,7)])
    tab2[7,2:3]<- unlist(Ricker.noerror.D$BUGSoutput$summary[4,c(3,7)])
  
    tab2[1,5:6]<- unlist(Ricker.noerror.IR$BUGSoutput$summary[1,c(3,7)])
    tab2[4,5:6]<- unlist(Ricker.noerror.IR$BUGSoutput$summary[2,c(3,7)])
    tab2[7,5:6]<- unlist(Ricker.noerror.IR$BUGSoutput$summary[4,c(3,7)])
  
    tab2[1,8:9]<- unlist(Ricker.noerror.SNF$BUGSoutput$summary[1,c(3,7)])
    tab2[4,8:9]<- unlist(Ricker.noerror.SNF$BUGSoutput$summary[2,c(3,7)])
    tab2[7,8:9]<- unlist(Ricker.noerror.SNF$BUGSoutput$summary[4,c(3,7)])
   
  #Poisson error
    tab2[2,2:3]<- unlist(Ricker.poiss.D$BUGSoutput$summary[1,c(3,7)])
    tab2[5,2:3]<- unlist(Ricker.poiss.D$BUGSoutput$summary[2,c(3,7)])
    tab2[8,2:3]<- unlist(Ricker.poiss.D$BUGSoutput$summary[4,c(3,7)])
  
    tab2[2,5:6]<- unlist(Ricker.poiss.IR$BUGSoutput$summary[1,c(3,7)])
    tab2[5,5:6]<- unlist(Ricker.poiss.IR$BUGSoutput$summary[2,c(3,7)])
    tab2[8,5:6]<- unlist(Ricker.poiss.IR$BUGSoutput$summary[4,c(3,7)])
  
    tab2[2,8:9]<- unlist(Ricker.poiss.SNF$BUGSoutput$summary[1,c(3,7)])
    tab2[5,8:9]<- unlist(Ricker.poiss.SNF$BUGSoutput$summary[2,c(3,7)])
    tab2[8,8:9]<- unlist(Ricker.poiss.SNF$BUGSoutput$summary[4,c(3,7)])

  #lognormal error
    tab2[3,2:3]<- unlist(Ricker.logN.D$BUGSoutput$summary[1,c(3,7)])
    tab2[6,2:3]<- unlist(Ricker.logN.D$BUGSoutput$summary[2,c(3,7)])
    tab2[9,2:3]<- unlist(Ricker.logN.D$BUGSoutput$summary[4,c(3,7)])
  
    tab2[3,5:6]<- unlist(Ricker.logN.IR$BUGSoutput$summary[1,c(3,7)])
    tab2[6,5:6]<- unlist(Ricker.logN.IR$BUGSoutput$summary[2,c(3,7)])
    tab2[9,5:6]<- unlist(Ricker.logN.IR$BUGSoutput$summary[4,c(3,7)])
  
  
    tab2[3,8:9]<- unlist(Ricker.logN.SNF$BUGSoutput$summary[1,c(3,7)])
    tab2[6,8:9]<- unlist(Ricker.logN.SNF$BUGSoutput$summary[2,c(3,7)])
    tab2[9,8:9]<- unlist(Ricker.logN.SNF$BUGSoutput$summary[4,c(3,7)])
  
  #Obs error
    tab2[10,1:3]<-unlist(Ricker.logN.D$BUGSoutput$summary[5,c(1,3,7)])
    tab2[10,4:6]<-unlist(Ricker.logN.IR$BUGSoutput$summary[5,c(1,3,7)])
    tab2[10,7:9]<-unlist(Ricker.logN.SNF$BUGSoutput$summary[5,c(1,3,7)])

  
# Format the table
  tab2c<-matrix(NA, 10,6)
  tab2c[c(1:3,7:10),1]<-as.character(round(tab2[c(1:3,7:10),1],2))
  tab2c[c(1:3,7:10),3]<-as.character(round(tab2[c(1:3,7:10),4],2))
  tab2c[c(1:3,7:10),5]<-as.character(round(tab2[c(1:3,7:10),7],2))
  tab2c[c(1:3,7:10),2]<-paste("(", round(tab2[c(1:3,7:10),2],2), ", ", round(tab2[c(1:3,7:10),3],2), ")", sep="")
  tab2c[c(1:3,7:10),4]<-paste("(", round(tab2[c(1:3,7:10),5],2), ", ", round(tab2[c(1:3,7:10),6],2), ")", sep="")
  tab2c[c(1:3,7:10),6]<-paste("(", round(tab2[c(1:3,7:10),8],2), ", ", round(tab2[c(1:3,7:10),9],2), ")", sep="")
  tab2c[4:6,1]<-as.character(round(tab2[4:6,1],3))
  tab2c[4:6,3]<-as.character(round(tab2[4:6,4],3))
  tab2c[4:6,5]<-as.character(round(tab2[4:6,7],3))
  tab2c[4:6,2]<-paste("(", round(tab2[4:6,2],3), ", ", round(tab2[4:6,3],2), ")", sep="")
  tab2c[4:6,4]<-paste("(", round(tab2[4:6,5],3), ", ", round(tab2[4:6,6],2), ")", sep="")
  tab2c[4:6,6]<-paste("(", round(tab2[4:6,8],3), ", ", round(tab2[4:6,9],2), ")", sep="")
  tab2c<-as.data.frame(tab2c)
  rownames(tab2c)<-c(paste(rep("a",3), c("none", "Poisson", "logN"), sep="."), paste(rep("b",3), c("none", "Poisson", "logN"), sep="."), 
                     paste(rep("sig_p",3), c("none", "Poisson", "logN"), sep="."), "sig_o.logN")
  colnames(tab2c)<-c("Denali.PE", "Denali.95%CI","IR.PE", "IR.95%CI","SNF.PE", "SNF.95%CI")
  tab2c
  
#' # Now run frequentist analysis 
#+ pvaclones, include=FALSE  
#  Models for population 1
  fit1a<-pva(wdat$Denali[-c(1:27, 56)], ricker("none", fixed=c(b=0)), 20,  n.iter=25000, progress.bar="none")
  fit1b<-pva(wdat$Denali[-c(1:27, 56)], ricker("none"), 20,  n.iter=25000) #ricker w/ density depend
  fit1c<-pva(wdat$Denali[-c(1:27, 56)], ricker("normal", fixed=c(b=0)), 20,  n.iter=25000)
  fit1d<-pva(wdat$Denali[-c(1:27, 56)], ricker("normal"), 20,  n.iter=25000) #ricker w/ density depend
  fit1e<-pva(wdat$Denali[-c(1:27, 56)], ricker("poisson", fixed=c(b=0)), 20,  n.iter=25000)
  fit1f<-pva(wdat$Denali[-c(1:27, 56)], ricker("poisson"), 20,  n.iter=25000) #ricker w/ density depend
   
# Models for population 2   
  fit2a<-pva(wdat$IsleRoyale, ricker("none", fixed=c(b=0)), 20,  n.iter=25000)
  fit2b<-pva(wdat$IsleRoyale, ricker("none"), 20,  n.iter=25000) #ricker w/ density depend
  fit2c<-pva(wdat$IsleRoyale, ricker("normal", fixed=c(b=0)), 20,  n.iter=25000)
  fit2d<-pva(wdat$IsleRoyale, ricker("normal"), 20,  n.iter=25000) #ricker w/ density depend
  fit2e<-pva(wdat$IsleRoyale, ricker("poisson", fixed=c(b=0)), 20,  n.iter=25000)
  fit2f<-pva(wdat$IsleRoyale, ricker("poisson"), 20,  n.iter=25000) #ricker w/ density depend

# Models for Population 3  
  fit3a<-pva(wdat$SNF[-c(1:8, 55:56)], ricker("none", fixed=c(b=0)), 20,  n.iter=25000)
  fit3b<-pva(wdat$SNF[-c(1:8, 55:56)], ricker("none"), 20,  n.iter=25000) #ricker w/ density depend
  fit3c<-pva(wdat$SNF[-c(1:8, 55:56)], ricker("normal", fixed=c(b=0)), 20,  n.iter=25000)
  fit3d<-pva(wdat$SNF[-c(1:8, 55:56)], ricker("normal"), 20,  n.iter=25000) #ricker w/ density depend
  fit3e<-pva(wdat$SNF[-c(1:8, 55:56)], ricker("poisson", fixed=c(b=0)), 20,  n.iter=25000)
  fit3f<-pva(wdat$SNF[-c(1:8, 55:56)], ricker("poisson"), 20,  n.iter=25000) #ricker w/ density depend
    
#+ pvaclonedummy, ref.label = 'pvaclones', eval = FALSE, echo = TRUE
#  Models for population 1
  fit1a<-pva(wdat$Denali[-c(1:27, 56)], ricker("none", fixed=c(b=0)), 20,  n.iter=25000, progress.bar="none")
  fit1b<-pva(wdat$Denali[-c(1:27, 56)], ricker("none"), 20,  n.iter=25000) #ricker w/ density depend
  fit1c<-pva(wdat$Denali[-c(1:27, 56)], ricker("normal", fixed=c(b=0)), 20,  n.iter=25000)
  fit1d<-pva(wdat$Denali[-c(1:27, 56)], ricker("normal"), 20,  n.iter=25000) #ricker w/ density depend
  fit1e<-pva(wdat$Denali[-c(1:27, 56)], ricker("poisson", fixed=c(b=0)), 20,  n.iter=25000)
  fit1f<-pva(wdat$Denali[-c(1:27, 56)], ricker("poisson"), 20,  n.iter=25000) #ricker w/ density depend
   
# Models for population 2   
  fit2a<-pva(wdat$IsleRoyale, ricker("none", fixed=c(b=0)), 20,  n.iter=25000)
  fit2b<-pva(wdat$IsleRoyale, ricker("none"), 20,  n.iter=25000) #ricker w/ density depend
  fit2c<-pva(wdat$IsleRoyale, ricker("normal", fixed=c(b=0)), 20,  n.iter=25000)
  fit2d<-pva(wdat$IsleRoyale, ricker("normal"), 20,  n.iter=25000) #ricker w/ density depend
  fit2e<-pva(wdat$IsleRoyale, ricker("poisson", fixed=c(b=0)), 20,  n.iter=25000)
  fit2f<-pva(wdat$IsleRoyale, ricker("poisson"), 20,  n.iter=25000) #ricker w/ density depend

# Models for Population 3  
  fit3a<-pva(wdat$SNF[-c(1:8, 55:56)], ricker("none", fixed=c(b=0)), 20,  n.iter=25000)
  fit3b<-pva(wdat$SNF[-c(1:8, 55:56)], ricker("none"), 20,  n.iter=25000) #ricker w/ density depend
  fit3c<-pva(wdat$SNF[-c(1:8, 55:56)], ricker("normal", fixed=c(b=0)), 20,  n.iter=25000)
  fit3d<-pva(wdat$SNF[-c(1:8, 55:56)], ricker("normal"), 20,  n.iter=25000) #ricker w/ density depend
  fit3e<-pva(wdat$SNF[-c(1:8, 55:56)], ricker("poisson", fixed=c(b=0)), 20,  n.iter=25000)
  fit3f<-pva(wdat$SNF[-c(1:8, 55:56)], ricker("poisson"), 20,  n.iter=25000) #ricker w/ density depend
        
#' # Summarize results in 2 Tables for manuscript
#' ## Table 1. Fit of Density Independent Model (Bayesian implementation only)
#' 
#' 
  tab1<-matrix(NA, 7,9)
  colnames(tab1)<-rep(c("PE", "LCL", "UCL"), 3)
  rownames(tab1)<-c(rep("a",3), rep("sigmap",3), "sigobs",rep("DIC",3))
  
# Pull off correct values from each model
# Point Estiamtes
    tab1[c(1,4),1]<-unlist(di.noerror.D$BUGSoutput$mean[c(1,3)])
    tab1[c(1,4),4]<- unlist(di.noerror.IR$BUGSoutput$mean[c(1,3)])
    tab1[c(1,4),7]<- unlist(di.noerror.SNF$BUGSoutput$mean[c(1,3)])
  
    tab1[c(2,5),1]<-unlist(di.poiss.D$BUGSoutput$mean[c(1,3)])
    tab1[c(2,5),4]<- unlist(di.poiss.IR$BUGSoutput$mean[c(1,3)])
    tab1[c(2,5),7]<- unlist(di.poiss.SNF$BUGSoutput$mean[c(1,3)])
  
    tab1[c(3,6),1]<-unlist(di.logN.D$BUGSoutput$mean[c(1,3)])
    tab1[c(3,6),4]<- unlist(di.logN.IR$BUGSoutput$mean[c(1,3)])
    tab1[c(3,6),7]<- unlist(di.logN.SNF$BUGSoutput$mean[c(1,3)])
  
# Conflimits  
  #no error  
    tab1[1,2:3]<- unlist(di.noerror.D$BUGSoutput$summary[1,c(3,7)])
    tab1[4,2:3]<- unlist(di.noerror.D$BUGSoutput$summary[3,c(3,7)])
  
    tab1[1,5:6]<- unlist(di.noerror.IR$BUGSoutput$summary[1,c(3,7)])
    tab1[4,5:6]<- unlist(di.noerror.IR$BUGSoutput$summary[3,c(3,7)])
  
    tab1[1,8:9]<- unlist(di.noerror.SNF$BUGSoutput$summary[1,c(3,7)])
    tab1[4,8:9]<- unlist(di.noerror.SNF$BUGSoutput$summary[3,c(3,7)])
  
  #Poisson error
    tab1[2,2:3]<- unlist(di.poiss.D$BUGSoutput$summary[1,c(3,7)])
    tab1[5,2:3]<- unlist(di.poiss.D$BUGSoutput$summary[3,c(3,7)])
  
    tab1[2,5:6]<- unlist(di.poiss.IR$BUGSoutput$summary[1,c(3,7)])
    tab1[5,5:6]<- unlist(di.poiss.IR$BUGSoutput$summary[3,c(3,7)])
  
    tab1[2,8:9]<- unlist(di.poiss.SNF$BUGSoutput$summary[1,c(3,7)])
    tab1[5,8:9]<- unlist(di.poiss.SNF$BUGSoutput$summary[3,c(3,7)])

  #lognormal error
    tab1[3,2:3]<- unlist(di.logN.D$BUGSoutput$summary[1,c(3,7)])
    tab1[6,2:3]<- unlist(di.logN.D$BUGSoutput$summary[3,c(3,7)])
  
    tab1[3,5:6]<- unlist(di.logN.IR$BUGSoutput$summary[1,c(3,7)])
    tab1[6,5:6]<- unlist(di.logN.IR$BUGSoutput$summary[3,c(3,7)])
  
    tab1[3,8:9]<- unlist(di.logN.SNF$BUGSoutput$summary[1,c(3,7)])
    tab1[6,8:9]<- unlist(di.logN.SNF$BUGSoutput$summary[3,c(3,7)])
  
  #Obs error
    tab1[7,1:3]<-unlist(di.logN.D$BUGSoutput$summary[4,c(1,3,7)])
    tab1[7,4:6]<-unlist(di.logN.IR$BUGSoutput$summary[4,c(1,3,7)])
    tab1[7,7:9]<-unlist(di.logN.SNF$BUGSoutput$summary[4,c(1,3,7)])

  
# Format the table
  tab1c<-matrix(NA, 7,6)
  tab1c[,1]<-as.character(round(tab1[,1],2))
  tab1c[,3]<-as.character(round(tab1[,4],2))
  tab1c[,5]<-as.character(round(tab1[,7],2))
  tab1c[,2]<-paste("(", round(tab1[,2],2), ", ", round(tab1[,3],2), ")", sep="")
  tab1c[,4]<-paste("(", round(tab1[,5],2), ", ", round(tab1[,6],2), ")", sep="")
  tab1c[,6]<-paste("(", round(tab1[,8],2), ", ", round(tab1[,9],2), ")", sep="")
  as.data.frame(tab1c)

  
  
#' ## Table 1. Fit of Density Independent Model (Bayesian implementation only)
#' 
#' 
  tab1<-matrix(NA, 7,9)
  colnames(tab1)<-rep(c("PE", "LCL", "UCL"), 3)
  rownames(tab1)<-c(rep("a",3), rep("sigmap",3), "sigobs")
  
# Pull off correct values from each model
# Point Estiamtes
    tab1[c(1,4),1]<-unlist(di.noerror.D$BUGSoutput$mean[c(1,3)])
    tab1[c(1,4),4]<- unlist(di.noerror.IR$BUGSoutput$mean[c(1,3)])
    tab1[c(1,4),7]<- unlist(di.noerror.SNF$BUGSoutput$mean[c(1,3)])
  
    tab1[c(2,5),1]<-unlist(di.poiss.D$BUGSoutput$mean[c(1,3)])
    tab1[c(2,5),4]<- unlist(di.poiss.IR$BUGSoutput$mean[c(1,3)])
    tab1[c(2,5),7]<- unlist(di.poiss.SNF$BUGSoutput$mean[c(1,3)])
  
    tab1[c(3,6),1]<-unlist(di.logN.D$BUGSoutput$mean[c(1,3)])
    tab1[c(3,6),4]<- unlist(di.logN.IR$BUGSoutput$mean[c(1,3)])
    tab1[c(3,6),7]<- unlist(di.logN.SNF$BUGSoutput$mean[c(1,3)])
  
# Conflimits  
  #no error  
    tab1[1,2:3]<- unlist(di.noerror.D$BUGSoutput$summary[1,c(3,7)])
    tab1[4,2:3]<- unlist(di.noerror.D$BUGSoutput$summary[3,c(3,7)])
  
    tab1[1,5:6]<- unlist(di.noerror.IR$BUGSoutput$summary[1,c(3,7)])
    tab1[4,5:6]<- unlist(di.noerror.IR$BUGSoutput$summary[3,c(3,7)])
  
    tab1[1,8:9]<- unlist(di.noerror.SNF$BUGSoutput$summary[1,c(3,7)])
    tab1[4,8:9]<- unlist(di.noerror.SNF$BUGSoutput$summary[3,c(3,7)])
  
  #Poisson error
    tab1[2,2:3]<- unlist(di.poiss.D$BUGSoutput$summary[1,c(3,7)])
    tab1[5,2:3]<- unlist(di.poiss.D$BUGSoutput$summary[3,c(3,7)])
  
    tab1[2,5:6]<- unlist(di.poiss.IR$BUGSoutput$summary[1,c(3,7)])
    tab1[5,5:6]<- unlist(di.poiss.IR$BUGSoutput$summary[3,c(3,7)])
  
    tab1[2,8:9]<- unlist(di.poiss.SNF$BUGSoutput$summary[1,c(3,7)])
    tab1[5,8:9]<- unlist(di.poiss.SNF$BUGSoutput$summary[3,c(3,7)])

  #lognormal error
    tab1[3,2:3]<- unlist(di.logN.D$BUGSoutput$summary[1,c(3,7)])
    tab1[6,2:3]<- unlist(di.logN.D$BUGSoutput$summary[3,c(3,7)])
  
    tab1[3,5:6]<- unlist(di.logN.IR$BUGSoutput$summary[1,c(3,7)])
    tab1[6,5:6]<- unlist(di.logN.IR$BUGSoutput$summary[3,c(3,7)])
  
    tab1[3,8:9]<- unlist(di.logN.SNF$BUGSoutput$summary[1,c(3,7)])
    tab1[6,8:9]<- unlist(di.logN.SNF$BUGSoutput$summary[3,c(3,7)])
  
  #Obs error
    tab1[7,1:3]<-unlist(di.logN.D$BUGSoutput$summary[4,c(1,3,7)])
    tab1[7,4:6]<-unlist(di.logN.IR$BUGSoutput$summary[4,c(1,3,7)])
    tab1[7,7:9]<-unlist(di.logN.SNF$BUGSoutput$summary[4,c(1,3,7)])

  
# Format the table
  tab1c<-matrix(NA, 7,6)
  tab1c[,1]<-as.character(round(tab1[,1],2))
  tab1c[,3]<-as.character(round(tab1[,4],2))
  tab1c[,5]<-as.character(round(tab1[,7],2))
  tab1c[,2]<-paste("(", round(tab1[,2],2), ", ", round(tab1[,3],2), ")", sep="")
  tab1c[,4]<-paste("(", round(tab1[,5],2), ", ", round(tab1[,6],2), ")", sep="")
  tab1c[,6]<-paste("(", round(tab1[,8],2), ", ", round(tab1[,9],2), ")", sep="")
  as.data.frame(tab1c)
  
#' ## Table 2. Fit of Ricker Model (Bayesian implementation only)
#' 
#' 
  tab2<-matrix(NA, 10,9)
  colnames(tab2)<-rep(c("PE", "LCL", "UCL"), 3)
  rownames(tab2)<-c(rep("a",3),rep("b",3), rep("sigmap",3), "sigobs")
  
# Pull off correct values from each model
# Point Estiamtes
    tab2[c(1,4,7),1]<- unlist(Ricker.noerror.D$BUGSoutput$mean[c(1:2,4)])
    tab2[c(1,4,7),4]<- unlist(Ricker.noerror.IR$BUGSoutput$mean[c(1:2,4)])
    tab2[c(1,4, 7),7]<- unlist(Ricker.noerror.SNF$BUGSoutput$mean[c(1:2,4)])
  
    tab2[c(2,5,8),1]<-unlist(Ricker.poiss.D$BUGSoutput$mean[c(1:2,4)])
    tab2[c(2,5, 8),4]<- unlist(Ricker.poiss.IR$BUGSoutput$mean[c(1:2,4)])
    tab2[c(2,5, 8),7]<- unlist(Ricker.poiss.SNF$BUGSoutput$mean[c(1:2,4)])
  
    tab2[c(3,6,9),1]<-unlist(Ricker.logN.D$BUGSoutput$mean[c(1:2,4)])
    tab2[c(3,6, 9),4]<- unlist(Ricker.logN.IR$BUGSoutput$mean[c(1:2,4)])
    tab2[c(3,6, 9),7]<- unlist(Ricker.logN.SNF$BUGSoutput$mean[c(1:2,4)])
  
# Conflimits  
  #no error  
    tab2[1,2:3]<- unlist(Ricker.noerror.D$BUGSoutput$summary[1,c(3,7)])
    tab2[4,2:3]<- unlist(Ricker.noerror.D$BUGSoutput$summary[2,c(3,7)])
    tab2[7,2:3]<- unlist(Ricker.noerror.D$BUGSoutput$summary[4,c(3,7)])
  
    tab2[1,5:6]<- unlist(Ricker.noerror.IR$BUGSoutput$summary[1,c(3,7)])
    tab2[4,5:6]<- unlist(Ricker.noerror.IR$BUGSoutput$summary[2,c(3,7)])
    tab2[7,5:6]<- unlist(Ricker.noerror.IR$BUGSoutput$summary[4,c(3,7)])
  
    tab2[1,8:9]<- unlist(Ricker.noerror.SNF$BUGSoutput$summary[1,c(3,7)])
    tab2[4,8:9]<- unlist(Ricker.noerror.SNF$BUGSoutput$summary[2,c(3,7)])
    tab2[7,8:9]<- unlist(Ricker.noerror.SNF$BUGSoutput$summary[4,c(3,7)])
   
  #Poisson error
    tab2[2,2:3]<- unlist(Ricker.poiss.D$BUGSoutput$summary[1,c(3,7)])
    tab2[5,2:3]<- unlist(Ricker.poiss.D$BUGSoutput$summary[2,c(3,7)])
    tab2[8,2:3]<- unlist(Ricker.poiss.D$BUGSoutput$summary[4,c(3,7)])
  
    tab2[2,5:6]<- unlist(Ricker.poiss.IR$BUGSoutput$summary[1,c(3,7)])
    tab2[5,5:6]<- unlist(Ricker.poiss.IR$BUGSoutput$summary[2,c(3,7)])
    tab2[8,5:6]<- unlist(Ricker.poiss.IR$BUGSoutput$summary[4,c(3,7)])
  
    tab2[2,8:9]<- unlist(Ricker.poiss.SNF$BUGSoutput$summary[1,c(3,7)])
    tab2[5,8:9]<- unlist(Ricker.poiss.SNF$BUGSoutput$summary[2,c(3,7)])
    tab2[8,8:9]<- unlist(Ricker.poiss.SNF$BUGSoutput$summary[4,c(3,7)])

  #lognormal error
    tab2[3,2:3]<- unlist(Ricker.logN.D$BUGSoutput$summary[1,c(3,7)])
    tab2[6,2:3]<- unlist(Ricker.logN.D$BUGSoutput$summary[2,c(3,7)])
    tab2[9,2:3]<- unlist(Ricker.logN.D$BUGSoutput$summary[4,c(3,7)])
  
    tab2[3,5:6]<- unlist(Ricker.logN.IR$BUGSoutput$summary[1,c(3,7)])
    tab2[6,5:6]<- unlist(Ricker.logN.IR$BUGSoutput$summary[2,c(3,7)])
    tab2[9,5:6]<- unlist(Ricker.logN.IR$BUGSoutput$summary[4,c(3,7)])
  
  
    tab2[3,8:9]<- unlist(Ricker.logN.SNF$BUGSoutput$summary[1,c(3,7)])
    tab2[6,8:9]<- unlist(Ricker.logN.SNF$BUGSoutput$summary[2,c(3,7)])
    tab2[9,8:9]<- unlist(Ricker.logN.SNF$BUGSoutput$summary[4,c(3,7)])
  
  #Obs error
    tab2[10,1:3]<-unlist(Ricker.logN.D$BUGSoutput$summary[5,c(1,3,7)])
    tab2[10,4:6]<-unlist(Ricker.logN.IR$BUGSoutput$summary[5,c(1,3,7)])
    tab2[10,7:9]<-unlist(Ricker.logN.SNF$BUGSoutput$summary[5,c(1,3,7)])

  
# Format the table
  tab2c<-matrix(NA, 10,6)
  tab2c[c(1:3,7:10),1]<-as.character(round(tab2[c(1:3,7:10),1],2))
  tab2c[c(1:3,7:10),3]<-as.character(round(tab2[c(1:3,7:10),4],2))
  tab2c[c(1:3,7:10),5]<-as.character(round(tab2[c(1:3,7:10),7],2))
  tab2c[c(1:3,7:10),2]<-paste("(", round(tab2[c(1:3,7:10),2],2), ", ", round(tab2[c(1:3,7:10),3],2), ")", sep="")
  tab2c[c(1:3,7:10),4]<-paste("(", round(tab2[c(1:3,7:10),5],2), ", ", round(tab2[c(1:3,7:10),6],2), ")", sep="")
  tab2c[c(1:3,7:10),6]<-paste("(", round(tab2[c(1:3,7:10),8],2), ", ", round(tab2[c(1:3,7:10),9],2), ")", sep="")
  tab2c[4:6,1]<-as.character(round(tab2[4:6,1],3))
  tab2c[4:6,3]<-as.character(round(tab2[4:6,4],3))
  tab2c[4:6,5]<-as.character(round(tab2[4:6,7],3))
  tab2c[4:6,2]<-paste("(", round(tab2[4:6,2],3), ", ", round(tab2[4:6,3],2), ")", sep="")
  tab2c[4:6,4]<-paste("(", round(tab2[4:6,5],3), ", ", round(tab2[4:6,6],2), ")", sep="")
  tab2c[4:6,6]<-paste("(", round(tab2[4:6,8],3), ", ", round(tab2[4:6,9],2), ")", sep="")
  as.data.frame(tab2c)

  
  
#' ###  Table 3: Sensitivity analysis, Frequentist fits to Density independent models
 tab3<-matrix(NA, 7,6)
  colnames(tab3)<-rep(c("Bayes", "Freq"), 3)
  rownames(tab3)<-c(rep("a",3),rep("sigmap",3),"Sigobs")
  
# population table with Bayesian values from table 1
  tab3[,c(1,3,5)]<-tab1[,c(1,4,7)]
  tab3[c(1,4),2]<-coef(fit1a)[c(1,3)]
  tab3[c(3,6,7),2]<-coef(fit1c)[c(1,3,4)]  
  tab3[c(2,5),2]<-coef(fit1e)[c(1,3)]  
  
  tab3[c(1,4),4]<-coef(fit2a)[c(1,3)]
  tab3[c(3,6,7),4]<-coef(fit2c)[c(1,3,4)]  
  tab3[c(2,5),4]<-coef(fit2e)[c(1,3)]  
  
  tab3[c(1,4),6]<-coef(fit3a)[c(1,3)]
  tab3[c(3,6,7),6]<-coef(fit3c)[c(1,3,4)]  
  tab3[c(2,5),6]<-coef(fit3e)[c(1,3)]  
  
#' ###  Table 4: Sensitivity analysis, Frequentist fits to Ricker models  
   tab4<-matrix(NA, 10,6)
  colnames(tab4)<-rep(c("Bayes", "Freq"), 3)
  rownames(tab4)<-c(rep("a",3),rep("b",3),rep("sigmap",3),"Sigobs")
  
# population table with Bayesian values from table 1
  tab4[,c(1,3,5)]<-tab2[,c(1,4,7)]
  tab4[c(1,4,7),2]<-coef(fit1b)[c(1:3)]
  tab4[c(3,6,9,10),2]<-coef(fit1d)[c(1:4)]  
  tab4[c(2,5,8),2]<-coef(fit1f)[c(1:3)]  
  
  tab4[c(1,4,7),4]<-coef(fit2b)[c(1:3)]
  tab4[c(3,6,9,10),4]<-coef(fit2d)[c(1:4)]  
  tab4[c(2,5,8),4]<-coef(fit2f)[c(1:3)]  
 
  tab4[c(1,4,7),6]<-coef(fit3b)[c(1:3)]
  tab4[c(3,6,9,10),6]<-coef(fit3d)[c(1:4)]  
  tab4[c(2,5,8),6]<-coef(fit3f)[c(1:3)]  
  
tab4b<-tab4
round(tab4b,3)