\documentclass[nojss]{jss}
\usepackage[utf8]{inputenc}
\usepackage{amsmath,booktabs}
%\VignetteIndexEntry{Using the runjags package}
\author{Matthew J. Denwood\\University of Copenhagen}
\Plainauthor{Matthew J. Denwood}
\title{\pkg{runjags}: An \proglang{R} Package Providing Interface Utilities, Model Templates, Parallel Computing Methods and Additional Distributions for MCMC Models in~\pkg{JAGS}}
\Plaintitle{runjags: An R Package Providing Interface Utilities, Model Templates, Parallel Computing Methods and Additional Distributions for MCMC Models in JAGS}
\Shorttitle{\pkg{runjags}: \pkg{JAGS} Interface Utilities and Additional Distributions}
\Abstract{
The \pkg{runjags} package provides a set of interface
functions to facilitate running Markov chain Monte Carlo models in
\pkg{JAGS} from within \proglang{R}. Automated calculation of
appropriate convergence and sample length diagnostics, user-friendly
access to commonly used graphical outputs and summary statistics,
and parallelized methods of running \pkg{JAGS} are
provided. Template model specifications can be generated using a
standard \pkg{lme4}-style formula interface to assist users less
familiar with the \proglang{BUGS} syntax. Automated simulation
study functions are implemented to facilitate model performance
assessment, as well as drop-$k$ type cross-validation studies, using
high performance computing clusters such as those provided by
\pkg{parallel}. A module extension for \pkg{JAGS} is also
included within \pkg{runjags}, providing the Pareto family of
distributions and a series of minimally-informative priors including
the DuMouchel and half-Cauchy priors. This vignette is taken from the
publication for the \pkg{runjags} package \citep{Denwood:2016ee}.
It outlines the primary functions of this package, and gives an
illustration of a simulation study to assess the performance of
equivalent model specifications.
}
\Keywords{MCMC, Bayesian, graphical models, interface utilities, \pkg{JAGS}, \proglang{BUGS}, \proglang{R}}
\Plainkeywords{MCMC, Bayesian, graphical models, interface utilities, JAGS, BUGS, R}
\Address{
Matthew J. Denwood\\
Department of Large Animal Sciences\\
Section for Animal Welfare and Disease Control\\
Faculty of Health and Medical Sciences\\
University of Copenhagen\\
Denmark\\
E-mail: \email{md@sund.ku.dk}\\
URL: \url{http://iph.ku.dk/english/employees/?pure=en/persons/487288/}
}
\begin{document}
\section{Introduction}
Over the last two decades, the increased availability of computing
power has led to a substantial increase in the availability and use of
Markov chain Monte Carlo (MCMC) methods for Bayesian estimation
\citep{Gilks:1998pt}. However, such methods have potential drawbacks
if used inappropriately, including difficulties in identifying
convergence \citep{Toft:2007vl, Brooks:1998ee} and the potential for
auto-correlation to decrease the effective sample size of the numerical
integration process \citep{Kass:1998fu}. Although writing MCMC
sampling algorithms such as the Metropolis-Hastings algorithm
\citep{Hastings:1970bi} is relatively straightforward, many users
employ software such as the Bayesian analysis Using Gibbs Sampling
(\proglang{BUGS}) software variants \pkg{WinBUGS} and \pkg{OpenBUGS}
\citep{Lunn00}. Just Another Gibbs Sampler
\citep[\pkg{JAGS};][]{Plummer2003} is a cross-platform alternative
with a direct interface to \proglang{R} using \pkg{rjags}
\citep{Plummer2013}, which can be easily extended with user-specified
modules supporting additional distributions and random number
generators \citep{Wabersich2014}. Each of these uses the \proglang{BUGS}
syntax to allow the user to define arbitrary models more easily, which
is attractive and attainable for researchers who are more familiar
with traditional modeling techniques. However, many of these less
experienced users may not be aware of the potential issues with MCMC
analysis, hence the prominent warning that ``MCMC sampling can be
dangerous'' in the \pkg{WinBUGS} user manual \citep{Lunn00}. Some of
this potential risk for inexperienced users can be reduced using a
wrapper for the model-fitting software that analyzes the model output
for common problems, such as failure to converge, parameter
auto-correlation and effective sample size, which may otherwise be
overlooked by the end user.
Bayesian statistical methods, such as those used by \proglang{BUGS}
and \pkg{JAGS}, also require prior belief to be incorporated into
the model. There are a number of different recommendations for an
appropriate choice of prior distribution under various different
circumstances, for example the half-Cauchy distribution has been
recommended as a reasonable choice for standard deviation parameters
within hierarchical models \citep{Gelman2006,Polson2011}, and
\citet{DuMouchel1994} gives an argument for the use of $\pi(\tau) =
\frac{s_0}{\left(s_0 + \tau\right)^2}$ as a prior for a variance
parameter $\tau$ in meta-analysis models. However, these are not
available as built-in distributions in \proglang{BUGS} or
\pkg{JAGS}.
This paper describes the \pkg{runjags} package \citep{runjags} for
\proglang{R} \citep{R2014} which can be used to automate MCMC fitting
and summarizing procedures for \pkg{JAGS} models and is available
from the Comprehensive \proglang{R} Archive Network (CRAN) at
\url{https://CRAN.R-project.org/package=runjags}. The functions are
designed to be user-friendly (particularly for those less experienced
with MCMC analysis), and provide a number of features to make the
recommended convergence and sample size checks more obvious to the end
user. The \pkg{runjags} package also provides additional
distributions to extend the core functionality of \pkg{JAGS},
including the half-Cauchy and DuMouchel distributions, as well as
functions implementing different types of simulation studies to assess
the performance of \pkg{JAGS} models. Section~\ref{sec:illustration}
gives a worked example of usage to assess the sensitivity of an
over-dispersed count observation model to various
minimally-informative prior distributions. Some prior familiarity
with the \proglang{BUGS} programming language and the underlying MCMC
algorithms is assumed. All code shown below is also included in an
\proglang{R} file in the supplementary material.
\section{Package functions}
\subsection{Preparation}
The core functionality of the \pkg{runjags} package allows a model
specified by the user to be run in \pkg{JAGS}, using the
\code{run.jags} function. The help file for this function gives an
overview of the core functionality of the \pkg{runjags} package and
provides links to other relevant functions. All functions require
installation of \pkg{JAGS}, which is an open source software package
available from \url{https://mcmc-jags.sourceforge.io/}.
Before running a model for the first time, it is advisable to check the installation of \pkg{JAGS} and set any desired global settings such as installation locations and warning message preferences using the \code{runjags.options} function. For example, the following will first test the \pkg{JAGS} installation, and then set function feedback from \pkg{runjags} and simulation updates from \pkg{JAGS} to be suppressed for future model runs in this \proglang{R} session:
%
\begin{Schunk}
\begin{Sinput}
R> testjags()
\end{Sinput}
\begin{Soutput}
You are using R version 3.3.0 (2016-05-03) on a unix machine, with the
X11 GUI
JAGS version 4.2.0 found successfully using the command
'/usr/local/bin/jags'
The rjags package is installed
\end{Soutput}
\begin{Sinput}
R> runjags.options(silent.runjags = TRUE, silent.jags = TRUE)
\end{Sinput}
\end{Schunk}
%
The help file for the \code{runjags.options} function gives a list of
other possible global options, and instructions on how to set these in
the \proglang{R} profile file for permanent use.
\subsection{Basic usage}
The \code{run.jags} function requires a valid model definition to the
\code{model} argument and a character string of monitored variables to
the \code{monitor} argument before a model can be run. The model can be
specified in an external text file, or as a character string within
\proglang{R}. The former is likely to be preferable for more complex
model formulations, but the latter eliminates the need for multiple
text files. Data will be necessary for most models, and it is highly
recommended to provide over-dispersed starting values for multiple
chains; the default settings give a warning if no initial values are
provided.
There are a number of ways to provide data and initial values,
depending on the preferences of the user. It is possible for the text
file containing the model to also contain data and initial value
``blocks'', in which case these will be automatically imported with the
model by \code{run.jags} and the number of chains is inferred from the
number of initial value lists found. This is also compatible with
standard \pkg{WinBUGS} or \pkg{OpenBUGS} text files, although the
addition of curly brackets is necessary to demarcate the data and
initial value blocks in the same way as for the model block. It is
also necessary to convert any \proglang{BUGS} arrays from row-major order
to column-major order, which is done automatically if the variables
are specified inside a list (as is the case for \proglang{BUGS}, but not
for \proglang{R}). To over-ride this setting within a specific data
or initial value block, the user can include \code{\#BUGSdata\#} to
ensure all arrays are converted from row- to column-major order,
\code{\#Rdata\#} to ensure none of the arrays are converted, and
\code{\#modeldata\#} to pass the data block directly to \pkg{JAGS} for
data transformation (see Section 7.0.4 of the \pkg{JAGS} user manual).
As a basic example, we can use the Salmonella example from
Chapter 6.5.2 of the \proglang{BUGS} book
\citep[\href{http://www.mrc-bsu.cam.ac.uk/software/bugs/the-bugs-project-the-bugs-book/bugs-book-examples/the-bugs-book-examples-chapter-6-6-5-2/}{\tt
http://www.mrc-bsu.cam.ac.uk/software/bugs/the-bugs-project-the-bugs-book/\linebreak bugs-book-examples/the-bugs-book-examples-chapter-6-6-5-2/}, with thanks to][for permission to reproduce their
model]{lunn2012bugs}. Simulation-specific options can be provided to
the \code{run.jags} function, which may include the required burn-in
period, sampling length and thinning interval. A basic model run with
a fixed burn-in period (default 4,000 iterations after 1,000 adaptive
iterations) and sampling period (default 10,000 iterations) can be
obtained as follows:
%
\begin{Schunk}
\begin{Sinput}
R> filestring <- "
+ The BUGS Book example Chapter 6.5.2
+ The following example has been modified only to include
+ curly brackets around the Data and Inits specifications
+
+ Poisson model...
+
+ model {
+ for (i in 1:6) {
+ for (j in 1:3) {
+ y[i,j] ~ dpois(mu[i])
+ }
+ log(mu[i]) <- alpha + beta * log(x[i] + 10) + gamma * x[i]
+ }
+ for (i in 1:6) {
+ y.pred[i] ~ dpois(mu[i])
+ }
+ alpha ~ dnorm(0, 0.0001)
+ beta ~ dnorm(0, 0.0001)
+ gamma ~ dnorm(0, 0.0001)
+ }
+
+ Data {
+ list(y = structure(.Data = c(15, 21, 29, 16, 18, 21, 16, 26, 33,
+ 27, 41, 60, 33, 38, 41, 20, 27, 42), .Dim = c(6, 3)),
+ x = c(0, 10, 33, 100, 333, 1000))
+ }
+
+ Inits {
+ list(alpha = 0, beta = 0, gamma = 0)
+ }
+ "
R> results <- run.jags(filestring, monitor = c("alpha", "beta", "gamma"))
\end{Sinput}
\begin{Soutput}
Warning message:
Convergence cannot be assessed with only 1 chain
\end{Soutput}
\end{Schunk}
%
A single chain was used for this model because only one set of initial
values was found in the example file, resulting in the warning message
regarding convergence assessment. The results of the simulation can
be examined using the default print method as follows:
%
\begin{Schunk}
\begin{Sinput}
R> results
\end{Sinput}
\begin{Soutput}
JAGS model summary statistics from 10000 samples (adapt+burnin = 5000):
Lower95 Median Upper95 Mean SD Mode
alpha 1.8018 2.1899 2.6135 2.1913 0.20439 2.1939
beta 0.20684 0.31538 0.41891 0.31414 0.053648 0.3157
gamma -0.0014678 -0.00099063 -0.00053861 -0.0009922 0.00023808 -0.001008
MCerr MC%ofSD SSeff AC.10 psrf
alpha 0.019523 9.6 110 0.78631 --
beta 0.0052936 9.9 103 0.81478 --
gamma 0.000019396 8.1 151 0.58568 --
Total time taken: 0.5 seconds
\end{Soutput}
\end{Schunk}
%
The results show similar inference to that provided by
\citet{lunn2012bugs}, although with additional information regarding
the effective sample size (\code{SSeff}), auto-correlation at a lag of
10 (\code{AC.10}), and the potential scale reduction factor
(\code{psrf}) of the Gelman-Rubin statistic \citep{Gelman1992} for
models with multiple chains (the latter is sometimes referred to as
`Rhat'). In this case, an insufficient number of
samples has been taken for this highly auto-correlated model (although
it is important to note that the auto-correlation is markedly reduced
if the `\code{glm}' module is loaded in \pkg{JAGS}). Displaying the
effective sample size with the summary information will alert the user
to the fact that additional steps should be taken before sensible
inference can be made.
The data can also be specified to \code{run.jags} using the
\code{data} argument, in which case it should take the format of a
named list, data frame, character string as produced by
\code{dump.format}, or a function (with no arguments) returning one of
these. Similarly, the initial values can be specified using the
\code{inits} argument as a list with length equal to the number of
chains, with each element specifying a named list, data frame or
character string for the initial values for that chain. The initial
values may also be specified as a function taking either no arguments
(as for the \code{data} argument) or one argument (specifying the
chain number), in which case an additional \code{n.chains} argument
will be required by \code{run.jags} to determine the number of chains
required.
\subsection{Alternative usage}
To facilitate a more streamlined function call within \proglang{R}, an
alternative method of specifying data and initial values is provided.
The model formulation may contain special inline comments including:
\code{\#data\#}, which indicates that the comma separated variable
names to the right of the statement are to be included in the
simulation as data, and \code{\#inits\#}, which indicates variables
for which initial values are to be provided. Any variables specified
by \code{\#data\#} and \code{\#inits\#} will be automatically
retrieved from a named list, data frame or environment passed to the
\code{data} and \code{inits} argument (or function returning one of
these), or from the global environment. Any variable names specified
in this way may also match a function returning an appropriate vector,
and in the case of initial values, this function may accept a single
argument indicating the chain for which the initial values are to be
used. Note that any variables specified by \code{\#data\#} or
\code{\#inits\#} will be ignored if a character string is provided to
the \code{data} or \code{inits} arguments, which may be useful for temporarily
over-riding the values specified in the model file. See the
\code{dump.format} function for a way to generate these. In addition
to \code{\#data\#} and \code{\#inits\#}, a number of optional inline
comments are supported as follows:
\begin{itemize}
\item \code{\#monitors\#} -- a comma-separated list of monitored
variables to use, which may include the special variables ``DIC''
\citep{Spiegelhalter:2002mq} and ``PED'' \citep{Plummer2008}, which
can be used to assess model fit;
\item \code{\#modules\#} -- a comma-separated list of any \pkg{JAGS} extension modules required, optionally also specifying the status (e.g., \code{\#modules\# glm on, dic on});
\item \code{\#factories\#} -- a comma-separated list of any \pkg{JAGS} factories and types required, optionally also specifying the status (e.g., \code{\#factories\# mix::TemperedMix sampler on});
\item \code{\#response\#} -- a single variable name specifying the response variable;
\item \code{\#residual\#} -- a single variable name specifying a variable that represents the residuals;
\item \code{\#fitted\#} -- a single variable name specifying a variable that represents the fitted value.
\end{itemize}
Each of these options can also be supplied directly to the relevant function call in \proglang{R}. An example of running a model using this style of model specification is as follows:
%
\begin{Schunk}
\begin{Sinput}
R> model <- "model {
+ for (i in 1:N) { #data# N
+ Y[i] ~ dnorm(true.y[i], precision) #data# Y
+ true.y[i] <- coef * X[i] + int #data# X
+ }
+ coef ~ dunif(-1000, 1000)
+ int ~ dunif(-1000, 1000)
+ precision ~ dexp(1)
+ #inits# coef, int, precision, .RNG.seed, .RNG.name
+ #monitor# coef, int, precision
+ }"
\end{Sinput}
\end{Schunk}
Simulate the data:
\begin{Schunk}
\begin{Sinput}
R> set.seed(1)
R> N <- 100
R> X <- seq(1, N, by = 1)
R> Y <- rnorm(N, 2 * X + 10, 1)
\end{Sinput}
\end{Schunk}
%
The following code specifies functions that return initial values
(including RNG seeds) for each chain. The use of \code{switch} within
these functions allows different initial values to be chosen for
chains one and two, ensuring that initial values are over-dispersed.
%
\begin{Schunk}
\begin{Sinput}
R> coef <- function(chain)
+ return(switch(chain, "1" = -10, "2" = 10))
R> int <- function(chain)
+ return(switch(chain, "1" = -10, "2" = 10))
R> precision <- function(chain)
+ return(switch(chain, "1" = 0.01, "2" = 100))
R> .RNG.seed <- function(chain)
+ return(switch(chain, "1" = 1, "2" = 2))
R> .RNG.name <- function(chain)
+ return(switch(chain, "1" = "base::Super-Duper",
+ "2" = "base::Wichmann-Hill"))
\end{Sinput}
\end{Schunk}
%
It is then possible to run the simulation specifying only the model
and the number of chains to use (the monitored variables, data and
initial values are specified in the model file and will be retrieved
form our \proglang{R} working environment):
%
\begin{Schunk}
\begin{Sinput}
R> results <- run.jags(model, n.chains = 2)
\end{Sinput}
\end{Schunk}
\label{sec:inline}
\subsection{Extending models}
The \code{autorun.jags} function can be used in the same way as
\code{run.jags}, but the burn-in period and sample length are
calculated automatically rather than being directly controlled by the
user. The \code{autorun.jags} function will continually extend a
simulation until convergence -- as assessed by the Gelman-Rubin
statistic \citep{Gelman1992} -- has been achieved for all monitored
variables, and will then extend the simulation further to compensate
for any observed auto-correlation. The automated assessment of
convergence should be verified graphically before making inference
from models fit to real data, but a fully automated analysis is useful
for simulated data and for reinforcing the importance of convergence
assessment for novice users. The following code will run the same
model as above, extending the model as necessary up to a maximum total
elapsed time of one hour:
%
\begin{Schunk}
\begin{Sinput}
R> results <- autorun.jags(model, n.chains = 2, max.time = "1hr")
\end{Sinput}
\end{Schunk}
%
Alternatively, an existing model may be extended by the user in order
to increase the sample size of the MCMC chains using either the
\code{extend.jags} or \code{autoextend.jags} function. For these
functions, the arguments \code{add.monitor}, \code{drop.monitor} and
\code{drop.chain} are provided in order to change the monitored variables
and number of chains being run. The \code{combine} argument controls
whether the old MCMC chains should be discarded, or combined with the
new chains. For example, the following code will manually extend the
existing simulation by 5000 iterations, and then extend the simulation
again with automatic control of convergence and sample size
diagnostics:
%
\begin{Schunk}
\begin{Sinput}
R> results <- extend.jags(results, sample = 5000)
R> results <- autoextend.jags(results)
\end{Sinput}
\end{Schunk}
%
In the second function call, the automated diagnostics run by
\code{autoextend.jags} determine that the simulation has converged and
already has an adequate sample size, so no additional samples are
taken. For more details on these functions including detailed
descriptions of the other arguments and additional examples, consult
the help pages for \code{run.jags} and \code{autorun.jags}.
Once a valid `\code{runjags}' class object has been obtained, the full
representation of the model, data and current status of the random
number generators can be saved to a file using the
\code{write.jagsfile} function. This allows a model to be run from
the last sampled values using the \code{run.jags} function at a later
time point, and it may also be instructive to use this function to
examine the format of a syntactically-valid and complete model file
that can be read directly using the \code{run.jags} function. It is
also possible to specify a value of 0 for the \code{sample} argument in the
original \code{run.jags} function call, and then subsequently use
\code{write.jagsfile} to produce a model file with the initial values
specified.
\subsection{Visualization methods}
The output of these functions is an object of class `\code{runjags}'.
This class is associated with a number of \proglang{S}3 methods, as
well as utility functions for combining multiple `\code{runjags}'
objects (\code{combine.jags}), and for conversion to and from objects
produced by the \pkg{rjags} package (\code{as.runjags} and
\code{as.jags}). Many of these allow a \code{vars} argument giving a
subset of monitored nodes (using partial matching), as well as a
\code{mutate} argument. This should specify a function (or a list with
first element a function and remaining elements arguments to this
function), and can be used to add new variables to the posterior
chains that are derived from the directly monitored variables in
\pkg{JAGS}. This allows the variables to be summarized or extracted
as part of the MCMC objects as if they had been calculated in
\pkg{JAGS}, but without the computational or storage overheads
associated with calculating them directly in \pkg{JAGS}. One possible
application for this is for pair-wise comparisons of different levels
within fixed effects using the supplied \code{contrasts.mcmc}
function.
The \code{print} method displays relevant overview information,
including summary statistics for monitored variables calculated and
stored by the \code{run.jags} function. The \code{summary} method
returns a summary table for the monitored variables, which is taken
from the stored values created by \code{run.jags} if available;
otherwise it will be recalculated during the function call.
Alternatively, summary statistics can be recalculated and stored in
the `\code{runjags}' object using the \code{add.summary} function.
There are a series of options available to these summary functions,
including \code{vars} and \code{mutate} as outlined above,
\code{confidence} which specifies a numeric vector of confidence
intervals to calculate, and \code{custom} which allows one or more
statistics calculated by a user-supplied function to be appended to
the summary statistics. Note that summary options may also be passed
to \code{run.jags} in order to control the summary statistics
calculated and appended to the `\code{runjags}' object.
The \code{plot} method produces a series of relevant plots for the
selected variables, including trace plots, empirical cumulative
distribution function plots, histograms, and a cross-correlation plot,
with additional options allowing auto-correlation plots and density
plots if desired. Further plot parameters can be specified using the
\code{col} and \code{separate.chains} arguments, as well as a named
list for each plot type which will be passed to the underlying
\pkg{lattice} functions \citep{lattice}. The primary intention with
these plots is to provide rapid access to commonly used convergence
diagnostics, and \code{plot} methods associated with `\code{mcmc}' or
`\code{mcmc.list}' objects may be more flexible and intuitive for
producing more specific graphical output from converged MCMC chains.
The \pkg{coda} package \citep{Plummer2006} provides such plotting
methods, as well as many of the underlying functions that calculate
the summaries given by \pkg{runjags}. A typical examination of a
simulation output (the default print method, and a plot output for
variable names partially matching the letter ``c'') could be obtained as
follows:
%
\begin{Schunk}
\begin{Sinput}
R> results
\end{Sinput}
\begin{Soutput}
JAGS model summary statistics from 30000 samples (chains = 2;
adapt+burnin = 5000):
Lower95 Median Upper95 Mean SD Mode MCerr MC%ofSD
coef 1.9922 1.9995 2.007 1.9996 0.0038115 1.9994 0.000073962 1.9
int 9.4388 9.8711 10.307 9.8692 0.22139 9.8768 0.0042942 1.9
precision 0.61225 0.82946 1.0685 0.83515 0.11731 0.826 0.0007013 0.6
SSeff AC.10 psrf
coef 2656 0.16954 1.002
int 2658 0.16942 1.0018
precision 27982 0.0020679 1
Total time taken: 4.3 seconds
\end{Soutput}
\begin{Sinput}
R> plot(results, vars = "c", layout = c(3, 3))
\end{Sinput}
\end{Schunk}
\begin{figure}[t!]
\centering
\includegraphics[width=\textwidth, trim = 0 20 0 30, clip]{traceplots.pdf}
\caption{A series of plots displayed by the \code{plot} method for the
`\code{runjags}' class, showing only parameters partially matched using the
letter ``c'' with plots shown in a 3 $\times$ 3 layout. Multiple
chains are shown using different colors as indicated by the key
plot.}
\end{figure}
The standard \code{plot} method presents the commonly required
information in an easily readable format (including model fit
statistics where available), but the same information can be returned
in the form of a numeric matrix using the \code{summary} method. To
extract additional information from the `\code{runjags}' object not
covered by these summary statistics, see the \code{extract} method.
\subsection{GLMM templates}
There are many available frameworks for fitting standard generalized
linear mixed models (GLMMs) in \proglang{R}, but new users to MCMC may
find that running relatively simple models in \pkg{JAGS} and comparing
the results to those obtained through other software packages allows
them to better understand the flexibility and syntax of \proglang{BUGS}
models. To this end, the \pkg{runjags} package provides a
\code{template.jags} function which generates model specification
files based on a formula syntax similar to that employed by the
well-known \pkg{lme4} package \citep{Bates2014, lme4}. After
generating the template model, the user is encouraged to examine the
model file and make whatever changes are necessary before running the
model using \code{run.jags}. For example, a basic generalized linear
model (from the help file for \code{glm}) can be compared to the
output of \pkg{JAGS} as follows:
%
\begin{Schunk}
\begin{Sinput}
R> counts <- c(18, 17, 15, 20, 10, 20, 25, 13, 12)
R> outcome <- gl(3, 1, 9)
R> treatment <- gl(3, 3)
R> d.AD <- data.frame(treatment, outcome, counts)
R> glm.D93 <- glm(counts ~ outcome + treatment, family = poisson())
R> template.jags(counts ~ outcome + treatment, data = d.AD,
+ family = "poisson")
\end{Sinput}
\begin{Soutput}
Your model template was created at "JAGSmodel.txt" - it is highly
advisable to examine the model syntax to be sure it is as intended
You can then run the model using run.jags("JAGSmodel.txt")
\end{Soutput}
\begin{Sinput}
R> jags.D93 <- run.jags("JAGSmodel.txt")
\end{Sinput}
\end{Schunk}
%
The results of these comparisons are not displayed here, but show how
the same inference is presented slightly differently in a Bayesian
framework. The \code{template.jags} function supports Gaussian,
(zero-inflated) binomial, (zero-inflated) Poisson and (zero-inflated)
negative binomial distributions, as well as linear and fixed effects,
2-way interactions and random intercept terms specified using the same
syntax as used in \pkg{lme4}. Additional distributions and link functions
can be introduced by manually editing the template model file. All
necessary data, initial values, monitored variables and modules are
saved to the model file using the previously described comment syntax,
and the template function also saves information about the response
variable, fitted estimates and residuals to the model file, allowing
\code{residuals} and \code{fitted} methods to be used with the objects
returned by \code{run.jags}.
\subsection[JAGS module]{\pkg{JAGS} module}
\label{sec:pareto}
In addition to the \proglang{R} code used to facilitate running \pkg{JAGS} models and summarizing results, the \pkg{runjags} package also provides a modular extension to the \pkg{JAGS} language, providing additional distributions. The module can be loaded using the following command:
%
\begin{Schunk}
\begin{Sinput}
R> load.runjagsmodule()
\end{Sinput}
\begin{Soutput}
module runjags loaded
\end{Soutput}
\end{Schunk}
%
This makes the module available to any \pkg{JAGS} model, including
those run using the \pkg{rjags} package. The available distributions
extend the Pareto Type I distribution provided within \pkg{JAGS} to
Pareto Types II, III and IV, as well as providing the generalized
Pareto distribution, the Lomax distribution (a special case of the
Pareto Type II distribution with $\mu=0$), and two distributions
advocated for use as ``minimally-informative'' priors for variance
parameters: the DuMouchel distribution \citep{DuMouchel1994}, and the
half-Cauchy distribution \citep{Gelman2006}. The usage, probability
density function (PDF) and lower bound for the support of each of the
distributions provided by the module are shown in
Table~\ref{tab:pareto}, and an example of how to use the distributions
in this module is given in Section~\ref{sec:illustration}.
One limitation of the module provided within \pkg{runjags} is that it
is only made available for the `\code{rjags}' and `\code{rjparallel}'
methods when loaded within \proglang{R}. However, a standalone
\pkg{JAGS} module containing the same functions for use with any
\pkg{JAGS} installation (independently of \proglang{R}) is available
from \url{http://runjags.sourceforge.net/}. This module is named
`\code{paretoprior}' to avoid naming conflicts with the internal
\pkg{runjags} module, and should install on a variety of platforms
using the standard `\code{./configure}', `\code{make}', `\code{make
install}' convention. Binary installers are also provided for some
platforms.
\begin{table}[t!]
{
\centering
\begin{tabular}{ l l c c }
\addlinespace
\hline
Name & Usage in \pkg{JAGS} & Density & Lower \\
\hline
Pareto I$^1$ & \parbox[t]{5.5cm}{\code{ dpar1(alpha, sigma)\\ $\alpha>0, \, \sigma>0$}} & \parbox[c][1.5cm]{5.5cm}{$$\alpha \, \sigma^\alpha \, x^{-\left(\alpha + 1\right)}$$ }& $\sigma$ \\
\addlinespace
Pareto II & \parbox[t]{5.5cm}{\code{ dpar2(alpha, sigma, mu) \\ $\alpha>0, \, \sigma>0$ }} & \parbox[c][1.5cm]{5.5cm}{$$\frac{\alpha}{\sigma} \left(\frac{\sigma+x-\mu}{\sigma}\right)^{-\left(\alpha + 1\right)}$$} & $\mu$ \\
\addlinespace
Pareto III & \parbox[t]{5.5cm}{\code{ dpar3(sigma, mu, gamma) \\ $\sigma>0, \, \gamma>0$ }} & \parbox[c][1.5cm]{5.5cm}{$$ \frac{\left(\frac{x-\mu}{\sigma}\right)^{\frac{1}{\gamma} -1} {\left(\frac{x-\mu}{\sigma}^{\frac{1}{\gamma}} + 1\right)}^{-2}}{\gamma \, \sigma} $$} & $\mu$ \\
\addlinespace
Pareto IV & \parbox[t]{5.5cm}{\code{ dpar4(alpha, sigma, mu, gamma) \\ $\alpha>0, \, \sigma>0, \, \gamma>0$ }} & \parbox[c][1.5cm]{5.5cm}{$$ \frac{\alpha \left(\frac{x-\mu}{\sigma}\right)^{\frac{1}{\gamma} -1} {\left(\frac{x-\mu}{\sigma}^{\frac{1}{\gamma}} + 1\right)}^{-\left(\alpha+1\right)}}{\gamma \, \sigma} $$} & $\mu$ \\
\addlinespace
Lomax$^2$ & \parbox[t]{5.5cm}{\code{ dlomax(alpha, sigma) \\ $\alpha>0, \, \sigma>0$ }} & \parbox[c][1.5cm]{5.5cm}{$$ \frac{\alpha}{\sigma} \, \left(1+\frac{x}{\sigma}\right)^{-\left(\alpha+1\right)} $$} & $0$ \\
\addlinespace
Gen. Par. & \parbox[t]{5.5cm}{\code{ dgenpar(sigma, mu, xi) \\ $\sigma>0$ }} & \parbox[c][1.5cm]{5.5cm}{$$\frac{1}{\sigma} \, \left(1+\xi \frac{x-\mu}\sigma\right)^{-\left(\frac{1}{\xi}+1\right)} $$} & $\mu$$^3$ \\
& & \parbox[c][1.5cm]{5.5cm}{For $\xi=0$:\;\; $\frac{1}{\sigma} \, e^{\frac{-\left(x-\mu\right)}{\sigma}} $} & \\
\addlinespace
DuMouchel & \parbox[t]{5.5cm}{\code{ dmouch(sigma) \\ $\sigma>0$ }} & \parbox[c][1.5cm]{5.5cm}{$$\frac{\sigma}{\left(x + \sigma\right)^2} $$} & $0$ \\
Half-Cauchy & \parbox[t]{5.5cm}{\code{ dhalfcauchy(sigma) \\ $\sigma>0$ }} & \parbox[c][1.5cm]{5.5cm}{$$\frac{2 \sigma}{\pi \left(x^2 + \sigma^2\right)} $$} & $0$ \\
\hline
\addlinespace
\end{tabular}
}
\caption{Distributions provided by the \pkg{JAGS} module included with the \pkg{runjags} package. The name, \pkg{JAGS} code with parameterization, PDF and lower bound of the distributions are shown. All distributions have an upper bound of $\infty$ unless otherwise stated. \newline
$^1$This is equivalent to the \code{dpar(alpha, c)} distribution and provided for naming consistency. \newline
$^2$This is referred to as the ``2nd kind Pareto'' distribution by \citet{D2009}; an alternative form for the PDF of this distribution is given by: $\frac{\alpha \, \sigma^\alpha}{\left(x + \sigma\right)^{\alpha+1}} $. \newline
$^3$The Generalized Pareto distribution also has an upper bound of $x \le \mu - \frac{\sigma}{\xi}$ for $\xi < 0$. }
\label{tab:pareto}
\end{table}
\subsection{Method options}
There are a number of different methods for calling \pkg{JAGS} from
within \proglang{R} using \pkg{runjags}, which can be controlled using
the \code{method} argument or by changing the global option using the
\code{runjags.options} function. The main difference between these is
that some allow multiple chains to be run in parallel using separate
\pkg{JAGS} models, with automatic pseudo-random number generation
handled by \pkg{runjags} where necessary. The \code{"interruptible"},
\code{"rjags"}, \code{"parallel"} or \code{"bgparallel"} methods are
recommended for most situations, but all possible methods and their
advantages and disadvantages are summarized in Table~\ref{tab:methods}.
Note that a pre-existing cluster created using the \pkg{parallel}
package can be used by specifying a \code{cl} argument, and a maximum
number of parallel simulations for these methods can optionally be
specified using a \code{n.sims} argument to the main function call (the
default will use a separate simulation per chain, but it is possible
to specify fewer simulations than chains). The model fit statistics
are not available with parallel methods because multiple chains within
the same model are required for calculation of DIC and PED, but these
can be obtained using the \code{extract} method which will extend the
simulation using a single simulation. The adaptation phase is always
explicitly controlled to allow MCMC simulations with the same
pseudo-random number seed to be reproducible regardless of the method
used to call \pkg{JAGS}.
\begin{table}[t!]
\centering
\begin{tabular}{ l p{6.5cm} p{5cm} }
\addlinespace
\hline
{Method name} & {Description} & {Method options} \\
\hline
\addlinespace
\code{"interruptible"}$^1$ & \pkg{JAGS} called using a shell, with output monitored and displayed within \proglang{R}. & -- \\
\addlinespace
\code{"rjags"}$^{1, 2}$ & \pkg{JAGS} called using the \pkg{rjags} package. & \code{by} and \code{progress.bar}: as for the \pkg{rjags} package. \\
\addlinespace
\code{"background"}$^{1~*}$ & \pkg{JAGS} called as a background process, with the \proglang{R} prompt returned to the user. & -- \\
\addlinespace
\code{"simple"}$^1$ & \pkg{JAGS} called directly using a shell. & -- \\
\addlinespace
\hline
\addlinespace
\code{"parallel"}$^{3}$ & Multiple \pkg{JAGS} instances called using separate shells to allow chain parallelization. & \code{n.sims}: the number of parallel simulations. \\
\addlinespace
\code{"bgparallel"}$^{3~*}$ & Multiple \pkg{JAGS} instances called using separate background processes to allow chain parallelization. & \code{n.sims}: the number of parallel simulations. \\
\addlinespace
\code{"rjparallel"}$^{3, 4}$ & Multiple \pkg{rjags} models run within \proglang{R} using a \pkg{parallel} cluster. & \code{cl}: a pre-created cluster to be used, and \code{n.sims}: the number of parallel simulations. \\
\addlinespace
\code{"snow"}$^1$ & Multiple \pkg{JAGS} instances called using separate shells set up using a \pkg{parallel} cluster. & \code{cl} and \code{n.sims}: as above, and \code{remote.jags}: the \pkg{JAGS} path on the cluster nodes. \\
\addlinespace
\hline
\addlinespace
\end{tabular}
\caption{Methods provided by the \pkg{runjags} package to run simulations in \pkg{JAGS}.
Availability of \pkg{JAGS} modules is as follows:\newline
$^1$Installed in \pkg{JAGS}. \newline
$^2$Loadable in the \proglang{R} session. \newline
$^3$Installed in \pkg{JAGS} (except DIC). \newline
$^4$Loadable in \proglang{R} code run remotely on the cluster nodes (except DIC).\newline
$^*$These methods are not compatible with \code{autorun.jags} and \code{autoextend.jags}.
}
\label{tab:methods}
\end{table}
The two background methods do not return a completed simulation, but
instead create a folder in the working environment where the
simulation results will be written once the \pkg{JAGS} process has
completed. For example, the following code will allow a \pkg{JAGS}
simulation to be run in the background using two processors in
parallel, and saving the results in a folder called `\code{mysimulation}' in
the current working directory:
%
\begin{Schunk}
\begin{Sinput}
R> info <- run.jags(model, n.chains = 2, method = "bgparallel",
+ keep.jags.files = "mysimulation", n.sims = 2)
\end{Sinput}
\begin{Soutput}
Starting the simulations in the background...
The JAGS processes are now running in the background
\end{Soutput}
\end{Schunk}
%
This returns the control of the terminal to the user, who can then
carry on working in \proglang{R} while waiting for the simulation to complete.
The default behavior on completion of the simulations is to alert the
user by emitting a beep from the speakers, but configuration using
\code{runjags.options} allows a shell script file to be executed
instead. The \code{info} variable in this code contains the name and
directory of the simulation, which is given to the user if the object
is printed. The results can be retrieved using either the folder name
or the variable returned by the function that started the simulation:
%
\begin{Schunk}
\begin{Sinput}
R> background.results <- results.jags("mysimulation")
\end{Sinput}
\end{Schunk}
If the simulation has not yet completed, the \code{results.jags}
function will display the \pkg{JAGS} output so that the user can gauge
how much longer the simulation will take. Further options for the
\code{results.jags} function include \code{recover.chains} which
allows the results of successful simulations to be read even if other
parallel simulations did not produce output, and \code{read.monitor}
which allows only a chosen subset of the monitored variables to be
read from the MCMC output. For all methods except \code{"rjags"} and
\code{"rjparallel"}, any calls to \code{run.jags} where the
\code{keep.jags.files} argument is specified will result in a folder
being created in the working directory that can be reimported using
\code{results.jags}. Any failed simulations created are also kept
using the same mechanism, and a message is displayed detailing how the
user can attempt to recover these simulations. These failed
simulation folders are automatically cleaned up when the \proglang{R}
session is terminated. The \code{failed.jags} function returns any
output captured from \pkg{JAGS} in such cases, and is helpful to debug
model code.
\subsection{Simulation studies}
One of the principle motivations behind the development of the
\pkg{runjags} package is to automate the analysis of simulated
data sets for the purposes of model validation. A common motivation
for this type of analysis is a drop-$k$ validation study, also known as
a leave-one-out cross-validation where $k=1$. This procedure re-fits
the same model to a single data set multiple times, with one or more of
the observed data points removed from each re-fit of the model. This
can either be a randomly selected group of a fixed number ``k'' of data
points, or each individual data point in turn. The goal is to evaluate
the ability of the model to predict each observation from the
explanatory variables, so that any unusual observations can be
identified. While it is possible to repeatedly use the
\code{autorun.jags} function to analyze multiple data sets, the higher
level \code{run.jags.study} and \code{drop.k} functions are provided
to automate much of this process. Large simulation studies are likely
to be computationally intensive, but are ideal candidates for
parallelization. For this reason, parallel computation is built
directly into these functions using the \pkg{parallel} package. This
can be used to parallelize the simulation locally, or to run the
simulation on any cluster set up using the \pkg{snow} package
\citep{Tierney2013}. This allows for the maximization of the
available computing power without requiring the end user to write any
additional code, and includes an initial check to ensure that the
model compiles and runs locally before beginning the parallelized
study.
A drop-$k$ study is implemented in \pkg{runjags} using the
\code{drop.k} function as follows. The `\code{runjags}' class object on which
the drop-$k$ analysis will be performed must first be obtained using
the \code{run.jags} function. Here, we will use the simple linear
regression model obtained in Section~\ref{sec:inline}, with the result
of \code{run.jags} contained in the variable \code{results}. The
\code{drop.k} function takes arguments \code{drop} (indicating the
data variables to remove between simulations), and \code{k}
(indicating the number of data points to drop for each simulation). In
this case, a drop-1 study is run with the number of simulations equal
to the number of data points. All individual simulations are run
using the underlying \code{autorun.jags} function; additional
arguments for \code{autorun.jags} can be passed through \code{drop.k}
as required. The initial values for each simulation are taken from
the parent simulation, including the observed values of the removed
data points to ensure that the model will compile. The drop-1 study
is run and the results displayed using the following syntax (limited
to the first five data-points for brevity):
%
\begin{Schunk}
\begin{Sinput}
R> assessment <- drop.k(results, drop = "Y[1:5]", k = 1)
R> assessment
\end{Sinput}
\begin{Soutput}
Values obtained from a drop-k study with a total of 5 simulations:
Target Median Mean Lower95%CI Upper95%CI Range95%CI Within95%CI
Y[1] 11.544 11.885 11.879 9.6268 14.055 4.4282 1
Y[2] 13.051 13.914 13.906 11.549 15.994 4.4445 1
Y[3] 15.828 15.888 15.878 13.559 18.059 4.5009 1
Y[4] 18.784 17.832 17.83 15.663 20.12 4.4572 1
Y[5] 21.398 19.821 19.821 17.657 22.086 4.4285 1
AutoCorr(Lag10)
Y[1] 0.010115
Y[2] 0.014556
Y[3] 0.0015873
Y[4] 0.0047826
Y[5] 0.0064177
Average time taken: 2.6 seconds (range: 2.5 seconds - 2.7 seconds)
Average adapt+burnin required: 5000 (range: 5000 - 5000)
Average samples required: 10506 (range: 10000 - 11292)
\end{Soutput}
\end{Schunk}
%
The results show the 95\% confidence interval (CI) for each data point
obtained from the corresponding simulation where this data point was
removed, which in this case indicates that the first five data-points
were predicted reasonably well. For drop-$k$ cross-validation with
``k'' greater than 1, the indicated number of data points will be
randomly removed from each simulation and the average values for the
corresponding summary statistics from each data point will be shown.
In this case, the argument \code{simulations} must also be provided.
Additional arguments to \code{autorun.jags} can also be provided to
the \code{drop.k} function. For example, the following syntax will
run 100 simulations with a random selection of 2 of the 5 first five
data-points removed from each:
%
\begin{Schunk}
\begin{Sinput}
R> assessment <- drop.k(results, drop = "Y[1:5]", k = 2, simulations = 100,
+ method = "simple", psrf.target = 1.1)
R> assessment
\end{Sinput}
\begin{Soutput}
Average values obtained from a drop-k study with a total of 100 simulations:
Target Av.Median Av.Mean Av.Lower95%CI Av.Upper95%CI Av.Range95%CI
Y[1] 11.544 11.879 11.874 9.6118 14.079 4.4667
Y[2] 13.051 13.894 13.899 11.667 16.13 4.4629
Y[3] 15.828 15.868 15.871 13.636 18.088 4.4516
Y[4] 18.784 17.847 17.853 15.627 20.036 4.4082
Y[5] 21.398 19.822 19.827 17.656 22.071 4.4152
Prop.Within95%CI Av.AutoCorr(Lag10) Simulations
Y[1] 1 0.010223 41
Y[2] 1 0.0066466 39
Y[3] 1 0.0032953 42
Y[4] 1 0.012031 38
Y[5] 1 0.0076121 40
Average time taken: 6.2 seconds (range: 3.4 seconds - 7.4 seconds)
Average adapt+burnin required: 5000 (range: 5000 - 5000)
Average samples required: 10645 (range: 10000 - 12182)
\end{Soutput}
\end{Schunk}
%
In the latter case, inference was made on each data point in several
different data sets, so the results present the mean values of each
summary statistic obtained from the multiple simulations.
The \code{drop.k} function is a wrapper for the \code{run.jags.study}
function, which can be used to perform various different types of
simulation studies. This function takes the following arguments: the
number of data sets to analyze, the model to use, a function to produce
data that will be provided to each simulation, and a named list of
``target'' variables with true values representing parameters to be
monitored and used to summarize the output of the simulation. Inline
\code{\#monitor\#} statements can be used as with \code{run.jags}, and
any target variables are also automatically monitored. Any variables
specified using the inline \code{\#data\#} statement will be retrieved
from the working environment as usual and will be common to all
simulations -- data which is intended to change between simulations
must therefore be provided using the \code{datafunction} argument instead.
Initial variables can be specified using \code{\#inits\#} in the model
file, but it is also necessary to pass a character string of all
variable names required to the \code{export.cluster} argument to ensure
these variables are visible on the cluster nodes. It may be
preferable to specify initial values as a function, to which the data
will be made available by \code{run.jags} at run time (this may be
required in cases where the choice of appropriate initial values
depends on the values in the data). An illustration of the
\code{run.jags.study} function is provided in
Section~\ref{sec:illustration}.
\section{Illustration of usage with a simulation study}
\label{sec:illustration}
Here we will consider a worked example of a simulation study analysis
using \pkg{runjags}, in order to assess the performance of two
equivalent model formulations with two different
``minimally-informative'' priors. The application is an over-dispersed
count model, the use of which is widespread in many biological fields
\citep{Bolker2009}, including parasitology
\citep{Wilson:1996vm,Wilson199733,Shaw:1998ty}, where Bayesian methods
of analysis have been shown to provide more robust inference than
traditional methods \citep{Denwood:2008rz,Denwood:2010fk}.
\subsection{Model formulation and assessment}
\label{sec:origmodel}
The gamma distribution is parameterized in \pkg{JAGS} and \proglang{BUGS} by the \code{shape} ($\alpha$) and \code{rate} ($\beta$) parameters, with the expectation given by $\frac{\alpha}{\beta}$ and variance given by $\frac{\alpha}{\beta^2}$. This distribution can be used to describe underlying variability in a Poisson observation, representing an unknown amount of over-dispersion between observations. In this situation the extra-Poisson coefficient of variation $cv$ is a useful measure of the variability of the underlying gamma distribution, and is a simple function of the \code{shape} parameter: $cv=\sqrt{\frac{1}{\alpha}}$
A candidate \pkg{JAGS} Model A (using inline data and monitor statements to be detected by \pkg{runjags}) is as follows:
%
\begin{Schunk}
\begin{Sinput}
R> ModelA <- "model {
+ for (i in 1:N) {
+ Count[i] ~ dpois(lambda[i])
+ lambda[i] ~ dgamma(shape, rate)
+ }
+ shape ~ dmouch(1)
+ mean ~ dmouch(1)
+ rate <- shape / mean
+
+ #data# N
+ #modules# runjags
+ #monitor# mean, shape
+ }"
\end{Sinput}
\end{Schunk}
%
This model allows each observed \code{Count} to follow a Poisson
distribution with \code{lambda} drawn from a gamma distribution with
\code{shape} parameter to be estimated, and \code{rate} parameter
calculated from the \code{shape} parameter and the \code{mean} of the
distribution, which is also to be estimated. The prior distribution
used for the \code{mean} and \code{shape} parameters is the DuMouchel
prior distribution as shown in Table~\ref{tab:pareto} -- this
distribution is provided by the \pkg{runjags} extension module which
can be loaded using the \code{\#modules\#} tag. Here we use the same
minimally-informative prior distribution for both \code{shape} and
\code{mean} parameters. The \code{\#data\#} statement is used to
include \code{N} as data that does not change between simulations.
The \code{Count} variable is also observed, but will vary between
simulations so it is not retrieved from \proglang{R} memory using \code{\#data\#}.
An alternative formulation of this same model could be provided using
a negative binomial distribution rather than a gamma mixture of
Poisson distributions, as represented in Model B:
%
\begin{Schunk}
\begin{Sinput}
R> ModelB <- "model {
+ for (i in 1:N) {
+ Count[i] ~ dnegbin(prob, shape)
+ }
+
+ shape ~ dmouch(1)
+ mean ~ dmouch(1)
+ prob <- shape / (shape + mean)
+
+ #data# N
+ #modules# runjags
+ #monitor# mean, shape
+ }"
\end{Sinput}
\end{Schunk}
%
In this model, the same priors are placed on the parameters
\code{shape} and \code{mean}, but the negative binomial distribution
is parameterized by a probability \code{p} in place of the parameter
\code{mean}. However, the gamma-Poisson and negative binomial
distributions are equivalent (see Appendix~\ref{sec:derivation}), and
these models share the same prior distributions for the two parameters
of interest. The two might therefore be expected to give equivalent
inference.
The posterior coverage and auto-correlation of these models can be
assessed using simulation studies, with data generated from a
distribution with a mean of 2, $cv$ of 1.1, and sample size of 20.
These values are chosen to exaggerate any model performance issues by
providing a comparatively small data set with a large number of zero
observations, and are similar to those typically found in veterinary
parasitological data sets \citep{denwood2010thesis}. The two
parameters of interest are the \code{mean} parameter which is directly
monitored in the model, and the $cv$ parameter which is a function of
the monitored \code{shape} parameter. Rather than calculate the $cv$
parameter in \pkg{JAGS}, this can be calculated more efficiently in
\proglang{R} using a mutate function:
%
\begin{Schunk}
\begin{Sinput}
R> getcv <- function(x)
+ return(list(cv = sqrt(1 / x[, "shape"])))
\end{Sinput}
\end{Schunk}
%
The model performance assessment can be automated using
\code{run.jags.study} by creating a function to return a pre-generated
simulated data set for each simulation:
%
\begin{Schunk}
\begin{Sinput}
R> N <- 20
R> S <- 1000
R> truemean <- 2
R> truecv <- 1.1
R> trueshape <- 1 / truecv^2
R> truerate <- trueshape / truemean
R> set.seed(1)
R> alldata <- lapply(1:S, function(x) {
+ return(rpois(N, rgamma(N, trueshape, rate = truerate)))
+ })
R> datafunction <- function(i) return(list(Count = alldata[[i]]))
\end{Sinput}
\end{Schunk}
%
In this case we specify the initial values as a function, illustrating
the potential to make use of the stochastically-generated data while
creating the initial values within the function:
%
\begin{Schunk}
\begin{Sinput}
R> initsfunction <- function(chain) {
+ stopifnot(data$N == 20)
+ stopifnot(chain %in% c(1, 2))
+ shape <- c(0.1, 10)[chain]
+ mean <- c(10, 0.1)[chain]
+ .RNG.seed <- c(1, 2)[chain]
+ .RNG.name <- c("base::Super-Duper", "base::Wichmann-Hill")[chain]
+ return(list(shape = shape, mean = mean, .RNG.seed = .RNG.seed,
+ .RNG.name = .RNG.name))
+ }
\end{Sinput}
\end{Schunk}
%
Finally, a \pkg{parallel} cluster with 10 nodes is set up on the local
machine, before the two simulation studies are run on this cluster
using the same data. The \code{run.jags.study} function will check
each of the models locally using a single randomly chosen data set to
ensure that the model is valid before it is passed to the cluster:
%
\begin{Schunk}
\begin{Sinput}
R> library("parallel")
R> cl <- makeCluster(10)
R> resultsA <- run.jags.study(S, ModelA, datafunction,
+ targets = list(mean = truemean, cv = truecv), cl = cl,
+ inits = initsfunction, n.chains = 2, mutate = getcv)
R> resultsB <- run.jags.study(S, ModelB, datafunction,
+ targets = list(mean = truemean, cv = truecv), cl = cl,
+ inits = initsfunction, n.chains = 2, mutate = getcv)
\end{Sinput}
\end{Schunk}
%
Each function call returns an object of class `\code{runjagsstudy}', with a
default \code{print} method that summarizes the results as for \code{drop.k}:
%
\begin{Schunk}
\begin{Sinput}
R> resultsA
\end{Sinput}
\begin{Soutput}
Average values obtained from a JAGS study with a total of 1000 simulations:
Target Av.Median Av.Mean Av.Lower95%CI Av.Upper95%CI Av.Range95%CI
mean 2 1.9689 2.0756 0.99815 3.354 2.3559
cv 1.1 1.0627 1.0909 0.52457 1.6983 1.1737
Prop.Within95%CI Av.AutoCorr(Lag10) Simulations
mean 0.93 0.043658 1000
cv 0.926 0.19093 1000
Average time taken: 5.2 seconds (range: 2.3 seconds - 12.4 seconds)
Average adapt+burnin required: 5550 (range: 5000 - 27000)
Average samples required: 10059 (range: 10000 - 21465)
\end{Soutput}
\begin{Sinput}
R> resultsB
\end{Sinput}
\begin{Soutput}
Average values obtained from a JAGS study with a total of 1000 simulations:
Target Av.Median Av.Mean Av.Lower95%CI Av.Upper95%CI Av.Range95%CI
mean 2 1.9632 2.069 1.0007 3.3294 2.3287
cv 1.1 1.0625 1.0891 0.52306 1.6905 1.1675
Prop.Within95%CI Av.AutoCorr(Lag10) Simulations
mean 0.925 0.015783 1000
cv 0.931 0.040415 1000
Average time taken: 4.7 seconds (range: 1.9 seconds - 10.1 seconds)
Average adapt+burnin required: 5099 (range: 5000 - 16000)
Average samples required: 10000 (range: 10000 - 10000)
\end{Soutput}
\end{Schunk}
%
The inference made from the two models indicates that they are
generally similar, except that the auto-correlation for both
parameters is reduced for Model B, meaning that on average fewer
samples were required for this model. As would be expected, the 95\%
confidence intervals for both parameters identified the true value
approximately 95\% of the time.
\subsection{Sensitivity to prior distributions}
\label{sec:priors}
\begin{table}[t!]
\centering
\begin{tabular}{ l l c c c c c c c}
\hline
Parameter & Priors: mean | shape & Mean & CI Range & Within CI & AC10 & Simulations \\
\hline
mean & dmouch | dmouch & 2.069 & 2.329 & 0.925 & 0.016 & 1000 \\
mean & dmouch | dgamma & 2.072 & 2.340 & 0.916 & 0.021 & 1000 \\
mean & dgamma | dmouch & 2.150 & 2.527 & 0.934 & 0.032 & 1000 \\
mean & dgamma | dgamma & 2.156 & 2.556 & 0.922 & 0.041 & 1000 \\
\addlinespace
\hline
\addlinespace
$cv$ & dmouch | dmouch & 1.089 & 1.167 & 0.931 & 0.040 & 1000 \\
$cv$ & dmouch | dgamma & 1.068 & 1.270 & 0.882 & 0.089 & 1000 \\
$cv$ & dgamma | dmouch & 1.093 & 1.173 & 0.933 & 0.041 & 1000 \\
$cv$ & dgamma | dgamma & 1.073 & 1.277 & 0.881 & 0.090 & 1000 \\
\hline
\end{tabular}
\caption{Average values for the inference on the mean parameter (true value 2) and $cv$ parameter (true value 1.1) obtained from a negative binomial MCMC model formulation using DuMouchel and gamma priors for the mean and shape parameters.}
\label{tab:priors}
\end{table}
The ability to incorporate prior information is an advantage of
Bayesian methods, but there is often a variety of potential
distributions that could be equally justifiable in a given situation.
The choice between these possibilities is known to affect the shape of
the posterior in some situations \citep{Lele2009}, particularly when
the information in the data is relatively sparse. In particular,
there are various different minimally-informative priors advocated for
use with variance parameters in hierarchical models, including the
\code{Gamma(0.001, 0.001)} distribution which is characterized by a
mean of one and a very large variance. The sensitivity of a model to
the choice of priors between this gamma prior and the DuMouchel prior
can be evaluated using the \code{run.jags.study} function, with a
total of four candidate sets of priors (using each combination of
DuMouchel and gamma distributions for the \code{mean} and \code{shape}
parameters). These were applied to the same 1,000 simulated data sets
using Model B and very similar \proglang{R} code to that given above.
The results of these four simulation studies are shown in
Table~\ref{tab:priors}. There are small but noticeable differences
between the inference made for both parameters using these prior
distributions. The bias and auto-correlation are both approximately
doubled for the \code{mean} parameter between DuMouchel and gamma
priors, and more substantial changes in bias and auto-correlation are
seen between priors for the $cv$ parameter. In addition, the 95\%
confidence intervals for the $cv$ parameter have less than 90\%
coverage when using the gamma prior, despite a slightly larger average
range of these confidence intervals relative to the DuMouchel prior.
\subsection{Discussion}
The results presented here demonstrate the utility of simulation
studies facilitated by the \pkg{runjags} package to evaluate the
relative performance of alternative model formulations and the effect
of prior distribution choices. In this case, the DuMouchel prior
out-performed the more standard gamma prior, and it also possesses
properties that are theoretically desirable for a
minimally-informative distribution, such as invariance to inverse
transformation, infinite variance and a mode of zero.
\citet{DuMouchel1994} proposed this prior for use with variance
parameters in hierarchical models, but it has also been used in
situations outside the meta-analysis application for which it was
originally devised \citep[see for example][]{Phillips2010,
Conti2011,Yin2013}. \citet{Christiansen1997} also used the same
distribution as a prior for a hierarchical regression model, and
\citet{Daniels1999} uses a uniform shrinkage prior which is equivalent
to the DuMouchel distribution. Although this connection is not stated
directly by \citet{DuMouchel1994}, the distribution is equivalent to a
Lomax distribution with $\tau=x$, $s_0=\sigma$ and $\alpha=1$, and
therefore to a Pareto type II distribution with $\tau=x$,
$s_0=\sigma$, $\alpha=1$ and $\mu=0$ (Table~\ref{tab:pareto}). The
choice of $\sigma$ dictates the median -- a value of 1 is advocated
since this also ensures invariance to the inverse transformation of
$\tau$, so this prior is equivalent in terms of variance and
precision. The half-Cauchy distribution has a similar form to the
DuMouchel distribution, and has also been suggested for use as a prior
for variance parameters \citep{Gelman2006,Polson2011}.
Although it is also possible to extend other variants of \proglang{BUGS},
\pkg{JAGS} is fully open source and written in \proglang{C++}, making
extension modules such as the one provided by \pkg{runjags} much
easier to implement. A very useful tutorial on writing and installing
a standalone \pkg{JAGS} module is provided by \cite{Wabersich2014},
but it is arguably more straightforward to implement a shared
\pkg{JAGS} library inside an \proglang{R} package. The configure
script provided inside the \pkg{runjags} package sets up the necessary
environmental variables for compilation on any platform, and can be
used as a template for creating additional extension modules within
\proglang{R} packages.
\section{Summary}
There are several advantages to using MCMC, but also some potential
disadvantages associated with failure to identify poor convergence and
high Monte Carlo error. The \pkg{runjags} package attempts to
partially safeguard against some of these difficulties by calculating
and automatically reporting convergence and sample length diagnostics
every time a \pkg{JAGS} model is run, and provides a more
user-friendly way to access commonly used visual convergence
diagnostics and summary statistics. Implementations of common GLMMs
are provided using a standard formula-style interface, in order to
encourage new users to explore the potential of MCMC inference without
having to generate the full code for the model themselves. A further
application of the \pkg{runjags} package is in implementing simulation
studies so that model formulations and prior specifications can be
validated using techniques such as drop-$k$ cross-validation
studies. Given that the inference made using \pkg{JAGS} and \proglang{BUGS}
can be sensitive to subtly different model specifications and prior
distributions, a user-friendly mechanism to perform these types of
analyses is potentially very useful.
\section*{Acknowledgments}
The author is grateful to the anonymous referees for their very useful
comments and suggestions, to Stefano Conti for useful discussions
regarding the Pareto family of distributions, to Vaetta Editing for
proofreading this manuscript, and to the authors of ``The \proglang{BUGS} Book''
\citep{lunn2012bugs} for kind permission to use the Salmonella
example.
\begin{thebibliography}{34}
\newcommand{\enquote}[1]{``#1''}
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{\texttt{#1}}
\providecommand{\urlprefix}{URL }
\expandafter\ifx\csname urlstyle\endcsname\relax
\providecommand{\doi}[1]{doi:\discretionary{}{}{}#1}\else
\providecommand{\doi}{doi:\discretionary{}{}{}\begingroup
\urlstyle{rm}\Url}\fi
\providecommand{\eprint}[2][]{\url{#2}}
\bibitem[{Bates \emph{et~al.}(2014)Bates, Maechler, Bolker, and
Walker}]{Bates2014}
Bates D, Maechler M, Bolker B, Walker S (2014).
\newblock \emph{\pkg{lme4}: Linear mixed-effects models using Eigen and S4}.
\newblock \proglang{R}~package version 1.1-7,
\urlprefix\url{http://CRAN.R-project.org/package=lme4}.
\bibitem[{Bolker \emph{et~al.}(2009)Bolker, Brooks, Clark, Geange, Poulsen,
Stevens, and White}]{Bolker2009}
Bolker BM, Brooks ME, Clark CJ, Geange SW, Poulsen JR, Stevens MHH, White JSS
(2009).
\newblock \enquote{{Generalized Linear Mixed Models: a Practical Guide for
Ecology and Evolution.}}
\newblock \emph{Trends in Ecology \& Evolution}, \textbf{24}(3), 127--35.
\newblock ISSN 0169-5347.
\newblock \urlprefix\url{http://dx.doi.org/10.1016/j.tree.2008.10.008}.
\bibitem[{Brooks and Roberts(1998)}]{Brooks:1998ee}
Brooks SP, Roberts GO (1998).
\newblock \enquote{{Assessing Convergence of Markov Chain Monte Carlo
Algorithms}.}
\newblock \emph{Statistics and Computing}, \textbf{8}, 319--333.
\bibitem[{Christiansen and Morris(1997)}]{Christiansen1997}
Christiansen CL, Morris CN (1997).
\newblock \enquote{{Hierarchical Poisson Regression Modeling}.}
\newblock \emph{Journal of the American Statistical Association}, \textbf{92},
618--632.
\newblock ISSN 0162-1459.
\newblock \doi{10.1080/01621459.1997.10474013}.
\newblock
\urlprefix\url{http://www.tandfonline.com/doi/abs/10.1080/01621459.1997.10474013}.
\bibitem[{Conti \emph{et~al.}(2011)Conti, Presanis, van Veen, Xiridou,
Donoghoe, {Rinder Stengaard}, and {De Angelis}}]{Conti2011}
Conti S, Presanis AM, van Veen MG, Xiridou M, Donoghoe MC, {Rinder Stengaard}
A, {De Angelis} D (2011).
\newblock \enquote{{Modeling of the HIV Infection Epidemic in the Netherlands:
a Multi-Parameter Evidence Synthesis Approach}.}
\newblock \emph{The Annals of Applied Statistics}, \textbf{5}(4), 2359--2384.
\newblock ISSN 1932-6157.
\newblock \urlprefix\url{http://dx.doi.org/10.1214/11-AOAS488}.
\bibitem[{Daniels(1999)}]{Daniels1999}
Daniels M (1999).
\newblock \enquote{{A prior for the variance in hierarchical models}.}
\newblock \emph{The Canadian Journal of Statistics}, \textbf{27}, 567--578.
\newblock ISSN 03195724.
\newblock \doi{doi: 10.2307/3316112}.
\newblock
\urlprefix\url{http://onlinelibrary.wiley.com/doi/10.2307/3316112/abstract}.
\bibitem[{Denwood(2016)}]{Denwood:2016ee}
Denwood, MJ (2016).
\newblock \enquote{{runjags: An R Package Providing Interface Utilities, Model Templates, Parallel Computing Methods and Additional Distributions for MCMC Models in JAGS}.}
\newblock \emph{Journal of Statistical Software}, \textbf{71(9)}, 1-25.
\newblock \doi{doi:10.18637/jss.v071.i09}.
\bibitem[{Denwood(2010)}]{denwood2010thesis}
Denwood MJ (2010).
\newblock \emph{{A Quantitative Approach to Improving the Analysis of Faecal
Worm Egg Count Data}}.
\newblock Doctoral thesis, University of Glasgow.
\newblock \urlprefix\url{http://theses.gla.ac.uk/1837/}.
\bibitem[{Denwood \emph{et~al.}(2010)Denwood, Reid, Love, Nielsen, Matthews,
McKendrick, and Innocent}]{Denwood:2010fk}
Denwood MJ, Reid SWJ, Love S, Nielsen MK, Matthews L, McKendrick IJ, Innocent
GT (2010).
\newblock \enquote{{Comparison of Three Alternative Methods for Analysis of
Equine Faecal egg Count Reduction Test Data}.}
\newblock \emph{Preventive Veterinary Medicine}, \textbf{93}(4), 316--23.
\newblock ISSN 1873-1716.
\newblock \urlprefix\url{http://dx.doi.org/10.1016/j.prevetmed.2009.11.009}.
\bibitem[{Denwood \emph{et~al.}(2008)Denwood, Stear, Matthews, Reid, Toft, and
Innocent}]{Denwood:2008rz}
Denwood MJ, Stear MJ, Matthews L, Reid SWJ, Toft N, Innocent GT (2008).
\newblock \enquote{{The Distribution of the Pathogenic Nematode
\textit{Nematodirus battus} in Lambs is Zero-Inflated}.}
\newblock \emph{Parasitology}, \textbf{135}(10), 1225--1235.
\newblock ISSN 1469-8161 (Electronic).
\newblock \urlprefix\url{http://dx.doi.org/10.1017/S0031182008004708}.
\bibitem[{DuMouchel(1994)}]{DuMouchel1994}
DuMouchel W (1994).
\newblock \enquote{{Hierarchical Bayes Linear Models for Meta-Analysis}.}
\newblock \emph{Technical Report~27}, {National Institute of Statistical
Sciences}.
\newblock
\urlprefix\url{http://www.niss.org/sites/default/files/pdfs/technicalreports/tr27.pdf}.
\bibitem[{Gelman(2006)}]{Gelman2006}
Gelman A (2006).
\newblock \enquote{{Prior Distributions for Variance Parameters in Hierarchical
Models}.}
\newblock \emph{Bayesian Analysis}, \textbf{1}(3), 515--533.
\bibitem[{Gelman and Rubin(1992)}]{Gelman1992}
Gelman A, Rubin DB (1992).
\newblock \enquote{{Inference from Iterative Simulation Using Multiple
Sequences}.}
\newblock \emph{Statistical Science}, \textbf{7}(4), 457--472.
\newblock \urlprefix\url{http://www.jstor.org/stable/2246093}.
\bibitem[{Gilks \emph{et~al.}(1998)Gilks, Richardson, and
Spiegelhalter}]{Gilks:1998pt}
Gilks WR, Richardson S, Spiegelhalter DJ (1998).
\newblock \emph{{Markov Chain Monte Carlo in Practice}}.
\newblock Chapman and Hall, Boca Raton, Fla.
\newblock ISBN 0412055511.
\newblock
\urlprefix\url{http://www.loc.gov/catdir/enhancements/fy0646/98033429-d.html}.
\bibitem[{Hastings(1970)}]{Hastings:1970bi}
Hastings WK (1970).
\newblock \enquote{{Monte Carlo Sampling Methods Using Markov Chains and Their
Applications}.}
\newblock \emph{Biometrika}, \textbf{57}(1), 97--109.
\newblock \urlprefix\url{http://dx.doi.org/10.1093/biomet/57.1.97}.
\bibitem[{Kass \emph{et~al.}(1998)Kass, Carlin, Gelman, and Neal}]{Kass:1998fu}
Kass RE, Carlin BP, Gelman A, Neal RM (1998).
\newblock \enquote{{Markov Chain Monte Carlo in Practice: a Roundtable
Discussion}.}
\newblock \emph{The American Statistician}, \textbf{52}(2), 93--100.
\bibitem[{Lele and Dennis(2009)}]{Lele2009}
Lele SR, Dennis B (2009).
\newblock \enquote{{Bayesian Methods for Hierarchical Models: are Ecologists
Making a Faustian Bargain?}}
\newblock \emph{Ecological Applications : a Publication of the Ecological
Society of America}, \textbf{19}(3), 581--4.
\newblock ISSN 1051-0761.
\newblock \urlprefix\url{http://www.ncbi.nlm.nih.gov/pubmed/19425420}.
\bibitem[{Lunn \emph{et~al.}(2012)Lunn, Jackson, Best, Thomas, and
Spiegelhalter}]{lunn2012bugs}
Lunn D, Jackson C, Best N, Thomas A, Spiegelhalter D (2012).
\newblock \emph{{The \proglang{BUGS} book: A practical introduction to Bayesian
analysis}}.
\newblock CRC press.
\bibitem[{Lunn \emph{et~al.}(2000)Lunn, Thomas, Best, and
Spiegelhalter}]{Lunn00}
Lunn DJ, Thomas A, Best N, Spiegelhalter D (2000).
\newblock \enquote{{\pkg{WinBUGS} - a Bayesian Modelling Framework: Concepts,
Structure, and Extensibility}.}
\newblock \emph{Statistics and Computing}, \textbf{10}(4), 325--337.
\newblock ISSN 0960-3174.
\newblock \urlprefix\url{http://dx.doi.org/10.1023/A:1008929526011}.
\bibitem[{Phillips \emph{et~al.}(2010)Phillips, Tam, Conti, Rodrigues, Brown,
Iturriza-Gomara, Gray, and Lopman}]{Phillips2010}
Phillips G, Tam CC, Conti S, Rodrigues LC, Brown D, Iturriza-Gomara M, Gray J,
Lopman B (2010).
\newblock \enquote{{Community Incidence of Norovirus-Associated Infectious
Intestinal Disease in England: Improved Estimates Using Viral Load for
Norovirus Diagnosis.}}
\newblock \emph{American Journal of Epidemiology}, \textbf{171}(9), 1014--22.
\newblock ISSN 1476-6256.
\newblock \urlprefix\url{http://dx.doi.org/10.1093/aje/kwq021}.
\bibitem[{Plummer(2003)}]{Plummer2003}
Plummer M (2003).
\newblock \enquote{{\proglang{JAGS} : A Program for Analysis of Bayesian
Graphical Models Using Gibbs Sampling}.}
\newblock In \emph{Proceedings of the 3rd International Workshop on Distributed
Statistical Computing (DSC 2003)}, pp. March 20--22,Vienna, Austria. ISSN
1609--395X.
\newblock ISSN 1609395X.
\newblock \doi{10.1.1.13.3406}.
\newblock
\urlprefix\url{http://www.ci.tuwien.ac.at/Conferences/DSC-2003/Drafts/Plummer.pdf}.
\bibitem[{Plummer(2008)}]{Plummer2008}
Plummer M (2008).
\newblock \enquote{{Penalized loss functions for Bayesian model comparison}.}
\newblock \emph{Biostatistics}, \textbf{9}, 523--539.
\newblock ISSN 14654644.
\newblock \doi{10.1093/biostatistics/kxm049}.
\bibitem[{Plummer(2016)}]{Plummer2013}
Plummer M (2016).
\newblock \emph{\pkg{rjags}: Bayesian Graphical Models using MCMC}.
\newblock \proglang{R}~package version 4-6,
\urlprefix\url{http://CRAN.R-project.org/package=rjags}.
\bibitem[{Plummer \emph{et~al.}(2006)Plummer, Best, Cowles, and
Vines}]{Plummer2006}
Plummer M, Best N, Cowles K, Vines K (2006).
\newblock \enquote{\pkg{CODA}: Convergence Diagnosis and Output Analysis for
MCMC.}
\newblock \emph{\proglang{R}~News}, \textbf{6}(1), 7--11.
\newblock \urlprefix\url{http://CRAN.R-project.org/doc/Rnews/}.
\bibitem[{Polson and Scott(2011)}]{Polson2011}
Polson NG, Scott JG (2011).
\newblock \enquote{{On the Half-Cauchy Prior for a Global Scale Parameter}.}
\newblock \emph{Cornell University Library: arXiv.org}.
\newblock \urlprefix\url{http://arxiv.org/abs/1104.4937}.
\bibitem[{{\proglang{R}~Core Team}(2015)}]{R2014}
{\proglang{R}~Core Team} (2015).
\newblock \emph{\proglang{R}: A Language and Environment for Statistical
Computing}.
\newblock \proglang{R}~Foundation for Statistical Computing, Vienna, Austria.
\newblock \urlprefix\url{http://www.R-project.org/}.
\bibitem[{Shaw \emph{et~al.}(1998)Shaw, Grenfell, and Dobson}]{Shaw:1998ty}
Shaw DJ, Grenfell BT, Dobson AP (1998).
\newblock \enquote{{Patterns of Macroparasite Aggregation in Wildlife Host
Populations.}}
\newblock \emph{Parasitology}, \textbf{117}, 597--610.
\newblock ISSN 0031-1820.
\bibitem[{Spiegelhalter \emph{et~al.}(2002)Spiegelhalter, Best, Carlin, and
van~der Linde}]{Spiegelhalter:2002mq}
Spiegelhalter DJ, Best NG, Carlin BP, van~der Linde A (2002).
\newblock \enquote{{Bayesian Measures of Model Complexity and Fit}.}
\newblock \emph{Journal of the Royal Statistical Society B}, \textbf{64}(4),
583--639.
\newblock ISSN 13697412.
\newblock \urlprefix\url{http://www.jstor.org/stable/3088806}.
\bibitem[{Tierney \emph{et~al.}(2013)Tierney, Rossini, Li, and
Sevcikova}]{Tierney2013}
Tierney L, Rossini AJ, Li N, Sevcikova H (2013).
\newblock \emph{{\pkg{snow}: Simple Network of Workstations}}.
\newblock \proglang{R}~package version 0.3-13,
\urlprefix\url{http://cran.r-project.org/package=snow}.
\bibitem[{Toft \emph{et~al.}(2007)Toft, Innocent, Gettinby, and
Reid}]{Toft:2007vl}
Toft N, Innocent GT, Gettinby G, Reid SWJ (2007).
\newblock \enquote{{Assessing the Convergence of Markov Chain Monte Carlo
Methods: an Example from Evaluation of Diagnostic Tests in Absence of a Gold
Standard.}}
\newblock \emph{Preventive Veterinary Medicine}, \textbf{79}(2-4), 244--256.
\newblock ISSN 0167-5877.
\newblock \urlprefix\url{http://dx.doi.org/10.1016/j.prevetmed.2007.01.003}.
\bibitem[{{Van Hauwermeiren} and Vose(2009)}]{D2009}
{Van Hauwermeiren} M, Vose D (2009).
\newblock \emph{{A Compendium of Distributions}}.
\newblock Vose Software, Ghent, Belgium.
\newblock \urlprefix\url{http://www.vosesoftware.com/content/ebook.pdf}.
\bibitem[{Wabersich and Vandekerckhove(2014)}]{Wabersich2014}
Wabersich D, Vandekerckhove J (2014).
\newblock \enquote{{Extending \proglang{JAGS}: a tutorial on adding custom
distributions to \proglang{JAGS} (with a diffusion model example).}}
\newblock \emph{Behavior research methods}, \textbf{46}(1), 15--28.
\newblock ISSN 1554-3528.
\newblock \doi{10.3758/s13428-013-0369-3}.
\newblock \urlprefix\url{http://www.ncbi.nlm.nih.gov/pubmed/23959766}.
\bibitem[{Wilson and Grenfell(1997)}]{Wilson199733}
Wilson K, Grenfell BT (1997).
\newblock \enquote{{Generalized Linear Modelling for Parasitologists}.}
\newblock \emph{Parasitology Today}, \textbf{13}(1), 33--38.
\newblock ISSN 0169-4758.
\newblock \urlprefix\url{http://dx.doi.org/10.1016/S0169-4758(96)40009-6}.
\bibitem[{Wilson \emph{et~al.}(1996)Wilson, Grenfell, and Shaw}]{Wilson:1996vm}
Wilson K, Grenfell BT, Shaw DJ (1996).
\newblock \enquote{{Analysis of Aggregated Parasite Distributions: a Comparison
of Methods}.}
\newblock \emph{Functional Ecology}, \textbf{10}, 592--601.
\bibitem[{Yin \emph{et~al.}(2013)Yin, Conti, Desai, Stafford, Slater, Gill, and
Simms}]{Yin2013}
Yin Z, Conti S, Desai S, Stafford M, Slater W, Gill ON, Simms I (2013).
\newblock \enquote{{The Geographic Relationship Between Sexual Health
Deprivation and the Index of Multiple Deprivation 2010: a Comparison of two
Indices.}}
\newblock \emph{Sexual Health}, \textbf{10}(2), 102--11.
\newblock ISSN 1448-5028.
\newblock \urlprefix\url{http://dx.doi.org/10.1071/SH12057}.
\end{thebibliography}
\newpage
\begin{appendix}
\section{Formulation of the negative binomial as a gamma-Poisson}
\label{sec:derivation}
The compound probability mass function of a Poisson distribution (with mean $\lambda$) integrated over a gamma distribution (with shape and scale parameters $\alpha$ and $\beta$ respectively) is given in Equation~\ref{eqn:nbgp1}.
%
\begin{align}
f(x; \alpha, \beta) &= \int_0^\infty \frac{\lambda^x}{x!} e^{-\lambda} \; . \; \beta^\alpha \frac{1}{\Gamma(\alpha)} \lambda^{\alpha-1} e^{-\beta \lambda} \; \mathrm{d} \lambda
\label{eqn:nbgp1}
\end{align}
%
Substituting $\alpha=r$ and $\beta=\frac{1-p}{p}$ into Equation~\ref{eqn:nbgp1} gives Equation~\ref{eqn:nbgp2a}, which can be re-written and simplified to Equation~\ref{eqn:nbgp2b}.
%
\begin{align}
f(x; r, p) &= \int_0^\infty \frac{\lambda^x}{x!} e^{-\lambda} \; . \; \left(\frac{1-p}{p}\right)^r \frac{1}{\Gamma(r)} \lambda^{r-1} e^{-\left(\frac{1-p}{p}\right)\lambda} \; \mathrm{d} \lambda
\label{eqn:nbgp2a} \\
&= \frac{\left(1-p\right)^r}{x! \; p^r \; \Gamma(r)} \; \int_0^\infty \lambda^{x+r-1} e^{-\lambda} e^{- \frac{\left(1-p\right)\lambda}{p}} \;\mathrm{d}\lambda \\
&= \frac{\left(1-p\right)^r}{x! \; p^r \; \Gamma(r)} \; \int_0^\infty \lambda^{x+r-1} e^{- \frac{\lambda}{p}} \;\mathrm{d}\lambda
\label{eqn:nbgp2b}
\end{align}
%
Substituting the gamma function $\frac{\Gamma(b+1)}{a^{b+1}} = \int_0^\infty t^b e^{-at} \mathrm{d}t$ for $a=\frac{1}{p}$, $b=x+r-1$ and $t=\lambda$ into Equation~\ref{eqn:nbgp2b} gives Equation~\ref{eqn:nbgp3a}.
%
\begin{align}
f(x; r, p) &= \frac{\left(1-p\right)^r}{x! \; p^r \; \Gamma(r)} \; \frac{\Gamma(x+r-1+1)}{\left(\frac{1}{p}\right)^{x+r-1+1}}
\label{eqn:nbgp3a} \\
&= \frac{\left(1-p\right)^r}{x! \; p^r \; \Gamma(r)} \; \Gamma(x+r) \, p^{x+r} \\
&= \frac{\Gamma(x+r)}{x!\;\Gamma(r)} \; p^x (1-p)^r
\label{eqn:nbgp3b}
\end{align}
%
Equation~\ref{eqn:nbgp3b} is the probability mass function of the negative binomial distribution defining the number of successes $x$ before $r$ failures with success probability $p$, which is therefore exactly equivalent to a gamma-Poisson compound distribution with mean $\frac{\alpha}{\beta}=\frac{pr}{1-p}$ and shape $\alpha=r$.
\end{appendix}
\vspace*{-1cm}
\end{document}