\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{graphicx}
\DeclareGraphicsExtensions{.pdf}
\newcommand{\code}[1]{{\tt #1}}
\bibliographystyle{chicago}

\title{Tadpole writeup}
\author{Ben Bolker}
\date{\today}
\begin{document}
\maketitle

\SweaveOpts{keep.source=TRUE}
<<echo=FALSE>>=
options(continue=" ")
if (require(cacheSweave,quietly=TRUE)) {
  setCacheDir("tadpole_cache")
}
require(ggplot2)
require(grid) ## for unit()
zmargin <- opts(panel.margin=unit(0,"lines"))
theme_update(theme_bw())
library(rbenchmark)
@ 

\tableofcontents

\section{Summary}

This is a relatively simple nonlinear fitting problem.  
With reasonable starting values (which are fairly easy to guess by
direct graphical examination of the data), none of the approaches
is problematic. 
This example is a good introduction to the basics of input (formatting
and data), running, and output (retrieving results) using
the different approaches.
\begin{itemize}
\item ADMB works fastest and has the most
reliable 'basic' (quadratic) confidence intervals.
\item With box constraints set, the standard optimizers in R also work fine,
and quickly.  
\item BUGS (here using JAGS)
is much slower but still, for this simple problem, pretty quick and robust,
although there is no clear advantage to using BUGS for this problem
(unless e.g. one wanted to introduce informative priors).
\end{itemize}

\section{Introduction}

The data are originally from \cite{VoneshBolker2005} (these data are also
described and analyzed in \cite{Bolker2008}), describing the
numbers of reed frog (\emph{Hyperolius spinigularis})
tadpoles killed by predators as a function of size in
a small-scale field trial. (\code{TBL} is total body length in mm,
\code{Kill} is the number killed out of 10 tadpoles exposed
to predation). Figure~\ref{fig:rfsp1} shows the data.

Our main interest is in a quantitative description of the
``window of vulnerability'' --- the unimodal pattern of 
proportion killed as a function of size.  In various contexts,
we can use this description either to describe and test differences
among treatments (e.g., does the window of vulnerability differ
by predator size, or with tadpoles exposed to different predator cues?)
or to project the effects of growth and mortality rates through a
life stage (see references above and \cite{mccoy_predicting_2011}
for more details and examples).


This is one of the most basic examples developed by the NCEAS nonlinear modeling
working group, with few of the more challenging components incorporated in
the more complex examples.  In particular, the data
are
\begin{itemize}
\item small (there were 3 trials for each of 6 size classes);
\item well-behaved (no missing data, no extreme outliers);
\item fully specified (no observation errors, latent variables, or random effects);
\item low-dimensional (a single predictor variable (TBL) and a single response (number killed);
\item easy to describe mechanistically (it is reasonable to assume a binomial distribution
  or some variant of it).
\end{itemize}

\section{Basics}

Load all the packages we'll need, up front
(not all of these are absolutely necessary, but
it will be most convenient to make sure you have
them installed now).

<<getlibs>>=
library(ggplot2)    ## pictures
theme_update(theme_bw()) ## I dislike the default gray background
library(bbmle)      ## MLE fitting in R (wrapper)
library(optimx)     ## additional R optimizers
library(MCMCpack)   ## for post-hoc MCMC in R
library(coda)       ## analysis of MCMC runs (R, BUGS, ADMB)
library(scapeMCMC)  ## prettier MCMC diagnostics
library(R2admb)     ## R interface to AD Model Builder
library(R2jags)     ## R interface to JAGS (BUGS dialect)
source("../R/tadpole_R_funs.R")
@ 

This version uses \Sexpr{R.version$version.string} and package versions:
<<echo=FALSE>>=
usedpkgs <- sort(c("ggplot2","bbmle","optimx","MCMCpack","coda","R2admb","R2jags"))
i1 <- installed.packages()
print(i1[usedpkgs,"Version"],quote=FALSE)
@ 

The data are very simple:

<<dat1>>=
(ReedfrogSizepred <- read.table("../DATA/tadpole.dat"))
@ 

\begin{figure}
  \begin{center}
<<fig1,echo=FALSE,fig=TRUE>>=
g1  <- ggplot(ReedfrogSizepred,
             aes(x=TBL,y=Kill/10))+
  geom_point()+stat_sum(aes(size=factor(..n..)))+
  geom_smooth()+
  theme_bw()+
  labs(size="n",x="Size (total body length)",
       y="Proportion killed")+
  coord_cartesian(ylim=c(-0.05,0.55))
startest <- stat_function(fun = function(x) { 0.45*((x/13)*exp(1-x/13)) },
                          lty=2,colour="red")
print(g1+startest + xlim(0,38)+
      geom_vline(xintercept=13,colour="darkgray",lty=2)+
      geom_hline(yintercept=0.45,colour="darkgray",lty=2))
##      geom_line(data=predfr,colour="purple",lty=2))
@ 
\end{center}
\caption{Proportions of reed frogs killed by predators,
  as a function of total body length in mm.
  Red: starting estimate. Blue/gray fill: nonparametric
  (loess) fit.}
\label{fig:rfsp1}
\end{figure}


We need a nonlinear function that is flexible enough
to increase, decrease, and change its shape.  We chose
\begin{equation}
P(\mbox{kill}) = c ((S/d) \exp(1-(S/d)))^g ,
\end{equation}
sometimes called the ``power-Ricker'' function (in its simplest
form, it is a variant of the Ricker function $y=x\exp(-x)$ raised
to a power: $y=\left(x\exp(-x)\right)^g$.  The $d$ and $c$ parameters
adjust the scale of the function in the $x$ and $y$ directions respectively.
A bit of calculation or numerical experimentation will show that this
function is equal to zero when $S=0$; declines to zero as $S$ gets large;
initially increases as $c (S/d)^g$ (e.g. linearly if $g=1$ (Ricker),
quadratically if $g=2$); and peaks when $S=d$ at a height $c$.
\footnote{A different form of ``generalized Ricker'' was used
in \cite{VoneshBolker2005}; as described by \cite{Bolker2008}), that turned
out to be a poor choice, because the same data could be described in
quite different ways by a single function (i.e., the likelihood
surface has multiple maxima.}.
Looking at the plot (a luxury we have in the low-dimensional case),
we can see that a reasonable starting set of estimates would be
\begin{itemize}
  \item $g=1$ (the Ricker model; linear increase near $S=0$);
  \item $c=0.45$, 
  \item $d=13$
\end{itemize}
(Figure~\ref{fig:rfsp1} illustrates the initial values).

The remainder of the model formulation describes the
stochastic part of the model. Although it is slightly questionable
whether tadpoles within a tank are identically vulnerable and
killed independently of each other, we will go ahead and
use a binomial model anyway:
\begin{equation}
  \mbox{Killed} \sim \mbox{Binomial}(P(\mbox{kill}),N)
\end{equation}

<<libs,echo=FALSE>>=
library(ggplot2) ## for pictures
@ 

\section{R}

\subsection{Fitting}
In R we'll use the \code{mle2} function from the \code{bbmle} package.  There is no
particular computational advantage to \code{mle2} over the various optimizers
in R (\code{optim}, \code{nlm}, and \code{nlminb} in base R, or the variants
found in the \code{optimx} package) --- it is just a wrapper that provides
various conveniences (profiling, confidence intervals, predictions, etc.) for
maximum likelihood estimation problems.

Here is the code to fit a binomial model with
\code{mle2} using these starting points:
<<mlefit>>=
library(bbmle)
tadpole_R_fit
mle2_fit <- tadpole_R_fit()
@ 

<<echo=FALSE>>=
## omit parscale in both cases -- not supported
mle2_fit_nlminb <- mle2(Kill~dbinom(c*((TBL/d)*exp(1-TBL/d))^g,size=10),
                        start=list(c=0.45,d=13,g=1),data=ReedfrogSizepred,
                        optimizer="nlminb",
                        lower=c(c=0.003,d=10,g=0),
                        upper=c(c=0.8,d=20,g=20))

library(optimx)
mle2_fit_bobyqa <- mle2(Kill~dbinom(c*((TBL/d)*exp(1-TBL/d))^g,size=10),
                            start=list(c=0.45,d=13,g=1),data=ReedfrogSizepred,
                            optimizer="optimx",method="bobyqa",
                            lower=c(c=0.003,d=10,g=0),
                            upper=c(c=0.8,d=20,g=20))

@ 
\emph{Notes}:
\begin{itemize}
\item It is not absolutely necessary to set bounds on the optimization in this case, but it
  is generally good practice and will help avoid many problems when the 
  optimizer wants to wander off to strange values on its way to its ultimate
  best solution (which may be perfectly sensible).
\item The tradeoff for setting bounds is that there is a more 
  limited set of optimizers available. I used \code{L-BFGS-B}, which is
  available within \code{optim}, but is finicky
  \begin{itemize}
  \item it may fail if your parameters are not scaled appropriately
    (which you can do by hand by redefining
    your parameters so they all have an expected magnitude between 1 and 10, or as above by setting
    \code{parscale} to the approximate magnitude of each parameter);
  \item it will fail if it encounters a non-finite value (\code{Inf}, \code{NA}, or \code{NaN})
    when calling the negative log-likelihood function; 
  \item it sometimes steps slightly
    over the specified lower/upper bounds when computing the derivatives of the function,
    so it is best to set these slightly in from any set of parameter
    values that will return a non-finite value (e.g. set the boundary for $c$ at 0.002 or
    above rather than exactly at zero so the probability always stays positive).
  \end{itemize}
  Alternatives to \code{L-BFGS-B} within R include \code{nlminb} (base R) 
  and \code{bobyqa} (derivative-free, in \code{optimx}) [neither supports
  \code{parscale}].
  In a quick test, \code{nlminb} and \code{L-BFGS-B} were approximately the same
  speed, while \code{bobyqa} (which may be more robust in some situations) was
  about 12 times slower \ldots
\end{itemize}
%\fixme{add examples (in footnote?) They are currently present in an Sweave chunk,
%but not printed \ldots}

Printing out the fit gives just the original function call, the coefficients, and the
log-likelihood:
<<>>=
mle2_fit
@ 

\code{summary(...)} gives estimates, approximate standard errors (based on a quadratic
approximation to the likelihood surface), $Z$ values (estimate/std. error), and
\emph{Wald $p$-values} against the null hypotheses that each parameter is 
separately equal to zero (based again on a quadratic approximation).
\code{coef(summary(...))} will extract the table printed here.

<<>>=
summary(mle2_fit)
@ 

Other accessor methods allow you to do inference (\code{AIC}, \code{logLik}, \code{deviance}, 
\code{anova} (Likelihood Ratio Test vs alternative models) \ldots)  %\fixme{how much should I say here?}

\code{mle2} allows you to generate predicted values, \emph{if} you have used
the formula interface rather than writing your own objective function:
<<predvals>>=
TBLvec = seq(9.5,36,length=100)
predfr <- data.frame(TBL=TBLvec,
                     Kill=predict(mle2_fit,
                       newdata=data.frame(TBL=TBLvec,
                         Exposed=10)))
@ 

\subsection{Profiling}
Compute likelihood profiles and profile confidence intervals:
<<profile>>=
mle2_profile <- profile(mle2_fit)
mle2_profile_confint <- confint(mle2_profile)
## or just confint(mle2_fit) if you don't want to do anything 
## else with the profile (such as plotting it or extracting 
## confidence intervals at multiple alpha levels)
@

\subsection{MCMC}

Using the \code{MCMCpack} package, which contains a general-purpose
Metropolis-Hastings sampler (\code{MCMCmetrop1R}), we can get a \emph{post hoc} MCMC chain by
extracting the negative log-likelihood function from the fitted model
and putting a wrapper around it to (1) pass the parameters as a list
rather than a vector (\code{do.call(...,as.list(p))}); (2) compute the (positive) log-likelihood
rather than the negative log-likelihood; (3) intercept bad (non-finite)
values and replace them with large negative values (the specific ``bad''
value given here might have to be adjusted for other applications).

We specify the starting values from our previous fit; the length
of the chain; and the thinning fraction.


<<mcmc>>=
tadpole_R_mcmc
mle2_mcmc <- tadpole_R_mcmc(mle2_fit)
colnames(mle2_mcmc) <- names(coef(mle2_fit))
@ 

The default \code{coda} function for trace plots shows the plots in a one-column, stacked
layout: this is generally good for longitudinal data like trace plots, but can be hard to
read for models (unlike this one) with lots of parameters. Use the \code{layout=c(m,n)} argument
for an $m \times n$ layout; \code{aspect="fill"} to allow the aspect ratio of the subplots
to be flexible); and \code{as.table=TRUE} if you want the parameters arranged top-to-bottom
rather than bottom-to-top.

<<traceplot1,fig=TRUE>>=
print(xyplot(mle2_mcmc,as.table=TRUE))
@ 

\code{scapeMCMC} offers a slightly prettier version of the trace plot
(although its versions do not work with multiple MCM chains):
<<traceplot2,fig=TRUE>>=
plotTrace(mle2_mcmc,layout=c(1,3))
@ 

<<>>=
gd <- geweke.diag(mle2_mcmc)                ## Z-scores
(gp <- 2*pnorm(abs(gd$z),lower.tail=FALSE)) ## p values
effectiveSize(mle2_mcmc)
@ 

We are aiming for an effective size of 1000 and non-significant
$p$ values from the Geweke diagnostic.

\section{AD Model Builder (via R2admb)}

Here is a \code{minimal} TPL (AD Model Builder definition) file:

\VerbatimInput[frame=single,numbers=left,samepage=true,fontsize=\small]{../ADMB/tadpole.tpl}

\begin{itemize}
\item Comments are written in C++ format: everything on a line after \code{//} is ignored. 
  \item lines 1--4 are the \code{PARAMETER} section; most of the parameters
    will get filled in automatically by \code{R2admb} based on the
    input parameters you specify, but you should include
    this section if you need to define any additional
    utility variables.
    In this case we define \code{prob} as a vector indexed from
    1 to \code{nobs} (we will specify \code{nobs},
    the number of observations, in our data list).
  \item most of the complexity of the \code{PROCEDURE} section
    (lines 7 and 11--14) has to do with making sure that the
    mortality probabilities do not exceed the range (0,1), which is
    not otherwise guaranteed by this model specification.  Line 7 defines
    a utility variable \code{fpen}; lines 11--14 use the built-in
    ADMB function \code{posfun} to adjust low probabilities
    up to 0.001 (line 11) and high probabilities down to 0.999
    (line 13), and add appropriate penalties to the negative
    log-likelihood to push the optimization away from these
    boundaries (lines 12 and 14).
    %\fixme{is this really necessary? leave it in anyway?}
  \item the rest of the \code{PROCEDURE} section
    simply computes the mortality
    probabilities as $c ((S/d) \exp(1-(S/d)))^g$
    as specified above (line 9) and computes the binomial
    log-likelihood on the basis of these
    probabilities (lines 16-18). Because this is a log-likelihood
    and we want to compute a negative log-likelihood, we \emph{subtract} it 
    from any penalty terms that have already accrued.
    The code is written in C++ syntax,
    using \code{=} rather than \verb+<-+ for assignment, \code{+=} to increment
    a variable and \code{-=} to decrement one. The power operator is 
    \code{pow(x,y)} rather than \verb+x^y+; elementwise multiplication of
    two vectors uses \code{elem\_prod} rather than \code{*}.
  \end{itemize}
    
<<admbfit_getruns,echo=FALSE>>=
## load in data from previously executed runs
zz <- load("../ADMB/tadpole_ADMB_fits.RData")
@ 

To run this model, we save it in a text
file called \code{tadpole.tpl};
run \verb+setup_admb()+ to locate the
AD Model Builder binaries and libraries on
our system; and run \verb+do_admb+ with appropriate
arguments.

<<admbfit_fun,echo=FALSE>>=
source("../ADMB/tadpole_ADMB_funs.R")
tadpole_ADMB_fit
@ 

The \code{data}, \code{params}, and \code{bounds} (parameter bounds)
arguments should be reasonably self-explanatory.
When \code{checkparam="write"} and \code{checkdata="write"} are
specified, \code{R2admb} attempts to write appropriate DATA
and PARAMETER sections into a modified TPL file, leaving the
results with the suffix \code{\_gen.tpl} at the end of the run.
(In this example, we could leave out all of the \code{DATA}
section and everything in
the \code{PARAMETER} sections except the
utility variable \code{prob}.)

%You might instead choose to compose the whole TPL file yourself,
%in which case you can add comments appropriately:
%\VerbatimInput[frame=single,numbers=left,fontsize=\small]{ReedfrogSizepred.tpl}

Now that we have fitted the model, here are some of
the things we can do with it:

\begin{itemize}
  \item Get basic information about the fit and coefficient estimates:
<<basic>>=
tfit_admb
@ 
\item Get vector of coefficients only:
<<coef>>=
coef(tfit_admb)
@
\item Get a coefficient table including standard errors and (approximate!!)
  $p$ values:
<<summary>>=
summary(tfit_admb)
@ 
(you can use \code{coef(summary(tfit\_admb))} to extract just the table).
@ 
\item Variance-covariance matrix of the parameters:
<<vcov>>=
vcov(tfit_admb)
@ 
Log-likelihood, deviance, AIC:
<<others>>=
c(logLik=logLik(tfit_admb),deviance=deviance(tfit_admb),
  AIC=AIC(tfit_admb))
@ 
\end{itemize}

\subsection{Profiling}

You can also ask ADMB to compute likelihood profiles for a model.
If you code it yourself in the TPL file you need to add variables
of type \code{likeprof\_number} to keep track of the values:
\code{R2admb} handles these details for you.  You just need to
specify \code{profile=TRUE} and give a list of the parameters
you want profiled.

<<profrun,eval=FALSE>>=
tfit_admb_prof <- do_admb("ReedfrogSizepred0",
               data=c(list(nobs=nrow(ReedfrogSizepred),
                 nexposed=rep(10,nrow(ReedfrogSizepred))),
                 ReedfrogSizepred),
               params=list(c=0.45,d=13,g=1),
               bounds=list(c=c(0,1),d=c(0,50),g=c(-1,25)),
               run.opts=run.control(checkparam="write",
                 checkdata="write"),
               profile=TRUE,
               profpars=c("c","d","g"),
               admb_errors="warn")
@ 

The profile information is stored in a list \verb+tfit_admb_prof$prof+ with entries
for each variable to be profiled.  Each entry in turn contains a list
with elements \code{prof} (a 2-column matrix containing the parameter value
and profile log-likelihood), \code{ci} (confidence intervals derived from
the profile), \code{prof\_norm} (a profile based on the normal approximation),
and \code{ci\_norm} (confidence intervals, ditto).

Let's compare ADMB's profiles to those generated from R:

<<mleprof,cache=TRUE>>=
m0prof <- profile(mle2_fit)
@ 

(A little bit of magic [hidden] gets everything into the same data frame
and expressed in the same scale that R uses for profiles, which is the square root
of the change in deviance ($-2L$) between the best fit and the profile:
this scale provides a quick graphical assessment of the profile
shape, because quadratic profiles will be {\sf V}-shaped on this scale.)

<<profcalcs2,echo=FALSE>>=
tmpf <- function(p,w="prof") {
  pp <- log(tfit_admb_prof$prof[[p]][[w]][,2])
  pp <- max(pp)-pp
  data.frame(param=p,z=sqrt(2*pp),
             par.vals.c=NA,par.vals.d=NA,
             par.vals.g=NA,focal=tfit_admb_prof$prof[[p]][[w]][,1])
}
proflist <- do.call(rbind,lapply(list("c","d","g"),tmpf))
profnlist <- do.call(rbind,lapply(list("c","d","g"),tmpf,w="prof_norm"))
pdat <- rbind(cbind(as.data.frame(m0prof),method="mle2"),
              cbind(proflist,method="ADMB"),
              cbind(profnlist,method="ADMB_norm"))
@ 

\SweaveOpts{width=10,height=6}
<<profpic,echo=FALSE,fig=TRUE>>=
print(ggplot(pdat,aes(x=focal,y=abs(z),group=method,colour=method))+geom_line()+
      geom_point(alpha=0.5)+
      facet_grid(.~param,scale="free_x")+
      ylim(0,3)+
      labs(x="",y=expression(Delta(sqrt(-2*L))))+
      geom_hline(yintercept=1.96,lty=2)+zmargin)
@ 
\SweaveOpts{width=7,height=7}

Notice that R evaluates the profile at a smaller number of locations,
using spline interpolation to compute confidence intervals.

\subsection{MCMC}

Another one of ADMB's features is that it can use Markov chain
Monte Carlo (starting at the maximum likelihood estimate and
using a candidate distribution based on the approximate sampling
distribution of the parameters) to get more information about the
uncertainty in the estimates.  This procedure is especially helpful
for complex models (high-dimensional or containing random effects)
where likelihood profiling becomes problematic.

To use MCMC, just add \code{mcmc=TRUE} and specify the parameters
you want to keep track of with \code{mcmcpars}:

<<admbfakemc,eval=FALSE>>=
tfit_admb_mcmc <- do_admb("ReedfrogSizepred",
              data=c(list(nobs=nrow(ReedfrogSizepred),
                nexposed=rep(10,nrow(ReedfrogSizepred))),
                ReedfrogSizepred),
                params=list(c=0.45,d=13,g=1),
                bounds=list(c=c(0,1),d=c(0,50),g=c(-1,25)),
                run.opts=run.control(checkparam="write",
                  checkdata="write"),
                mcmc=TRUE,
                mcmc.opts=mcmc.control(mcmcpars=c("c","d","g")))
## clean up leftovers:
unlink(c("reedfrogsizepred0.tpl",
         "reedfrogsizepred0_gen.tpl",
         "reedfrogsizepred0"))
@ 

The output of MCMC is stored in two ways. 

(1) ADMB internally computes a histogram of the MCMC sampled
density; this is stored in a list element called \verb+$hist+,
as an object of class \code{admb\_hist}.  It has its own plot method:

<<printhist,fig=TRUE>>=
print(plot(tfit_admb_mcmc$hist))
@ 

(2) Alternatively, the full set of samples is stored (as
a data frame) in list element \verb+$mcmc+.  If you load
the \code{coda} package, you can convert this into an object
of class \code{mcmc}, and then use the various methods
implemented in \code{coda} to analyze it.


<<codaload>>=
library(coda)
mmc <- as.mcmc(tfit_admb_mcmc$mcmc)
@ 

Highest posterior density (i.e. Bayesian credible) intervals:
<<HPDinterval>>=
HPDinterval(mmc)
@ 

Density plots:
<<admbdensplots,fig=TRUE>>=
print(densityplot(mmc,layout=c(3,1)))
@ 

Trace plots:
<<admbtraceplots,fig=TRUE>>=
print(xyplot(mmc,as.table=TRUE))
@ 

(You don't need to use \code{print} to see these plots in 
an interactive session --- it's just required for generating documents.)

See the \code{coda} and \code{scapeMCMC} packages for more information on quantitative and graphical diagnostics for MCMC.

The jaggedness of this trace plot indicates that the chain is not actually
mixing very well. ADMB does offer some options for tuning the chain: see
the ADMB documentation for more details. 
If you have an ADMB-compiled executable, then
(e.g. \code{./tadpole -? | grep mcmc} at a Unix or MacOS terminal command line
will give you a terse list of the possibilities:
\begin{verbatim}
 -mcdiag         use diagonal covariance matrix for mcmc with diagonal values 1
 -mcmc [N]       perform markov chain monte carlo with N simulations
 -mcmult N       multiplier N for mcmc default
 -mcr            resume previous mcmc
 -mcrb  N        reduce the amount of correlation in the covariance matrix 1<=N<=9
 -mcnoscale      don't rescale step size for mcmc depending on acceptance rate
 -nosdmcmc       turn off mcmc histogram calcs to make mcsave run faster
 -mcgrope N      use probing strategy for mcmc with factor N
 -mcseed N       seed for random number generator for markov chain monte carlo
 -mcscale N       rescale step size for first N evaluations
 -mcsave N       save the parameters for every N'th simulation
 -mceval         Go throught the saved mcmc values from a previous mcsave
\end{verbatim}

\section{BUGS}

The BUGS input file is very simple:

\VerbatimInput[frame=single,numbers=left,samepage=true,fontsize=\small]{../BUGS/tadpole_bugs.txt}

\begin{itemize}
\item BUGS is not vectorized --- all calculations on vectors must be written out
  with explicit \code{for} loops
\item BUGS uses different names and (more dangerously) different
  parameterizations for distributions even when they have the same names (see \url{http://tinyurl.com/bugsparams}):
  the most common trap is that it parameterizes the normal distribution by the \emph{precision} (inverse variance)
  rather than the standard deviation
\item the priors used here are a bit unusual -- we set them as uniform priors with the
same range as the box constraints on the parameters used in the previous fits.
Bayesians would more typically use a Normal distribution with a large variance
(small precision, as specified in BUGS) for parameters that need not be positive,
and a log-Normal, Gamma or other positive distribution with a large variance
for vague priors on parameters that must be positive.  There is much discussion
of priors: if you're feeling lost, start with \cite{McCarthy2007} for a relatively
gentle discussion, and then continue with other references
\citep{Gelman+1996,Lambert2005}.
Gelman and co-authors
\citep{GelmanHill2006,gelman_prior_2006}
warn about situations where the traditional Gamma prior for Normal precisions
can be problematic, such as in estimating variances from small samples.
The debate over priors continues in the blogosphere (see e.g. \url{http://tinyurl.com/gelmanprior}
and other posts about priors on Andrew Gelman's blog).
\item it is general better practice to run multiple chains (i.e. to run BUGS from several
different starting points); it makes slightly stronger tests of convergence, such as
the Gelman-Rubin statistic, possible (see e.g. the Weeds example).  (This setting
goes in the R code to call BUGS, not in the BUGS model file itself.)
\end{itemize}

Using JAGS:

<<jags1>>=
jags_fit <- jags(model="../BUGS/tadpole_bugs.txt",
                 data=c(list(N=10),as.list(subset(ReedfrogSizepred,
                          select=-c(Exposed)))),
                 parameters.to.save=c("c","d","g"),
                 n.chains=2,
                 inits=list(list(c=0.45,d=13,g=1),
                   list(c=0.4,d=10,g=2)),
                 progress.bar="none")
jags_fit
jags_fit_mcmc <- as.mcmc(jags_fit)
@ 

<<jagssum>>=
summary(jags_fit_mcmc)
@ 

Since \verb+jags_fit_mcmc+ is an \code{mcmc} object, we can run the same plots
and diagnostics as on the previous MCMC runs: \code{xyplot}, \code{densityplot},
\code{geweke.diag}, \code{effectiveSize}, \code{HPDinterval}, etc..

<<jags_traceplot,fig=TRUE>>=
source("../../TOOLS/misc_funs.R")  ## get dropdev() tool to ignore deviance
print(xyplot(dropdev(jags_fit_mcmc),as.table=TRUE))
@ 

This is better behaved than the ADMB trace plot, although we would
like it to look still more like white noise \ldots

\section{Simulation results}

(See figures.)
<<loadsims,echo=FALSE>>=
setwd("../../TOOLS")  ## move up to projects
source("tools.R")
setwd("..")
p <- get_proj("tadpole")
setwd("tadpole/WRITEUP")  ## restore original directory
S <- with(p,simsum(sim,true,clip=5))
B <- basesum(p)
@ 

\begin{figure}
  \centering
<<timefig,echo=FALSE,fig=TRUE>>=
print(S$plots$time)
@   
  \caption{Timings}
  \label{fig:sim_timings}
\end{figure}

\begin{figure}
  \centering
<<paramfig,echo=FALSE,fig=TRUE>>=
print(S$plots$params)
@   
  \caption{Distribution of parameter estimates}
  \label{fig:sim_params}
\end{figure}


\begin{figure}
  \centering
<<fullparamfig,echo=FALSE,fig=TRUE>>=
print(S$plots$fullparams)
@   
  \caption{Full distribution of parameters+CIs}
  \label{fig:sim_fullparams}
\end{figure}


\begin{figure}
  \centering
<<coverfig,echo=FALSE,fig=TRUE>>=
print(S$plots$cisum)
@   
  \caption{Confidence interval summaries}
  \label{fig:cisum}
\end{figure}

\bibliography{../../nonlin}

\end{document}