\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}
\usepackage{Sweave}
\begin{document}
\maketitle


%Below are some of the basic R-packages we use to run the models and retrieve results into the R environment.
\begin{Schunk}
\begin{Sinput}
> 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
\end{Sinput}
\end{Schunk}

% Now load in the workspaces for the relevant objects

\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}
% latex table generated in R 2.15.2 by xtable 1.7-0 package
% Wed Mar 20 16:12:48 2013
\begin{table}[!htp]
\begin{center}
\caption{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.}
\label{table:runtimes}
\begin{tabular}{rrr}
  \hline
ADMB & ADMBmcmc & JAGS \\ 
  \hline
13.0 & 61.6 & 174.2 \\ 
   \hline
\end{tabular}
\end{center}
\end{table}\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:


\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 13 seconds to run.

\begin{center}
\begin{tiny}
% latex table generated in R 2.15.2 by xtable 1.7-0 package
% Wed Mar 20 16:12:48 2013
\begin{table}[!htp]
\begin{center}
\begin{tabular}{rrrrr}
  \hline
 & Estimate & Std. Error & z value & Pr($>$$|$z$|$) \\ 
  \hline
recrate & 3.11 & 0.01 & 350.29 & 0.00 \\ 
  theta.1 & 0.19 & 0.06 & 3.03 & 0.00 \\ 
  theta.2 & 0.25 & 0.06 & 3.91 & 0.00 \\ 
  z.juv1.1 & 1.24 & 0.12 & 10.77 & 0.00 \\ 
  z.juv1.2 & 0.56 & 0.13 & 4.50 & 0.00 \\ 
  z.juv1.3 & 0.52 & 0.10 & 5.12 & 0.00 \\ 
  z.juv2.1 & 0.28 & 0.16 & 1.76 & 0.08 \\ 
  z.juv2.2 & 0.29 & 0.15 & 2.03 & 0.04 \\ 
  z.juv2.3 & 0.47 & 0.11 & 4.44 & 0.00 \\ 
  z.adult.1 & 0.12 & 0.12 & 0.93 & 0.35 \\ 
  z.adult.2 & 0.22 & 0.11 & 2.04 & 0.04 \\ 
  z.adult.3 & 0.45 & 0.07 & 6.08 & 0.00 \\ 
   \hline
\end{tabular}
\caption{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}
\label{table:admbMLE}
\end{center}
\end{table}\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 61.6 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
\includegraphics{skate-admb_geweke1}

\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
\includegraphics{skate-admb_geweke2}
\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
\includegraphics{skate-ADMB_trace}
\caption{Traceplots of select ADMB-RE mcmc model parameters}
\label{admbmcmc:trace}

\end{figure}


% Below is the latex table results
\begin{center}
\begin{tiny}
% latex table generated in R 2.15.2 by xtable 1.7-0 package
% Wed Mar 20 16:12:51 2013
\begin{table}[!htbp]
\begin{center}
\begin{tabular}{rrrrrr}
  \hline
 & Estimate & Std.Error & 50\% & 2.5\% & 97.5\% \\ 
  \hline
recrate & 3.11 & 0.77 & 3.36 & 2.15 & 5.42 \\ 
  theta.1 & 0.19 & 0.04 & 0.21 & 0.14 & 0.29 \\ 
  theta.2 & 0.25 & 0.04 & 0.27 & 0.19 & 0.36 \\ 
  z.1.1 & 1.24 & 0.39 & 1.32 & 0.57 & 2.06 \\ 
  z.1.2 & 0.28 & 0.29 & 0.32 & 0.03 & 0.90 \\ 
  z.1.3 & 0.12 & 0.15 & 0.15 & 0.01 & 0.44 \\ 
  z.2.1 & 0.56 & 0.31 & 0.64 & 0.10 & 1.25 \\ 
  z.2.2 & 0.29 & 0.27 & 0.37 & 0.03 & 0.93 \\ 
  z.2.3 & 0.22 & 0.17 & 0.23 & 0.02 & 0.56 \\ 
  z.3.1 & 0.52 & 0.24 & 0.55 & 0.12 & 1.03 \\ 
  z.3.2 & 0.47 & 0.24 & 0.51 & 0.09 & 1.07 \\ 
  z.3.3 & 0.45 & 0.16 & 0.46 & 0.17 & 0.86 \\ 
   \hline
