\documentclass[11pt]{article} \usepackage[american]{babel} \usepackage[utf8]{inputenc} \newcommand{\R}{{\sf R}} \newcommand{\Splus}{{\sf S-PLUS}} \newcommand{\fixme}[1]{\textbf{FIXME: #1}} %\newcommand{\fixme}[1]{} \usepackage{url} \usepackage{alltt} \usepackage{fancyvrb} % with VerbatimFootnotes, allow verb in footnotes \usepackage{listings} \usepackage{verbatim} \usepackage{hyperref} \usepackage{natbib} \usepackage{listings} \usepackage{graphicx} \usepackage{textcomp} %had to add this due to correct 'LaTeX Font Warning: Font shape `TS1/aett/m/sl' undefined" \DeclareGraphicsExtensions{.pdf} \newcommand{\code}[1]{{\tt #1}} \bibliographystyle{chicago} \title{Comparison of ADMB-RE and JAGS \\ Bayesian State-space length disaggregated population model: Estimating total mortality by decade for winter skate (\textit{Leucoraja ocellata})} \author{Trevor Davies, Dalhousie University, Halifax, Canada \\ Steven J.D. Martell, International Pacific Halibut Commission, Seattle WA} \date{\today} \begin{document} \maketitle \SweaveOpts{keep.source=TRUE} %Below are some of the basic R-packages we use to run the models and retrieve results into the R environment. <<>>= library(reshape2) ## for 'melt' function with ggplot2 library(ggplot2) ## pictures library(coda) ## analysis of MCMC runs (JAGS \& ADMB) library(R2admb) ## R interface to AD Model Builder library(R2jags) ## R interface to JAGS (BUGS dialect) library(xtable) ## Makes latex tables from R dataframes \& objects library(gdata) ## Rename vars package @ % Now load in the workspaces for the relevant objects <>= options(continue=" ") rm(list=ls()) load("../ADMB/tfit_admb.RData") # read in admb model object load("../ADMB/fit.RData") # read in admb summarized results results_admb <- results load("../ADMBmcmc/tfit_admbmcmc.RData") # read in admb mcmc model object load("../ADMBmcmc/fit.RData") # read in admb mcmc summarized results results_admbmcmc <- results load('../BUGS/tfit_jags.RData') # read in jags model object load('../BUGS/fit.RData') # read in jags summarized results results_jags <- results data1 <- read.csv('../DATA/skate_bugs.dat',sep='\t',header=T) # these are the raw data @ \section{Summary} The model described here is a fully Bayesian state-space model implemented in both JAGS \citep{plummer03jags} and AD Model Builder (ADMB; \citep{Fournieretal12}) with random effects. The goal of the model was to obtain decadal mortality estimates of three different size classes of winter skates (\textit{Leucoraja ocellata}) on the eastern Scotian Shelf. The three time series of length disaggregated stages (Figure \ref{thedata}; Juvenile1, Juvenile2, and adult) are largely non-informative to the model parameters: catchability, recruitment rate, and stage transition probability. Consequently, it was necessary to include highly informative priors for these parameters to obtain reasonable mortality estimates. We did not implement this model in R due to the model containing 28 hyper parameters and 105 random effects which would have resulted in excessively long run times. The full model description and alternative model formulations are fully described in \cite{swain:09:1347}. \begin{itemize} \item ADMB-RE, ADMB in MCMC mode, and JAGS all gave similar results \item ADMB-RE was able to obtain maximum likelihood estimate (MLEs) much more; quickly than JAGS was able to converge to a stationary; \item ADMB-RE run in MCMC was unable to converge to stationary; posteriors even when twice the burn-in of the JAGS model was implemented \item Chain mixing appears to be much better in JAGS compared to ADMB-RE; however, the effect of this on parameter estimates is unknown. \item Coding of JAGS model was easier compared to ADMB-re \end{itemize} \begin{center} \begin{tiny} <>= time1 <- data.frame(ADMB=round(as.numeric(results_admb$time),1),ADMBmcmc=round(as.numeric(results_admbmcmc$time),1),JAGS=round(as.numeric(results_jags$time[3]),1)) cap <- 'Run times of ADMB-RE (Maximum Likelihood Estimate only; MLE); ADMB-RE with MCMC mode on (MCMC), and JAGS run times in seconds. Both JAGS and ADMB-RE were run using a single chain (a requirement of ADMB-RE) for 175 000 iterations, discarding the first 100 000 for burn-in, and sampling every 30th sample to reduce autocorrelation between samples.' label1 <- 'table:runtimes' table.x <- xtable(time1,caption=cap, label=label1) digs <- matrix(rep(1, nrow(time1) * (ncol(time1))) ,nrow=nrow(time1),ncol=ncol(time1)) digs <- cbind(0,(digs)) # adds column to cover parameter names digits(table.x) <- digs print(table.x,type='latex',include.rownames=FALSE,sanitize.text.function = function(x){x},caption.placement="top",table.placement='!htp') # USe long table because @ \end{tiny} \end{center} \section{Introduction} The model implemented here is a stage-structured state-space model that was previously used to explore trends in natural mortality for winter skate (\textit{Leucoraja ocellata}) in the eastern Scotian Shelf. Data on relative abundance for each of three length-classes (or stages) were obtained from bottom trawl surveys in July and September since 1970 (Figure \ref{thedata}). Age-disaggregated data are not available for this species. \section{The Data} Three time series of relative numbers of individuals by size class are as follows: <>= par(mfrow=c(3,1),mar=c(0,4,0,4),oma=c(3.5,1.5,3.5,1.5)) with(data1,plot(year,Juv1,xlab="",ylab="",ylim=c(0,max(Juv1)),xlim=c(min(year),max(year)),mgp=c(1.25,.25,.25),tcl=-0.2,axes=FALSE,type='b')) # box() axis(2,tick = TRUE, line =0,tcl=-0.2,las=1,cex.axis=1,mgp=c(1.25,.25,1)) mtext('Juvenile Stage 1',side=2,line=4) with(data1,plot(year,Juv2,xlab="",ylab="",ylim=c(0,max(Juv2)),xlim=c(min(year),max(year)),mgp=c(1.25,.25,.25),tcl=-0.2,axes=FALSE,type='b')) # box() axis(2,tick = TRUE, line =0,tcl=-0.2,las=1,cex.axis=1,mgp=c(1.25,.25,1)) mtext('Juvenile Stage 2',side=2,line=4) mtext('(millions of skates)',side=2,line=2.25,cex=1.25) with(data1,plot(year,Adult,xlab="",ylab="",ylim=c(0,max(Adult)),xlim=c(min(year),max(year)),mgp=c(1.25,.25,.25),tcl=-0.2,axes=FALSE,type='b')) # box() axis(2,tick = TRUE, line =0,tcl=-0.2,las=1,cex.axis=1,mgp=c(1.25,.25,1)) axis(1,tick = TRUE, line =0,tcl=-0.2,las=0,cex.axis=1,mgp=c(1.25,.25,1)) mtext('Adult Stage',side=2,line=4) mtext('Year',side=1,line=2.5) @ \begin{figure}[!htbp] \centering \includegraphics[width=5in]{skate-thedata.pdf} \caption{Survey abundance indices of three length classes of winter skate on the easter Scotian Shelf. Indices were scaled up to absolute abundance by catchability coefficients that were estimated outside the model. Note different scales on the y-axis. Data from Swain et al. (2009).} \label{thedata} \end{figure} \section{The Models} The model described here is a fully Bayesian state-space model and was implemented in both AD Model Builder with random effects (ADMB-RE) \citep{Fournieretal12} and JAGS \citep{plummer03jags} The population dynamics model for this species is divided into three stages, where the transition from small to larger stages occurs via a transition probability ($\theta_i$), where $i$ is the index for the stage. The state transition equation for the first stage is given by eq. \ref{model:juv1}, where $N_{1,t}$ is the number of individuals in stage 1 in year $t$, $\theta_1$ is the transition probability from stage 1 to stage 2, $r$ is the annual recruitment rate, $N_{3,t-a}$ is the number of stage-3 individuals in year lagged by $a$ years to represent the time required to recruit to the first stage. The total mortality rate $Z_{1,d}$ by decade $d$ which is the main parameter of interest. \begin{equation} N_{1,t} = \left[N_{1,t-1}(1-\theta_1) + \frac{1}{2} (rN_{3,t-a)}\right]e^{-Z_{1,d}}e^{\eta_{1,t}} \label{model:juv1} \end{equation} The state transition equations for stage two and three are given in eqs \ref{model:juv2} and \ref{model:adult}. \begin{equation} N_{2,t} = \left[N_{2,t-1}(1-\theta_2) + N_{1,t-1} \theta_1 \right] e^{-Z_{2,d}}e^{\eta_{2,t}} \label{model:juv2} \end{equation} \begin{equation} N_{3,t} = \left[N_{3,t-1} + N_{2,t-1} \theta_2 \right] e^{-Z_{3,d}}e^{\eta_{3,t}} \label{model:adult} \end{equation} \noindent where $N_{i,t}$ is the skate abundance (millions) for stage \textit{i} in year \textit{t}, $\theta_{1:2}$ is the transition probability between the first to second juvenile stage and second juvenile and adult age respectively, \textit{r} is the recruitment rate, \textit{a} accounts for the 5 year lag from the adult stage to entering the first juvenile stage, \textit{Z} is the instantaneous mortality rate for each stage and each decadal period, \textit{d}, and $\eta_{i,t}$ is an independent normal random variable for each stage and year with a mean of zero and process variance $\sigma^2_{i}$ to account for stochasticity in the population dynamics for each stage, $\sim LN(0, \sigma^2_i)$. All parameters are held constant with the exception of $Z_i$ which is allowed to vary decadally. The observation model relates the unobserved states, $N_{i,t}$ of the population model to the indices of abundance $y_{i,t}$, that are observed with error. \begin{equation} y_{i,t} = q_{i} N_{i,t}e^{\epsilon_{i,t}} \label{model:obs} \end{equation} \noindent where $y_{i,t}$ is the scaled index of relative abundance of stage \textit{i} at time \textit{t}; $q_i$ is the catchability coefficient for each stage \textit{i} which describes the effectiveness of each unit of fishing effort; and $\epsilon_{i,t}$ is a normal random variable with a mean of zero and variance $\tau_{i}^2$ to account for error in the observations of index \textit{i}, $\epsilon_{i,t} \sim LN(0, \tau_{i}^2)$. The time series data is largely uninformatitve to many of the parameters in the model, therefore these data are supplied in the form of priors. The bayesian approach requires priors be used to supply additional external information to the model. Flat priors were used when prior information was not known about certain model parameters. Priors used in these analyses are summarized in Table \ref{table1:priors}. \begin{table}[!htp] \caption{Summary of specified priors for Bayesian state-space model} \begin{center} \begin{tabular}{ll} \hline Parameter & Prior \\ \hline $recrate$ (Recruitment rate) & log normal(1.18, 4)I(0.01,1000)\\ $theta.1$ (Transition probability 1) & beta(20, 80) \\ $theta.2$ (Transition probability 2) & beta(25, 75)\\ All $z$ (Total mortality) & uniform(0, 3) \\ $q$ (Catchability coefficient) & uniform(0,100) \\ \hline \end{tabular} \end{center} \label{table1:priors} \end{table} \section{AD Model Builder with random effects} ADMB has the benefit of being able to quickly obtain MLE of the parameters of interest compared to JAGS which must obtain the full stationary distribution profiles of all model parameters. This can result in a large time saving advantage of using ADMB over JAGS when analyzing large dataset and implimenting complex models. This feature of ADMB can be particularly useful during the model formulation and debugging phases of model development and full MCMC profiles can then be obtained once the model has reached a satisfactory level of development. Although this advantage is not a serious issue with the model explored here, this flexibility can potentially offset the generally larger time investment of formulating complex models in ADMB over those coded in JAGS. We ran the ADMB -RE model in both MLE and MCMC modes. To implement MCMC in ADMB you use the arguments -mcmc2 \textit{N} -mcsave \textit{N2} commands to run the model with a chain of length \textit{N} and save the output every \textit{N2} steps. We use the argument -mcmc2 to generate a chain on the joint vector of hyper-parameters and random effects. A key difference between JAGS and ADMB is that ADMB does not begin sampling the MCMC chains until a maximum likelihood estimate has been identified and therefore presumably a shorter burn-in period is required. The length of burn-in, however, is dependent upon model complexity and identifiability of parameter estimates. We therefore, ran ADMB for the same number of iterations and used the same burn-in as the JAGS model. \newpage \subsection{ADMB-RE Results: Maximum Likelihood estimates} The ADMB-RE MLE model obtained parameter estimates the most quickly and only required \Sexpr{with(time1,ADMB)} seconds to run. \begin{center} \begin{tiny} <>= keepers.name <-c("recrate", "theta.1", "theta.2", "z.1", "z.2", "z.3", "z.4", "z.5", "z.6", "z.7", "z.8", "z.9") admbmle.sum1 <- coef(summary(tfit_admb))[keepers.name,] row.names(admbmle.sum1) <- c("recrate","theta.1", "theta.2", "z.juv1.1", "z.juv1.2" , "z.juv1.3" ,"z.juv2.1", "z.juv2.2", "z.juv2.3", "z.adult.1", "z.adult.2", "z.adult.3") cap <- 'Maximum likelihood estimates (MLEs) and associated standard errors of the recruitment rate (recrate), transition probability between size stages (theta 1 and 2) and mortality estimates for all stages and time periods. Subscripts to mortality parameter estimates ($Z_{s,d}$), are subscripted by stage and decade' label1 <- 'table:admbMLE' print(xtable(admbmle.sum1,caption=cap, label=label1), type="latex",table.placement="!htp") @ \end{tiny} \end{center} \newpage \subsection{ADMB Results: MCMC} The ADMB MCMC model took substantially longer than the ADMB-RE MLE model and finished in \Sexpr{with(time1,ADMBmcmc)} seconds. Because ADMB only runs one MCMC chain we used Geweke's diagnostic to assess the convergence to a stationary distribution of the key model parameters (Figs \ref{admb:geweke1} \& \ref{admb:geweke2}). The existence of some Geweke z-scores outside 2 standard deviations of zero suggests that there may be some convergence issues with some of the parameters and the model may require a longer burn-in period to achieve satisfactory model convergence. The Geweke's convergence diagnostic is available in the R coda package \citep{R:coda} and convergence failure is identified by a large number of Z-scores falling outside the 95\% confidence intervals. \begin{figure}[!htbp] \centering <>= sam.keeps <- 2500 admb.mcmc1 <- tail(tfit_admbmcmc$mcmc,sam.keeps) # allows same burn in as JAGS admb.mcmc2 <- subset(admb.mcmc1,select=c('recrate','theta.1','theta.2')) geweke.plot(as.mcmc(admb.mcmc2)) @ \caption{Geweke diagnostoc of select model paramters of ADMB-RE MCMC. Dashed lines identify 2 standard deviation boundaries around zero.} \label{admb:geweke1} \end{figure} \begin{figure}[!htbp] \centering <>= old.zs <- c(paste('z',seq(1,9,3),sep='.'),paste('z',seq(2,9,3),sep='.'),paste('z',seq(3,9,3),sep='.')) new.zs <- c("z.1.1","z.1.2", "z.1.3","z.2.1","z.2.2", "z.2.3","z.3.1","z.3.2","z.3.3") all.parms <- c("recrate","theta.1","theta.2",old.zs) all.parms2 <- c("recrate","theta.1","theta.2",new.zs) admb.mcmc2 <- subset(admb.mcmc1,select=all.parms) admb.mcmc3 <- rename.vars(admb.mcmc2,from=(old.zs),to=new.zs, info = FALSE) geweke.plot(as.mcmc(subset(admb.mcmc3,select=new.zs))) @ \caption{Geweke diagnostic of select model paramaters of ADMB-RE MCMC. Dashed lines identify 2 standard deviation boundaries around zero.} \label{admb:geweke2} \end{figure} We examined trace plots of samples of select model parameters as another measure to assess model convergence (Fig. \ref{admbmcmc:trace}). Examination of the trace plots suggests that there is substantial autocorrelation between samples even though we specified a thinning for every 30 samples. This suggests that further burn-in and thinning may be necessary to achieve model convergence. We did, however, increase burn-in to 225 000 iterations and substantial autocorrelation was still present in the trace plots. \begin{figure}[!htbp] \centering <>= admb.mcmc4 <- subset(admb.mcmc1,select=c("recrate","theta.1","theta.2")) xyplot(as.mcmc(admb.mcmc4)) @ \caption{Traceplots of select ADMB-RE mcmc model parameters} \label{admbmcmc:trace} \end{figure} % Below is the latex table results \begin{center} \begin{tiny} <>= table1 <-t(apply(admb.mcmc3,2,quantile,c(0.5,0.025,0.975))) table2 <- data.frame(Estimate=tfit_admbmcmc$coefficients[all.parms],Std.Error=tfit_admbmcmc$se[all.parms]) table3 <- cbind(table2,table1) rownames(table3) <- rownames(table1) cap <- 'Parameter estimates of select parameters, associated standard errors, and quantiles of ADMB MCMC results. Subscripts to total mortality parameter estimates estimates ($Z_{s,d}$), are subscripted by stage and decade' label1 <- 'table:admbmcmc' print(xtable(table3,caption=cap, label=label1), type="latex",table.placement="!htbp") @ \end{tiny} \end{center} \newpage \section{BUGS (using JAGS)} The JAGS model required \Sexpr{with(time1,JAGS)} seconds to run. There are a variety of methods to gauge convergence of an MCMC chain \citep{cowles:96:883}. Although JAGS can run multiple chains which allows for a variety of convergence diagnostics, ADMB can only run single chains and we therefore also ran the JAGS model with a single chain to maintain consistency between the two platforms. We also used the Geweke's convergence diagnostic and trace plots to gauge convergence. Examination of the Geweke's diagnostic plots shows much better convergence of the JAGS model parameters (Figures \ref{jags:geweke1} \& ref{jags:geweke2} compared to the ADMB-RE model for the majority of parameters examined. Further, trace plots \ref{jags:trace} shows good mixing of the select model parameters which is indicative that a stationary posterior distribution had been reached. \begin{figure}[!htbp] \centering <>= jags.df <- as.data.frame(tfit_jags_mcmc) jags.parms1 <- c("deviance", "recrate", "theta[1]", "theta[2]") sub.jags.df1 <- subset(jags.df,select=jags.parms1) geweke.plot(as.mcmc(sub.jags.df1)) @ \caption{Geweke plots of select model parameters for model run in JAGS. Dashed lines are 95\% confidence regions.} \label{jags:geweke1} \end{figure} \begin{figure}[htp!] \centering <>= jags.zs1 <- c("z[1,1]","z[1,2]", "z[1,3]","z[2,1]", "z[2,2]", "z[2,3]", "z[3,1]", "z[3,2]", "z[3,3]") sub.jags.df2 <- subset(jags.df,select=jags.zs1) sub.jags.df2 <- rename.vars(sub.jags.df2,from=jags.zs1,to=new.zs, info = FALSE) geweke.plot(as.mcmc(sub.jags.df2)) @ \caption{Geweke plots of decadal mortality estimates [stage,decade] for model run in JAGS. Dashed lines are 95\% confidence regions. } \label{jags:geweke2} \end{figure} % Below are the Jags results. \begin{figure}[htp!] \centering <>= xyplot(as.mcmc(sub.jags.df1)) @ \caption{Traceplots of select JAGS model parameters after burning of 100000 iterations (single chain)} \label{jags:trace} \end{figure} \newpage \subsection{JAGS Results} Summary of JAGS results are as follows: \begin{center} \begin{tiny} <>= temp1.sd <- ((summary(tfit_jags_mcmc)[[1]])) temp2.sd <- temp1.sd[c(jags.parms1,jags.zs1),] jags.new.rows <- c("deviance", "recrate", "theta.1", "theta.2", "z.1.1","z.2.1", "z.3.1","z.1.2","z.2.2", "z.3.2","z.1.3","z.2.3","z.3.3") row.names(temp2.sd) <- jags.new.rows #Quants temp1.q <- ((summary(tfit_jags_mcmc)[[2]])) temp2.q <- temp1.q[c(jags.parms1,jags.zs1),] row.names(temp2.q) <- jags.new.rows temp3.q <- subset(as.data.frame(temp2.q),select=c("50%","2.5%","97.5%")) all.jags <- round(cbind(temp2.sd,temp3.q),2) cap <- 'Select Parameter and associated standard deviation from JAGS. Subscripts to total mortality parameter estimates estimates ($Z_{s,d}$), are subscripted by stage and decade' label1 <- 'table:jagsmcmc1' print(xtable(all.jags,caption=cap, label=label1), type="latex",table.placement="!htp") @ \end{tiny} \end{center} \newpage \section{Comparison of ADMB-RE MCMC and JAGS bayesian posteriors} Posterior distributions of select model parameters of the ADMB-RE MCMC and JAGS models are compared in Fig \ref{jags_admb_compare}. Posterior distributions between the two modelling platforms were similar, however, the suspect convergence diagnostics of the ADMB-RE MCMC model may be an explanation for the descrepencies between the two approaches. <>= admb.params <- c("recrate", "theta.1", "theta.2","z.1", "z.2", "z.3", "z.4", "z.5", "z.6", "z.7", "z.8", "z.9") admb.mcmc <- tfit_admbmcmc$mcmc[,admb.params] jags.params <- c("recrate", "theta[1]", "theta[2]", "z[1,1]", "z[2,1]", "z[3,1]","z[1,2]", "z[2,2]", "z[3,2]", "z[1,3]", "z[2,3]", "z[3,3]") jags.mcmc <- tfit_jags_mcmc[,jags.params] new.zs <- c("z.1.1","z.1.2", "z.1.3","z.2.1","z.2.2", "z.2.3","z.3.1","z.3.2","z.3.3") new.parms <- c("recrate","theta.1","theta.2") admb.mcmc <-rename.vars(admb.mcmc,from=admb.params,to=c(new.parms,new.zs),info = FALSE) colnames(jags.mcmc) <-c(new.parms,new.zs) jags.wide <- data.frame(melt(as.data.frame(jags.mcmc)),model='JAGS') admb.wide <- data.frame(melt(admb.mcmc),model='ADMBmcmc') ##Priors - need to do the recrate this way to deal with facet scaling sim.sum1 <- nrow(subset(admb.wide,variable=='recrate')) recrate.prior1 <- rlnorm(n=(2*sim.sum1),meanlog=1.18,sdlog=log(4)) recrate.prior2 <- recrate.prior1[recrate.prior1 <=7] recrate.prior3 <- sample(recrate.prior2,sim.sum1) recrate.prior <- data.frame(variable='recrate',value=recrate.prior3,model='prior') theta1.prior <- data.frame(variable='theta.1',value=rbeta(n=sim.sum1,20,80),model='prior') theta2.prior <- data.frame(variable='theta.2',value=rbeta(n=sim.sum1,25,75),model='prior') z.priors <- data.frame(variable=rep(new.zs,each=10*sim.sum1),value=runif(9*10*sim.sum1,min=0,max=3),model='prior') prior.wide <- rbind(recrate.prior,theta1.prior,theta2.prior,z.priors) all.density <- rbind(jags.wide,admb.wide,prior.wide) ggplot(all.density, aes(x=value,fill=model)) + geom_density(alpha=.3) + facet_wrap(~variable, ncol=4,scales="free") + theme_bw() + ylab("Density") +theme(axis.text.x = element_text(size =6, colour = 'black')) @ \begin{figure}[!htp] \centering \includegraphics[width=5in]{skate-jags_admb_compare.pdf} \caption{Comparison of posterior and prior densities of select model parameters from ADMB-RE Markov Chain Monte Carlo (MCMC) and JAGS analysis.} \label{jags_admb_compare} \end{figure} Time series predictions of the three modelling approaches indicate all three give similar estimates on population abundance through time (Fig. \ref{time_series_compare}). <>= # Observation darta data.juv1 <- with(data1, data.frame(year=year,Juv1=log(Juv1))) data.juv2 <- with(data1, data.frame(year=year,Juv2=log(Juv2))) data.adult <- with(data1, data.frame(year=year,Adult=log(Adult))) ### JAGS # Take conf of everything jags.mcmc1 <- apply(tfit_jags_mcmc,2,quantile,c(0.5,0.025,0.975)) jags.juv1.names1 <- paste("logIpred[1,",seq(1,35,1),']',sep='') jags.juv2.names1 <- paste("logIpred[2,",seq(1,35,1),']',sep='') jags.adult.names1 <- paste("logIpred[3,",seq(1,35,1),']',sep='') jags.names1 <- c(jags.juv1.names1,jags.juv2.names1,jags.adult.names1) jags.wide <- data.frame(year=rep(seq(1970,2004,1),3),model='JAGS',t(jags.mcmc1[,jags.names1]),stage=rep(c('Juv1','Juv2','Adult'),each=35)) ### MLE ADMB - JUV1 admb.mle.conf <- as.matrix(results_admb$confint.quad) admb.mle.conf2 <- admb.mle.conf[rownames(admb.mle.conf)=='Ipred',] admb.wide <- data.frame(year=rep(1970:2004,3),model='ADMBmle',median=tfit_admb$coefficients[names(tfit_admb$coefficients)=='Ipred'],quant.025=admb.mle.conf2[,1],quant.975=admb.mle.conf2[,2],stage=rep(c('Juv1','Juv2','Adult'),each=35)) names(jags.wide) <- names(admb.wide) #fix jags names ### MCMC ADMB Ipredmcmc <- read.csv('../ADMBmcmc/Ipred_posteriors.dat',header=FALSE,sep=' ')[,-1] # Read one line colnames(Ipredmcmc) <- paste('X',seq(1970,2004,1)) stage2 <- (stage=rep(rep(c('Juv1','Juv2','Adult'),2),nrow(Ipredmcmc/3))) Ipredmcmc2 <- data.frame(stage=stage2,Ipredmcmc) Juv1admbmcmc <- data.frame(year=seq(1970,2004,1),model='ADMBmcmc',as.data.frame(t(apply(tail(subset(Ipredmcmc2, stage=='Juv1')[,-1],2500),2,quantile,c(0.5,0.025,0.975)))),stage='Juv1') Juv2admbmcmc <- data.frame(year=seq(1970,2004,1),model='ADMBmcmc',as.data.frame(t(apply(tail(subset(Ipredmcmc2, stage=='Juv2')[,-1],2500),2,quantile,c(0.5,0.025,0.975)))),stage='Juv2') adultadmbmcmc <- data.frame(year=seq(1970,2004,1),model='ADMBmcmc',as.data.frame(t(apply(tail(subset(Ipredmcmc2, stage=='Adult')[,-1],2500),2,quantile,c(0.5,0.025,0.975)))),stage='Adult') admb.mcmc.wide <- rbind(Juv1admbmcmc,Juv2admbmcmc,adultadmbmcmc) colnames(admb.mcmc.wide) <- colnames(jags.wide) # Put all data together a <- rbind(jags.wide,admb.wide,admb.mcmc.wide) data2 <-data.frame(year=rep(1970:2004,3),melt(data1[,2:4])) a2 <- data.frame(a,data1=data2[,3]) # Change plot order a2 <- within(a2, stage <- factor(stage, levels = c('Juv1','Juv2','Adult'))) #### GGPLOT colors1 <- c('red','yellow','blue') ggplot(data=a2, aes(year, median,colour=model,fill=model)) + geom_ribbon(aes(ymin = quant.025, ymax = quant.975), colour=NA, alpha = .25) + facet_grid(stage ~. , scales="free_y") + geom_point(aes(x=year,y=log(data1)),colour='black') + geom_line(size=1.25,alpha=1) + theme_bw() + ylab("Log abundance / catchability") + guides(colour=guide_legend(title="model")) + scale_colour_manual(values = colors1)+ scale_fill_manual(values = colors1) @ \begin{figure}[htp!] \centering \includegraphics[width=6in]{skate-timeseries_compare.pdf} \caption{Comparison of the three model fits by JAGS, ADMB-RE maximum likelihood and ADMB-RE MCMC. Fit lines and 95\% confidence limits are estimated abundance for each stage divided by catchability (logged). } \label{time_series_compare} \end{figure} \newpage \section{Simulation results} <>= load("../SIMS/allsim.RData") source("../../TOOLS/compile-funs.R") source("../../TOOLS/sim_plotfuns.R") true.data <- read.table("../SIMS/true.dat",header=TRUE,as.is=TRUE) @ We used a simulation approach to assess the ability of each modelling platform to accurately and precisely estimate model parameters. We simulated \Sexpr{length(sims)} datasets using fixed parameter values with error added in the form of both process and observation error (Table \ref{table:skatesim}). Data was simulated in \textbf{R} and identical data sets were given to ADMB and JAGS for each simulation. ADMB-RE was run in MLE mode only for these tests. JAGS appeared to be marginally better at estimating the mortality parameters in comparison to ADMB (Figure \ref{simsplot1}). The JAGS model, however, appeared to have a large bias in both transition probability parameters and consistently overestimated recruitment rate (Fig. \ref{simsplot2}). \begin{figure}[!htbp] \centering <>= new.names <- c("z.juv1.1", "z.juv1.2" , "z.juv1.3" ,"z.juv2.1", "z.juv2.2", "z.juv2.3", "z.adult.1", "z.adult.2", "z.adult.3","theta.1", "theta.2",'Recruitment','Model') true.data$parameter <- row.names(true.data) true.data2 <- subset(true.data,parameter %in% new.names) true.data2 <- rename.vars(true.data2,from='parameter',to='variable',info=FALSE) ADMBnames <- c("z.1", "z.2", "z.3", "z.4", "z.5", "z.6", "z.7", "z.8", "z.9", "theta.1", "theta.2", "recrate") L.ADMB <- lapply(sims, function(X) with(X,ADMB$coef[ADMBnames])) admb.sim <- data.frame(do.call("rbind", L.ADMB),model='ADMBmle') BUGSnames <- c("z1", "z2", "z3", "z4", "z5", "z6", "z7", "z8", "z9","theta1", "theta2","recrate") L.BUGS <- lapply(sims, function(X) with(X,BUGS$coef[BUGSnames])) bugs.sim <- data.frame(do.call("rbind", L.BUGS),model='JAGS') admb.sim <- rename.vars(admb.sim,from=names(admb.sim),to=new.names,info=FALSE) bugs.sim <- rename.vars(bugs.sim,from=names(bugs.sim),to=new.names,info=FALSE) sim.sum <- rbind(admb.sim,bugs.sim) sim.sum.wide <- melt(sim.sum) ## use wildcard z.plots <- with(sim.sum.wide, subset(sim.sum.wide, subset = grepl(glob2rx("z.*"), variable))) true.data.z <- with(true.data2, subset(true.data2, subset = grepl(glob2rx("z.*"), variable))) other.plots <- subset(sim.sum.wide,!variable %in% unique(z.plots$variable)) true.data.other <- subset(true.data2,!variable %in% unique(true.data.z$variable)) x <- ggplot(z.plots, aes(x=Model, y=value, fill=Model)) + geom_boxplot()+stat_summary(fun.y=mean, geom="point", shape=5, size=4)+ facet_wrap(~variable, ncol=3) + theme_bw() + ylab("estimated value") + geom_hline(aes(yintercept = value,group=variable), data=true.data.z, colour="blue", linetype="dashed") print(x) @ \caption{Simulation results of select parameters from obtained using JAGS and ADMB maximum likelihood estimation. Subscripts to total mortality parameter estimates estimates ($Z_{s,d}$), are subscripted by stage and decade. Diamonds are mean of the estimates and dashed blue line is the true value of each parameter.}. \label{simsplot1} \end{figure} \begin{figure}[!htbp] \centering <>= x1 <- ggplot(other.plots, aes(x=Model, y=value, fill=Model)) + geom_boxplot()+stat_summary(fun.y=mean, geom="point", shape=5, size=4)+ facet_wrap(~variable, ncol=3, scales="free") + theme_bw() + ylab("estimated value") + geom_hline(aes(yintercept = value,group=variable), data=true.data.other, colour="blue", linetype="dashed")+theme(axis.text.x=element_text(size=7)) print(x1) @ \caption{Simulation results of select parameters from obtained using JAGS and ADMB maximum likelihood estimation. Diamonds are mean of the estimates and dashed blue line is the true value of each parameter.} \label{simsplot2} \end{figure} \newpage \section{Data} \tiny <>= cap <- 'Winter skate (\\textit{Leucoraja ocellata}) abundance data from Swain et al. (2009), \\textit{Ecol. App}. Data are millions of individuals (both sexes)' label1 <- 'table:skatedata' data1<-read.table("../DATA/skate_bugs.dat",skip=1) colnames(data1)<-c("Year", "Juvenile stage 1", "Juvenile stage 2", "Adult stage") table.x <- xtable(data1,caption=cap,align='llccc',label=label1) digs <- matrix(c(rep(0, nrow(data1)),rep(3,nrow(data1)*3)),nrow=nrow(data1),ncol=ncol(data1)) digs <- cbind(0,(digs)) # adds column to cover parameter names digits(table.x) <- digs print(table.x,type='latex',include.rownames=FALSE,sanitize.text.function = function(x){x},caption.placement="top",table.placement='!htbp') # USe long table because @ \normalsize <>= cap <- 'Fixed parameter values used to generate simulated data for comparison between software platforms.' label1 <- 'table:skatesim' q <- with(true.data, subset(true.data, subset = grepl(glob2rx("q*"), parameter)))$value r <- c(with(true.data, subset(true.data, subset = grepl(glob2rx("Recr*"), parameter)))$value,NA,NA) theta <- c(NA,with(true.data, subset(true.data, subset = grepl(glob2rx("theta*"), parameter)))$value) sdpro <- with(true.data, subset(true.data, subset = grepl(glob2rx("sdpro*"), parameter)))$value sdobs <- with(true.data, subset(true.data, subset = grepl(glob2rx("sdobs*"), parameter)))$value z1970 <- subset(true.data, parameter %in% c('z.juv1.1','z.juv2.1','z.adult.1'))$value z1980 <- subset(true.data, parameter %in% c('z.juv1.2','z.juv2.2','z.adult.2'))$value z1990 <- subset(true.data, parameter %in% c('z.juv1.3','z.juv2.3','z.adult.3'))$value Nadult <- c(NA,NA,subset(true.data, parameter %in% c('N.adult.1969'))$value) N69 <- c(subset(true.data, parameter %in% c('N.juv1.1969','N.juv2.1969'))$value,NA) temp1 <- data.frame(q,r,theta,sdpro,sdobs,z1970,z1980,z1990,Nadult,N69) temp2 <- t(temp1) par.names<- c("catchability scaling parameter (\\textit{q})","recruitment rate (\\textit{r})","stage transition prob. ($\\theta$)","process error std. dev. ($\\sigma$)","observation error std. dev. ($\\tau$)","inst. total mortality ($Z_{\\textrm{1970-1979}}$)","inst. total mortality ($Z_{\\textrm{1980-1989}}$)","inst. total mortality ($Z_{\\textrm{1990-2004}}$)","adult abundance$_{1965-1969}$ (millions)","juvenile abundance$_{1969}$ (millions)") table.x <- xtable(data.frame(Parameter=par.names,temp2),caption=cap,align='llccc',label=label1) colnames(table.x) <- c("Parameter","Juvenile stage 1", "Juvenile stage 2", "Adult stage") digs <- matrix(1,nrow=10,ncol=5) digs[2,4] <- 0 digs[2,5] <- 0 digs[3,3] <- 0 digs[9,3] <- 0 digs[9,4] <- 0 digs[10,5] <- 0 digits(table.x) <- digs print(table.x,type='latex',include.rownames=FALSE,sanitize.text.function = function(x){x},caption.placement="top",table.placement='!htbp') @ \newpage \section{ADMB Code} The AD Model Builder template code (\code{skate.tpl})\\ %\VerbatimInput[frame=single,numbers=left,samepage=false,fontsize=\tiny]{../ADMB/skate.tpl} \lstinputlisting[basicstyle=\tiny,language=C++,breaklines=TRUE,numbers=left,frame=single]{../ADMB/skate.tpl} \newpage \section{BUGS Code} The JAGS (BUGS) template code (\code{skate\_bugs.txt})\\ \lstinputlisting[basicstyle=\tiny,language=S,breaklines=TRUE,numbers=left,frame=single]{../BUGS/skate_bugs.txt} \newpage \section{Simulation Code} R simulation code skate\_sims\_run.R) \lstinputlisting[basicstyle=\tiny,language=R,breaklines=TRUE,numbers=left,frame=single]{../SIMS/sim.R} \newpage \bibliography{../../nonlin} \end{document}