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
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:
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
## r.bar SE.r.bar Sig2.r n
## Denali -0.003835 0.03912 0.04133 27
## IsleRoyale -0.01452 0.03501 0.06742 55
## SNF 0.005758 0.02548 0.02921 45
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
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)))
## Warning: Removed 38 rows containing missing values (geom_path).
## Warning: Removed 38 rows containing missing values (geom_point).
Yearly growth rates over time
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]))
## Warning: Removed 41 rows containing missing values (geom_path).
## Warning: Removed 41 rows containing missing values (geom_point).
Distribution of yearly growth rates
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)))
## Warning: Removed 29 rows containing non-finite values (stat_density).
## Warning: Removed 1 rows containing non-finite values (stat_density).
## Warning: Removed 11 rows containing non-finite values (stat_density).
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')
\[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" )
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 102
##
## Initializing model
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" )
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 148
##
## Initializing model
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" )
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 142
##
## Initializing model
\[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" )
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 143
##
## Initializing model
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" )
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 288
##
## Initializing model
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" )
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 238
##
## Initializing model
\[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" )
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 147
##
## Initializing model
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" )
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 292
##
## Initializing model
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" )
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 242
##
## Initializing model
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)
\[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" )
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 125
##
## Initializing model
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" )
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 177
##
## Initializing model
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" )
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 172
##
## Initializing model
\[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" )
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 171
##
## Initializing model
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" )
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 345
##
## Initializing model
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" )
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 285
##
## Initializing model
\[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" )
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 175
##
## Initializing model
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" )
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 349
##
## Initializing model
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" )
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 289
##
## Initializing model
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)
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
## Denali IR SNF
## DI.None -4.831 11.060 -27.83
## DI.Poisson 206.570 335.664 321.59
## DI.LogN 152.345 778.672 124.02
## Ricker.None -7.407 8.983 -29.06
## Ricker.Poisson 209.212 337.258 321.86
## Ricker.LogN 225.919 958.944 225.79
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.
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
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]))
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
## Denali.PE Denali.95%CI IR.PE IR.95%CI SNF.PE
## a.none 0 (-0.08, 0.09) -0.01 (-0.09, 0.06) 0.01
## a.Poisson 0 (-0.08, 0.08) -0.01 (-0.06, 0.04) 0.01
## a.logN 0 (-0.08, 0.08) -0.01 (-0.08, 0.05) 0.01
## sig_p.none 0.22 (0.17, 0.29) 0.27 (0.22, 0.32) 0.18
## sig_p.Poisson 0.2 (0.14, 0.28) 0.19 (0.13, 0.26) 0.1
## sig_p.logN 0.2 (0.14, 0.29) 0.23 (0.16, 0.3) 0.12
## sig_o.logN 0.09 (0.03, 0.17) 0.11 (0.04, 0.19) 0.09
## SNF.95%CI
## a.none (-0.05, 0.06)
## a.Poisson (-0.03, 0.04)
## a.logN (-0.03, 0.04)
## sig_p.none (0.14, 0.22)
## sig_p.Poisson (0.07, 0.15)
## sig_p.logN (0.08, 0.17)
## sig_o.logN (0.05, 0.14)
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
## Denali.PE Denali.95%CI IR.PE IR.95%CI SNF.PE
## a.none 0.37 (0.02, 0.73) 0.17 (-0.02, 0.36) 0.19
## a.Poisson 0.39 (0.02, 0.79) 0.11 (-0.06, 0.31) 0.09
## a.logN 0.41 (0.02, 0.88) 0.14 (-0.05, 0.33) 0.11
## b.none -0.004 (-0.008, 0) -0.008 (-0.016, 0) -0.003
## b.Poisson -0.004 (-0.008, 0) -0.006 (-0.014, 0) -0.001
## b.logN -0.004 (-0.009, 0) -0.007 (-0.015, 0) -0.002
## sig_p.none 0.2 (0.15, 0.28) 0.26 (0.21, 0.32) 0.17
## sig_p.Poisson 0.19 (0.13, 0.27) 0.2 (0.14, 0.27) 0.11
## sig_p.logN 0.19 (0.12, 0.27) 0.24 (0.18, 0.3) 0.13
## sig_o.logN 0.1 (0.03, 0.21) 0.1 (0.03, 0.18) 0.09
## SNF.95%CI
## a.none (-0.02, 0.41)
## a.Poisson (-0.07, 0.28)
## a.logN (-0.06, 0.31)
## b.none (-0.006, 0)
## b.Poisson (-0.004, 0)
## b.logN (-0.005, 0)
## sig_p.none (0.14, 0.22)
## sig_p.Poisson (0.07, 0.16)
## sig_p.logN (0.08, 0.18)
## sig_o.logN (0.04, 0.14)
# 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
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))
## Error: length of 'dimnames' [1] not equal to array extent
# 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)
## V1 V2 V3 V4 V5 V6
## 1 0 (-0.08, 0.09) -0.01 (-0.09, 0.06) 0.01 (-0.05, 0.06)
## 2 0 (-0.08, 0.08) -0.01 (-0.06, 0.04) 0.01 (-0.03, 0.04)
## 3 0 (-0.08, 0.08) -0.01 (-0.08, 0.05) 0.01 (-0.03, 0.04)
## 4 0.22 (0.17, 0.29) 0.27 (0.22, 0.32) 0.18 (0.14, 0.22)
## 5 0.2 (0.14, 0.28) 0.19 (0.13, 0.26) 0.1 (0.07, 0.15)
## 6 0.2 (0.14, 0.29) 0.23 (0.16, 0.3) 0.12 (0.08, 0.17)
## 7 0.09 (0.03, 0.17) 0.11 (0.04, 0.19) 0.09 (0.05, 0.14)
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)
## V1 V2 V3 V4 V5 V6
## 1 0 (-0.08, 0.09) -0.01 (-0.09, 0.06) 0.01 (-0.05, 0.06)
## 2 0 (-0.08, 0.08) -0.01 (-0.06, 0.04) 0.01 (-0.03, 0.04)
## 3 0 (-0.08, 0.08) -0.01 (-0.08, 0.05) 0.01 (-0.03, 0.04)
## 4 0.22 (0.17, 0.29) 0.27 (0.22, 0.32) 0.18 (0.14, 0.22)
## 5 0.2 (0.14, 0.28) 0.19 (0.13, 0.26) 0.1 (0.07, 0.15)
## 6 0.2 (0.14, 0.29) 0.23 (0.16, 0.3) 0.12 (0.08, 0.17)
## 7 0.09 (0.03, 0.17) 0.11 (0.04, 0.19) 0.09 (0.05, 0.14)
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)
## V1 V2 V3 V4 V5 V6
## 1 0.37 (0.02, 0.73) 0.17 (-0.02, 0.36) 0.19 (-0.02, 0.41)
## 2 0.39 (0.02, 0.79) 0.11 (-0.06, 0.31) 0.09 (-0.07, 0.28)
## 3 0.41 (0.02, 0.88) 0.14 (-0.05, 0.33) 0.11 (-0.06, 0.31)
## 4 -0.004 (-0.008, 0) -0.008 (-0.016, 0) -0.003 (-0.006, 0)
## 5 -0.004 (-0.008, 0) -0.006 (-0.014, 0) -0.001 (-0.004, 0)
## 6 -0.004 (-0.009, 0) -0.007 (-0.015, 0) -0.002 (-0.005, 0)
## 7 0.2 (0.15, 0.28) 0.26 (0.21, 0.32) 0.17 (0.14, 0.22)
## 8 0.19 (0.13, 0.27) 0.2 (0.14, 0.27) 0.11 (0.07, 0.16)
## 9 0.19 (0.12, 0.27) 0.24 (0.18, 0.3) 0.13 (0.08, 0.18)
## 10 0.1 (0.03, 0.21) 0.1 (0.03, 0.18) 0.09 (0.04, 0.14)
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)]
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)
## Bayes Freq Bayes Freq Bayes Freq
## a 0.374 0.304 0.167 0.167 0.193 0.194
## a 0.392 0.259 0.113 0.092 0.093 0.068
## a 0.406 0.285 0.141 0.161 0.113 0.088
## b -0.004 -0.003 -0.008 -0.008 -0.003 -0.003
## b -0.004 -0.003 -0.006 -0.005 -0.001 -0.001
## b -0.004 -0.003 -0.007 -0.008 -0.002 -0.001
## sigmap 0.203 0.186 0.259 0.248 0.173 0.163
## sigmap 0.188 0.151 0.197 0.176 0.112 0.094
## sigmap 0.188 0.176 0.236 0.244 0.128 0.109
## Sigobs 0.099 0.050 0.096 0.038 0.090 0.088