title: "originalDMModel" author: "Yoichiro Kanno" date: "August 5, 2014" output: html_document
Shenandoah National Park brook trout data 72 sites in 1996-2010: not all sites were sampled every year & there is a mix of 1-pass and 3-pass surveys ========================================================================
setwd('G:/Conte/Broad spatial modelling/VA/analysis/Dail Madsen 6/ver 6.8/') rm(list=ls()) getwd() library(reshape2) library(rjags) library(plyr) library(ggplot2) library(knitr) library(arm) library(boot) load.module("glm")
load("data/Shen Dail Madsen data and results ver 6.8.RData")
# Data structure nSites = 72; nYears = 15; nAges = 2; nPasses = 3 # Bundle data Dat <- list(nSites=nSites, nYears=nYears, nAges=nAges, y=y, elev=elev, Jdate=Jdate, width=width, fall.flow=fall.flow, winter.flow=winter.flow, spring.flow=spring.flow, spring.temp=spring.temp) # INITIAL VALUES ## Provide first year abundance for S matrix SNew <- array(NA, dim=c(nSites,nYears-1,nAges)) SNew[,1,] <- 100 ## Set intial values Init <- function() list(S=SNew, alpha=2, beta=0.001,lambda=rep(200,nAges), gamma.b=rep(0,5), omega.mean=rep(0.5, nAges), omega.b=array(0, c(5,nAges)), p.mean=rep(0.5,nAges), p.b=rep(0,2))
export.list <- list("Dat", "Init", "SNew", "nSites", "nYears", "nAges", "nPasses", "model", "n.burn") out.yo.coda <- runJagsModel(model = "src/Yoichiro-DM-Model.R", data=Dat, inits=Init, export.list, n.chains = 3, n.burn=1000, params, n.iter = 1000, thin = 3, jags.out = FALSE) CL <- makeCluster(3) clusterExport(cl=CL, export.list, envir = environment()) clusterSetRNGStream(cl=CL, iseed = 2345642) system.time(out <- clusterEvalQ(CL, { library(rjags) load.module('glm') jm <- jags.model(paste0("src/Yoichiro-DM-Model.R"), data=Dat, inits=Init, n.adapt=10000, n.chains=1) fm <- coda.samples(jm, params, n.iter = 20000, thin = 50) return(as.mcmc(fm)) })) out1 <- mcmc.list(out) stopCluster(CL) # Sequential StageBurnin <- jags.model(paste0("src/Yoichiro-DM-Model.R"), Dat, Init, n.chains=3, n.adapt=100) # coda.samples ## Parameters to be monitored ParsStage <- c("omega.mean","omega.b","gamma.b","lambda", "Ntotal","p.mean","p.b","alpha","beta", "N", "p") ## run Niter=20000 out1 <- coda.samples(StageBurnin, ParsStage, n.iter=Niter, thin=100) summary(out1) plot(out1) # jags.samples ## Parameters to be monitored ParsStage <- c("omega.mean","omega.b","gamma.b","lambda", "Ntotal","p.mean","p.b","alpha","beta","N","p") ## run out2 <- jags.samples(StageBurnin, ParsStage, n.iter=Niter, thin=100) str(out2)
## gelman gelman.diag(out1) ## DIC chains=as.matrix(StageBurnin) dev=chains[,"deviance"] dic=mean(dev)+0.5*var(dev) dic.samples(StageBurnin, n.iter=1000, thin=1, type="pD")
library(popbio) # when adult abundance is lower (N = 50 per 100m) N=50 A1 <- matrix(c(0, exp(1.002-0.0075*N), 0.319, 0.453), nrow=2, byrow=TRUE) lambda(A1) A1 sensitivity(A1, zero=TRUE) elasticity(A1) # when adult abundance is higher (N = 100 per 100m) N=100 A2 <- matrix(c(0, exp(1.002-0.0075*N), 0.319, 0.453), nrow=2, byrow=TRUE) lambda(A2) A2 sensitivity(A2, zero=TRUE) elasticity(A2)
################################# ## Graphs below: please ignore ## #################################
biplot <- as.data.frame(cbind(c(out2$omega.mean[1,,1], out2$omega.mean[1,,2], out2$omega.mean[1,,3]), c(out2$omega.mean[2,,1], out2$omega.mean[2,,2], out2$omega.mean[2,,3]), c(out2$omega.b[1,1,,1], out2$omega.b[1,1,,2], out2$omega.b[1,1,,3]), c(out2$omega.b[2,1,,1], out2$omega.b[2,1,,2], out2$omega.b[2,1,,3]), c(out2$omega.b[3,1,,1], out2$omega.b[3,1,,2], out2$omega.b[3,1,,3]), c(out2$omega.b[4,1,,1], out2$omega.b[4,1,,2], out2$omega.b[4,1,,3]), c(out2$omega.b[5,1,,1], out2$omega.b[5,1,,2], out2$omega.b[5,1,,3]), c(out2$omega.b[1,2,,1], out2$omega.b[1,2,,2], out2$omega.b[1,2,,3]), c(out2$omega.b[2,2,,1], out2$omega.b[2,2,,2], out2$omega.b[2,2,,3]), c(out2$omega.b[3,2,,1], out2$omega.b[3,2,,2], out2$omega.b[3,2,,3]), c(out2$omega.b[4,2,,1], out2$omega.b[4,2,,2], out2$omega.b[4,2,,3]), c(out2$omega.b[5,2,,1], out2$omega.b[5,2,,2], out2$omega.b[5,2,,3]), c(out2$alpha[,,1], out2$alpha[,,2], out2$alpha[,,3]), c(out2$beta[,,1], out2$beta[,,2], out2$beta[,,3]), c(out2$gamma.b[1,,1], out2$gamma.b[1,,2], out2$gamma.b[1,,3]), c(out2$gamma.b[2,,1], out2$gamma.b[2,,2], out2$gamma.b[2,,3]), c(out2$gamma.b[3,,1], out2$gamma.b[3,,2], out2$gamma.b[3,,3]), c(out2$gamma.b[4,,1], out2$gamma.b[4,,2], out2$gamma.b[4,,3]), c(out2$gamma.b[5,,1], out2$gamma.b[5,,2], out2$gamma.b[5,,3]), c(out2$p.mean[1,,1], out2$p.mean[1,,2], out2$p.mean[1,,3]), c(out2$p.mean[2,,1], out2$p.mean[2,,2], out2$p.mean[2,,3]), c(out2$p.b[1,,1], out2$p.b[1,,2], out2$p.b[1,,3]), c(out2$p.b[2,,1], out2$p.b[2,,2], out2$p.b[2,,3]))) names(biplot) <- c("omega.mean1","omega.mean2", "omega.b1.1","omega.b2.1","omega.b3.1","omega.b4.1","omega.b5.1", "omega.b1.2","omega.b2.2","omega.b3.2","omega.b4.2","omega.b5.2", "alpha","beta","gamma.b1","gamma.b2","gamma.b3","gamma.b4","gamma.b5", "p.mean1","p.mean2","p.b1","p.b2") ## Pair plot # function for pearson correlation (http://www2.warwick.ac.uk/fac/sci/moac/people/students/peter_cock/r/iris_plots/) panel.pearson <- function(x, y, ...) { horizontal <- (par("usr")[1] + par("usr")[2]) / 2; vertical <- (par("usr")[3] + par("usr")[4]) / 2; text(horizontal, vertical, format(cor(x,y), digits=2)) } pairs(biplot, main="Pairs plot of mcmc estimates at each iteration", upper.panel=panel.pearson) ## Look at omega terms only biplot2 <- as.data.frame(cbind(c(out2$omega.mean[1,,1], out2$omega.mean[1,,2], out2$omega.mean[1,,3]), c(out2$omega.mean[2,,1], out2$omega.mean[2,,2], out2$omega.mean[2,,3]), c(out2$omega.b[1,1,,1], out2$omega.b[1,1,,2], out2$omega.b[1,1,,3]), c(out2$omega.b[2,1,,1], out2$omega.b[2,1,,2], out2$omega.b[2,1,,3]), c(out2$omega.b[3,1,,1], out2$omega.b[3,1,,2], out2$omega.b[3,1,,3]), c(out2$omega.b[4,1,,1], out2$omega.b[4,1,,2], out2$omega.b[4,1,,3]), c(out2$omega.b[5,1,,1], out2$omega.b[5,1,,2], out2$omega.b[5,1,,3]), c(out2$omega.b[1,2,,1], out2$omega.b[1,2,,2], out2$omega.b[1,2,,3]), c(out2$omega.b[2,2,,1], out2$omega.b[2,2,,2], out2$omega.b[2,2,,3]), c(out2$omega.b[3,2,,1], out2$omega.b[3,2,,2], out2$omega.b[3,2,,3]), c(out2$omega.b[4,2,,1], out2$omega.b[4,2,,2], out2$omega.b[4,2,,3]), c(out2$omega.b[5,2,,1], out2$omega.b[5,2,,2], out2$omega.b[5,2,,3]))) names(biplot2) <- c("omega.mean1","omega.mean2", "omega.b1.1","omega.b2.1","omega.b3.1","omega.b4.1","omega.b5.1", "omega.b1.2","omega.b2.2","omega.b3.2","omega.b4.2","omega.b5.2") pairs(biplot2, main="Pairs plot of mcmc estimates at each iteration", upper.panel=panel.pearson)
N.est <- p.est <- array(NA, dim=c(nSites,nYears,nAges)) for(i in 1:nSites){ for(t in 1:nYears){ for(a in 1:nAges){ N.est[i,t,a] <- median(out2$N[i,t,a,1:200,1:3]) p.est[i,t,a] <- median(out2$p[i,t,a,1:200,1:3]) } } } y.est <- array(NA, dim=c(nSites,nYears,nAges,nPasses)) y.est[,,,1] <- N.est*p.est y.est[,,,2] <- N.est*(1-p.est)*p.est y.est[,,,3] <- N.est*(1-p.est)*(1-p.est)*p.est ## YOY, 1st pass yoy1.obs <- adply(y[,,1,1], c(1,2), mean) yoy1.est <- adply(y.est[,,1,1], c(1,2), mean) yoy1.fit <- cbind(yoy1.obs, yoy1.est$V1) names(yoy1.fit) <- c("siteID", "yearID", "yoy1.obs", "yoy1.est") ggplot(yoy1.fit, aes(x=yoy1.obs, y=yoy1.est)) + geom_point() + geom_abline(intercept = 0, slope = 1, colour = "red", size = 1.1, linetype="dashed") + facet_wrap( ~ yearID, scales="free", ncol=5) + xlab("Observed YOY count in 1st pass") + ylab("Predicted YOY count in 1st pass") + labs(title = "Model fit: YOY 1st pass") + theme(axis.title.y = element_text(size = rel(1.5), angle=90), axis.title.x = element_text(size = rel(1.5)), axis.text = element_text(size = rel(1)), plot.title = element_text(size = rel(1.5), colour="blue"), legend.background = element_rect(colour = "black"), legend.text = element_text(size = 15), legend.title = element_blank()) ## YOY, 2nd pass yoy2.obs <- adply(y[,,1,2], c(1,2), mean) yoy2.est <- adply(y.est[,,1,2], c(1,2), mean) yoy2.fit <- cbind(yoy2.obs, yoy2.est$V1) names(yoy2.fit) <- c("siteID", "yearID", "yoy2.obs", "yoy2.est") ggplot(yoy2.fit, aes(x=yoy2.obs, y=yoy2.est)) + geom_point() + geom_abline(intercept = 0, slope = 1, colour = "red", size = 1.1, linetype="dashed") + facet_wrap( ~ yearID, scales="free", ncol=5) + xlab("Observed YOY count in 2nd pass") + ylab("Predicted YOY count in 2nd pass") + labs(title = "Model fit: YOY 2nd pass") + theme(axis.title.y = element_text(size = rel(1.5), angle=90), axis.title.x = element_text(size = rel(1.5)), axis.text = element_text(size = rel(1)), plot.title = element_text(size = rel(1.5), colour="blue"), legend.background = element_rect(colour = "black"), legend.text = element_text(size = 15), legend.title = element_blank()) ## YOY, 3rd pass yoy3.obs <- adply(y[,,1,3], c(1,2), mean) yoy3.est <- adply(y.est[,,1,3], c(1,2), mean) yoy3.fit <- cbind(yoy3.obs, yoy3.est$V1) names(yoy3.fit) <- c("siteID", "yearID", "yoy3.obs", "yoy3.est") ggplot(yoy3.fit, aes(x=yoy3.obs, y=yoy3.est)) + geom_point() + geom_abline(intercept = 0, slope = 1, colour = "red", size = 1.1, linetype="dashed") + facet_wrap( ~ yearID, scales="free", ncol=5) + xlab("Observed YOY count in 3rd pass") + ylab("Predicted YOY count in 3rd pass") + labs(title = "Model fit: YOY 3rd pass") + theme(axis.title.y = element_text(size = rel(1.5), angle=90), axis.title.x = element_text(size = rel(1.5)), axis.text = element_text(size = rel(1)), plot.title = element_text(size = rel(1.5), colour="blue"), legend.background = element_rect(colour = "black"), legend.text = element_text(size = 15), legend.title = element_blank()) ## Adult, 1st pass adult1.obs <- adply(y[,,2,1], c(1,2), mean) adult1.est <- adply(y.est[,,2,1], c(1,2), mean) adult1.fit <- cbind(adult1.obs, adult1.est$V1) names(adult1.fit) <- c("siteID", "yearID", "adult1.obs", "adult1.est") ggplot(adult1.fit, aes(x=adult1.obs, y=adult1.est)) + geom_point() + geom_abline(intercept = 0, slope = 1, colour = "red", size = 1.1, linetype="dashed") + facet_wrap( ~ yearID, scales="free", ncol=5) + xlab("Observed adult count in 1st pass") + ylab("Predicted adult count in 1st pass") + labs(title = "Model fit: adult 1st pass") + theme(axis.title.y = element_text(size = rel(1.5), angle=90), axis.title.x = element_text(size = rel(1.5)), axis.text = element_text(size = rel(1)), plot.title = element_text(size = rel(1.5), colour="blue"), legend.background = element_rect(colour = "black"), legend.text = element_text(size = 15), legend.title = element_blank()) ## Adult, 2nd pass adult2.obs <- adply(y[,,2,2], c(1,2), mean) adult2.est <- adply(y.est[,,2,2], c(1,2), mean) adult2.fit <- cbind(adult2.obs, adult2.est$V1) names(adult2.fit) <- c("siteID", "yearID", "adult2.obs", "adult2.est") ggplot(adult2.fit, aes(x=adult2.obs, y=adult2.est)) + geom_point() + geom_abline(intercept = 0, slope = 1, colour = "red", size = 1.1, linetype="dashed") + facet_wrap( ~ yearID, scales="free", ncol=5) + xlab("Observed adult count in 2nd pass") + ylab("Predicted adult count in 2nd pass") + labs(title = "Model fit: adult 2nd pass") + theme(axis.title.y = element_text(size = rel(1.5), angle=90), axis.title.x = element_text(size = rel(1.5)), axis.text = element_text(size = rel(1)), plot.title = element_text(size = rel(1.5), colour="blue"), legend.background = element_rect(colour = "black"), legend.text = element_text(size = 15), legend.title = element_blank()) ## Adult, 3rd pass adult3.obs <- adply(y[,,2,3], c(1,2), mean) adult3.est <- adply(y.est[,,2,3], c(1,2), mean) adult3.fit <- cbind(adult3.obs, adult3.est$V1) names(adult3.fit) <- c("siteID", "yearID", "adult3.obs", "adult3.est") ggplot(adult3.fit, aes(x=adult3.obs, y=adult3.est)) + geom_point() + geom_abline(intercept = 0, slope = 1, colour = "red", size = 1.1, linetype="dashed") + facet_wrap( ~ yearID, scales="free", ncol=5) + xlab("Observed adult count in 3rd pass") + ylab("Predicted adult count in 3rd pass") + labs(title = "Model fit: adult 3rd pass") + theme(axis.title.y = element_text(size = rel(1.5), angle=90), axis.title.x = element_text(size = rel(1.5)), axis.text = element_text(size = rel(1)), plot.title = element_text(size = rel(1.5), colour="blue"), legend.background = element_rect(colour = "black"), legend.text = element_text(size = 15), legend.title = element_blank())
yearTotalAbundMean <- adply(out2$Ntotal, c(1,2), mean) yearTotalAbundCI <- adply(out2$Ntotal, c(1,2), quantile, probs=c(.025, .975)) year2 <- rep(1996:2010, 2) size2 <- rep(c("YOY","adult"), each=15) yearTotalAbund <- cbind(yearTotalAbundMean, yearTotalAbundCI[,3:4], year2, size2) names(yearTotalAbund) <- c("year","size","total.mean","total.lower","total.upper","year2","size2") yearTotalAbund$year2 <- as.numeric(as.character(yearTotalAbund$year2)) ggplot(yearTotalAbund, aes(x=year2, y=total.mean, colour=size2)) + geom_errorbar(aes(ymin=total.lower, ymax=total.upper), width=0.5, size=1.2) + geom_line(size=1.2) + geom_point(size=1.2) + xlab("Year") + ylab("Estimated total abundance") + labs(title = "Yearly total abundance (with 95 % CI) over 72 sites in Shenandoah") + theme(axis.title.y = element_text(size = rel(1.5), angle=90), axis.title.x = element_text(size = rel(1.5)), axis.text = element_text(size = rel(1.2)), plot.title = element_text(size = rel(1.5), colour="blue"), legend.background = element_rect(colour = "black"), legend.text = element_text(size = 15), legend.title = element_blank())
yearTotalAbundMean <- adply(out2$Ntotal, c(1,2), mean) yearTotalAbundCI <- adply(out2$Ntotal, c(1,2), quantile, probs=c(.025, .975)) year2 <- rep(1996:2010, 2) size2 <- rep(c("YOY","adult"), each=15) yearTotalAbund <- cbind(yearTotalAbundMean, yearTotalAbundCI[,3:4], year2, size2) names(yearTotalAbund) <- c("year","size","total.mean","total.lower","total.upper","year2","size2") yearTotalAbund$year2 <- as.numeric(as.character(yearTotalAbund$year2)) ggTot <- ggplot(yearTotalAbund, aes(x=year2, y=total.mean, colour=size2)) + geom_errorbar(aes(ymin=total.lower, ymax=total.upper), width=0.5, size=1.2) + geom_line(size=1.2) + geom_point(size=1.2) + scale_colour_manual(values = c("#003300","brown")) + xlab("Year") + ylab("Estimated total abundance") + #labs(title = "Yearly total abundance (with 95 % CI) over 72 sites in Shenandoah") + theme_bw() + theme(panel.border=element_rect(colour='black'), panel.background=element_rect(fill='#FDF5B7', colour='black'), panel.grid.major=element_line(colour=NA), panel.grid.minor=element_line(colour=NA), plot.background=element_rect(fill="#FDF5B7"), axis.title.y = element_text(size = rel(1.8), angle=90), axis.title.x = element_text(size = rel(1.8)), axis.text = element_text(size = rel(1.5)), plot.title = element_text(size = rel(1.8), colour="black", face="bold"), legend.position="none") #legend.background = element_rect(colour = "black"), #legend.text = element_text(size = 15), #legend.title = element_blank()) ggTot dpiIn <- 600 ggsave( file=paste(getwd(),'/ShenTotalBen.png',sep=''), plot=ggTot, dpi=dpiIn , width=8, height=5, units='in',scale=2 )
## make df lag.df <- adply(out2$Ntotal, c(1,2), mean) names(lag.df) <- c("year","age","totalN") lag.df.wide <- dcast(lag.df, year ~ age, value.var="totalN") names(lag.df.wide) <- c("year","yoyTotal","adultTotal") ## adult at year t vs. yoy at year t+1 lag.df.wide2 <- as.data.frame(cbind(lag.df.wide$adultTotal[1:(nYears-1)], lag.df.wide$yoyTotal[2:nYears])) library(splines);library(MASS) ggplot(lag.df.wide2, aes(x=V1,y=V2)) + geom_point(size=4) + stat_smooth(method="lm", formula= y ~ ns(x,2), se=FALSE, size=1.2, linetype="dashed") + xlab("Total adult abundance at year t") + ylab("Total YOY abundance at year t+1") + labs(title = "Based on analysis of 72 sites in 1996-2010 in Shenandoah") + theme_bw() + theme(axis.title.y = element_text(size = rel(1.5), angle=90), axis.title.x = element_text(size = rel(1.5)), axis.text = element_text(size = rel(1.2)), plot.title = element_text(size = rel(1.5), colour="blue"), panel.border=element_rect(colour='black'), panel.background=element_rect(fill='white', colour='black'), panel.grid.major=element_line(colour=NA), panel.grid.minor=element_line(colour=NA)) ## yoy at year t vs. adult at year t+1 lag.df.wide3 <- as.data.frame(cbind(lag.df.wide$yoyTotal[1:(nYears-1)], lag.df.wide$adultTotal[2:nYears])) ggplot(lag.df.wide3, aes(x=V1,y=V2)) + geom_point(size=4) + stat_smooth(method="lm", se=FALSE, size=1.2, linetype="dashed") + xlab("Total YOY abundance at year t") + ylab("Total adult abundance at year t+1") + labs(title = "Based on analysis of 72 sites in 1996-2010 in Shenandoah") + theme_bw() + theme(axis.title.y = element_text(size = rel(1.5), angle=90), axis.title.x = element_text(size = rel(1.5)), axis.text = element_text(size = rel(1.2)), plot.title = element_text(size = rel(1.5), colour="blue"), panel.border=element_rect(colour='black'), panel.background=element_rect(fill='white', colour='black'), panel.grid.major=element_line(colour=NA), panel.grid.minor=element_line(colour=NA))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.