\end{tabular}
\caption{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}
\label{table:admbmcmc}
\end{center}
\end{table}\end{tiny}
\end{center}

\newpage
\section{BUGS (using JAGS)}
The JAGS model required 174.2 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
\begin{Schunk}
\begin{Sinput}
> 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))
\end{Sinput}
\end{Schunk}
\includegraphics{skate-jags_geweke1}
\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
\includegraphics{skate-jags_geweke2}
\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

\begin{Schunk}
\begin{Sinput}
>  xyplot(as.mcmc(sub.jags.df1))
\end{Sinput}
\end{Schunk}
\includegraphics{skate-jags_trace}
\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}
% latex table generated in R 2.15.2 by xtable 1.7-0 package
% Wed Mar 20 16:12:58 2013
\begin{table}[!htp]
\begin{center}
\begin{tabular}{rrrrrrrr}
  \hline
 & Mean & SD & Naive SE & Time-series SE & 50\% & 2.5\% & 97.5\% \\ 
  \hline
deviance & 319.08 & 14.77 & 0.30 & 0.29 & 320.48 & 285.85 & 343.96 \\ 
  recrate & 4.49 & 2.24 & 0.04 & 0.04 & 4.07 & 1.57 & 9.94 \\ 
  theta.1 & 0.21 & 0.04 & 0.00 & 0.00 & 0.21 & 0.14 & 0.30 \\ 
  theta.2 & 0.25 & 0.04 & 0.00 & 0.00 & 0.25 & 0.17 & 0.34 \\ 
  z.1.1 & 1.48 & 0.54 & 0.01 & 0.01 & 1.47 & 0.43 & 2.59 \\ 
  z.2.1 & 0.72 & 0.38 & 0.01 & 0.01 & 0.68 & 0.10 & 1.59 \\ 
  z.3.1 & 0.61 & 0.30 & 0.01 & 0.01 & 0.58 & 0.11 & 1.28 \\ 
  z.1.2 & 0.43 & 0.30 & 0.01 & 0.01 & 0.37 & 0.03 & 1.13 \\ 
  z.2.2 & 0.42 & 0.26 & 0.01 & 0.01 & 0.40 & 0.04 & 1.04 \\ 
  z.3.2 & 0.46 & 0.22 & 0.00 & 0.00 & 0.44 & 0.07 & 0.96 \\ 
  z.1.3 & 0.14 & 0.09 & 0.00 & 0.00 & 0.13 & 0.01 & 0.37 \\ 
  z.2.3 & 0.22 & 0.13 & 0.00 & 0.00 & 0.20 & 0.02 & 0.51 \\ 
  z.3.3 & 0.43 & 0.13 & 0.00 & 0.00 & 0.42 & 0.19 & 0.71 \\ 
   \hline
\end{tabular}
\caption{Select Parameter and associated standard deviation from JAGS. Subscripts to total mortality parameter estimates estimates ($Z_{s,d}$), are subscripted by stage and decade}
\label{table:jagsmcmc1}
\end{center}
\end{table}\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.



\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}).


\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}


We used a simulation approach to assess the ability of each modelling platform to accurately and
precisely estimate model parameters.  We simulated
25 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  
\includegraphics{skate-simsplot1}
\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
\includegraphics{skate-simsplot2}
\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
% latex table generated in R 2.15.2 by xtable 1.7-0 package
% Wed Mar 20 16:13:10 2013
\begin{table}[!htbp]
\begin{center}
\caption{Winter skate (\textit{Leucoraja ocellata}) abundance data from Swain et al. (2009), \textit{Ecol. App}.  Data are millions of individuals (both sexes)}
\label{table:skatedata}
\begin{tabular}{lccc}
  \hline
Year & Juvenile stage 1 & Juvenile stage 2 & Adult stage \\ 
  \hline
1970 & 0.226 & 0.034 & 3.880 \\ 
  1971 & 1.005 & 0.069 & 1.949 \\ 
  1972 & 0.300 & 0.139 & 3.750 \\ 
  1973 & 2.370 & 0.940 & 6.429 \\ 
  1974 & 1.247 & 0.239 & 1.622 \\ 
  1975 & 3.022 & 1.362 & 5.140 \\ 
  1976 & 0.837 & 0.494 & 0.764 \\ 
  1977 & 8.099 & 1.973 & 2.422 \\ 
  1978 & 1.341 & 0.881 & 1.397 \\ 
  1979 & 16.862 & 4.335 & 2.732 \\ 
  1980 & 2.880 & 0.851 & 2.491 \\ 
  1981 & 7.078 & 3.619 & 1.425 \\ 
  1982 & 2.887 & 2.364 & 3.407 \\ 
  1983 & 3.203 & 0.879 & 0.445 \\ 
  1984 & 4.409 & 0.587 & 0.778 \\ 
  1985 & 4.903 & 1.762 & 2.107 \\ 
  1986 & 3.268 & 0.456 & 1.171 \\ 
  1987 & 2.049 & 0.358 & 1.122 \\ 
  1988 & 6.347 & 2.992 & 3.994 \\ 
  1989 & 0.719 & 0.592 & 0.633 \\ 
  1990 & 1.825 & 1.808 & 2.357 \\ 
  1991 & 11.914 & 2.654 & 2.175 \\ 
  1992 & 5.842 & 2.163 & 2.663 \\ 
  1993 & 14.844 & 0.624 & 0.311 \\ 
  1994 & 10.202 & 0.323 & 0.357 \\ 
  1995 & 2.021 & 0.757 & 0.464 \\ 
  1996 & 2.170 & 2.102 & 1.203 \\ 
  1997 & 4.742 & 2.034 & 0.957 \\ 
  1998 & 3.559 & 0.174 & 0.342 \\ 
  1999 & 1.700 & 1.444 & 0.353 \\ 
  2000 & 0.351 & 2.021 & 1.127 \\ 
  2001 & 0.957 & 0.359 & 0.155 \\ 
  2002 & 2.360 & 0.166 & 0.155 \\ 
  2003 & 5.584 & 0.657 & 0.706 \\ 
  2004 & 0.059 & 0.277 & 0.620 \\ 
   \hline
\end{tabular}
\end{center}
\end{table}\normalsize

% latex table generated in R 2.15.2 by xtable 1.7-0 package
% Wed Mar 20 16:13:10 2013
\begin{table}[!htbp]
\begin{center}
\caption{Fixed parameter values used to generate simulated data for comparison between software platforms.}
\label{table:skatesim}
\begin{tabular}{lccc}
  \hline
Parameter & Juvenile stage 1 & Juvenile stage 2 & Adult stage \\ 
  \hline
catchability scaling parameter (\textit{q}) & 1.0 & 1.0 & 1.0 \\ 
  recruitment rate (\textit{r}) & 3.2 &  &  \\ 
  stage transition prob. ($\theta$) &  & 0.2 & 0.2 \\ 
  process error std. dev. ($\sigma$) & 0.3 & 0.4 & 0.2 \\ 
  observation error std. dev. ($\tau$) & 1.0 & 0.9 & 0.8 \\ 
  inst. total mortality ($Z_{\textrm{1970-1979}}$) & 1.4 & 0.4 & 0.1 \\ 
  inst. total mortality ($Z_{\textrm{1980-1989}}$) & 0.7 & 0.4 & 0.2 \\ 
  inst. total mortality ($Z_{\textrm{1990-2004}}$) & 0.6 & 0.4 & 0.4 \\ 
  adult abundance$_{1965-1969}$ (millions) &  &  & 4.0 \\ 
  juvenile abundance$_{1969}$ (millions) & 1.0 & 0.2 &  \\ 
   \hline
\end{tabular}
\end{center}
\end{table}\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}
