--- a +++ b/doc/manual/jags_user_manual.tex @@ -0,0 +1,1883 @@ +\documentclass[11pt, a4paper, titlepage]{report} +\usepackage{amsmath} +\usepackage{natbib} +\usepackage{a4wide} +\usepackage{url} +\usepackage{multirow} +\usepackage{amsfonts} +\newcommand{\release}{2.2.0} +\newcommand{\JAGS}{\textsf{JAGS}} +\newcommand{\rjags}{\textsf{rjags}} +\newcommand{\BUGS}{\textsf{BUGS}} +\newcommand{\WinBUGS}{\textsf{WinBUGS}} +\newcommand{\R}{\textsf{R}} +\newcommand{\CODA}{\textsf{coda}} +\begin{document} + +\title{JAGS Version \release\ user manual} +\author{Martyn Plummer} +\maketitle + +\tableofcontents + +\chapter{Introduction} + +JAGS is Just Another Gibbs Sampler. It is a program for the analysis +of Bayesian models using Markov Chain Monte Carlo (MCMC) which is not +wholly unlike +\WinBUGS\ (\url{http://www.mrc-bsu.cam.ac.uk}). \JAGS\ was written +with three aims in mind: to have an engine for the \BUGS\ language +that runs on Unix; to be extensible, allowing users to write their own +functions, distributions, and samplers; and to be a platform for +experimentation with ideas in Bayesian modelling. + +\JAGS\ is designed to work closely with the \R\ language and +environment for statistical computation and graphics +(\url{http://www.r-project.org}). You will find it useful to install +the \CODA\ package for \R\ to analyze the output. You can also use the +\rjags\ package to work directly with \JAGS\ from within R. + +\JAGS\ is licensed under the GNU General Public License +version 2. You may freely modify and redistribute it under certain +conditions (see the file \texttt{COPYING} for details). + +\chapter{What is new in JAGS \release} + +\section{New features of JAGS} + +Modules can now be unloaded. You can inspect the list of currently +loaded modules using the \verb+LIST MODULES+ command. You can inspect +the factory objects that are currently loaded using the +\verb+LIST FACTORIES+ command and optionally turn on or off factory +objects using \verb+SET FACTORY+. + +The CODA format has been extended to allow more general monitors +to be written out to file in a rectangular format. This replaces +the former S-style dump format which is no longer available. + +\section{New modules} + +The \textsf{lecuyer} model provides random number generators with with +multiple independent streams developed by Lecyer {\em et al} +\cite{lecuyer02}. + +The \textsf{glm} module provides samplers to do block updating of the +parameters in generalized linear models +\cite{AlbertChib93,Gamerman97,HolmesHeld06,Fruhwirth-Schnatter09}. + +\section{New features in existing modules} + +\begin{itemize} +\item The \texttt{dsum} distribution in the \textsf{bugs} module is now more +flexible. It may take a variable number of arguments, and the +arguments may be either discrete-valued or real-valued (but you cannot +mix discrete-valued and real-valued arguments). +\item The \textsf{bugs} modules now includes density, distribution, and +quantile functions for many of the univariate distributions. See +table \ref{table:bugs:dpq}, section \ref{section:functions:bugs}. +\item The \textsf{mix} module contains an improved sampler for normal mixture +models defined with the \verb+dnormmix+ distribution using tempered +transtions \cite{Neal94,Celeux99} that provides much better label +switching. It also now works when the mixture probabilities are given +a Dirichlet distribution. +\end{itemize} + +\section{Old features that no longer work} + +\begin{itemize} +\item The model no longer automatically creates a \verb+deviance+ + node. You can still monitor the model deviance in the same way as + before with the command +\begin{verbatim} +monitor deviance +\end{verbatim} +and you can write the monitored values out to file using the CODA +command. However, you must now load the \verb+dic+ module first. +\item There are no longer any ``default'' monitors. When you set a monitor +you must explicitly state the node name or other summary statistic +from your model that you wish to monitor. +\item The \verb+inverse+ function only works on symmetric positive-definite +matrices. The previous version, which worked on general square matrices, +used a QR decomposition, but for medium to large matrices this would create +a slight asymmetry in the inverse, which caused errors when the inverted +matrix was used as a precision matrix in a distribution. +\end{itemize} + +\chapter{Running a model in \JAGS} + +\JAGS\ is designed for inference on Bayesian models using Markov Chain +Monte Carlo (MCMC) simulation. Running a model refers to generating +samples from the posterior distribution of the model parameters. This +takes place in five steps: +\begin{enumerate} +\item Defining the model +\item Compiling +\item Initializing +\item Adapting and burning-in +\item Monitoring +\end{enumerate} +The next stages of analysis are done outside of \JAGS: convergence +diagnostics, model criticism, and summarizing the samples must be done +using other packages more suited to this task. There are several +\R\ packages designed for analyzing MCMC output, and \JAGS\ can be +used from within \R\ using the \rjags\ package. + +\section{Definition} + +There are two parts to the definition of a model in \JAGS: a +description of the model and the definition of the data. + +\subsection{Model definition} + +The model is defined in a text file using a dialect of the +\BUGS\ language. The model definition consists of a series of +relations inside a block delimited by curly brackets \verb+{+ and +\verb+}+ and preceded by the keyword \verb+model+. Here is the standard +linear regression example: + +\begin{verbatim} +model { + for (i in 1:N) { + Y[i] ~ dnorm(mu[i], tau) + mu[i] <- alpha + beta * (x[i] - x.bar) + } + x.bar <- mean(x) + alpha ~ dnorm(0.0, 1.0E-4) + beta ~ dnorm(0.0, 1.0E-4) + sigma <- 1.0/sqrt(tau) + tau ~ dgamma(1.0E-3, 1.0E-3) +} +\end{verbatim} + +Each relation defines a node in the model in terms of other nodes that +appear on the right hand side. These are referred to as the parent +nodes. Taken together, the nodes in the model (together with the +parent/child relationships represented as directed edges) form a +directed acyclic graph. The very top-level nodes in the graph, with no +parents, are constant nodes, which are defined either in the model +definition ({\em e.g.} \verb+1.0E-3+), or in the data file ({\em + e.g.} \verb+x[1]+). + +Relations can be of two types. A {\em stochastic relation} (\verb+~+) +defines a stochastic node, representing a random variable in the +model. A {\em deterministic relation} (\verb+<-+) defines a +deterministic node, the value of which is determined exactly by the +values of its parents. + +Nodes defined by a relation are embedded in named arrays. Array names may +contain letters, numbers, decimal points and underscores, but they must +start with a letter. The node array \verb+mu+ is a vector of length +$N$ containing $N$ nodes (\verb+mu[1]+, $\ldots$, \verb+mu[N]+). The +node array \verb+alpha+ is a scalar. \JAGS\ follows the S language +convention that scalars are considered as vectors of length 1. Hence +the array \verb+alpha+ contains a single node \verb+alpha[1]+. + +Deterministic nodes do not need to be embedded in node arrays. The +node \verb+Y[i]+ could equivalently be defined as +\begin{verbatim} +Y[i] ~ dnorm(alpha + beta * (x[i] - x.bar), tau) +\end{verbatim} +In this version of the model definition, the node previously defined +as \verb+mu[i]+ still exists, but is not accessible to the user as it +does not have a name. This ability to hide deterministic nodes by +embedding them in other expressions underscores an important fact: +only the stochastic nodes in a model are really +important. Deterministic nodes are merely a syntactically convenient +way of describing the relations between, or transformations of, the +stochastic nodes. + +\subsection{Data} +\label{section:data} + +The data are defined in a separate file from the model definition, in +the format created by the \texttt{dump()} function in \R. The +simplest way to prepare your data is to read them into \R\ and then +dump them. Only numeric vectors, matrices and arrays are +allowed. More complex data structures such as factors, lists and data +frames cannot be parsed by \JAGS\, nor can non-numeric vectors. Any +\R\ attributes of the data (such as names and dimnames) are stripped +when they are read into \JAGS. + +The data may contain missing values, but you cannot supply partially +missing values for a multivariate node. In \JAGS\ a node is either +completely observed, or completely unobserved. The unobserved nodes +are referred to as the {\em parameters} of the model. The data file +therefore defines the parameters of the model by omission. + +Here are the data for the \verb+LINE+ example: +\begin{verbatim} +`x` <- +c(1, 2, 3, 4, 5) +#R-style comments, like this one, can be embedded in the data file +`Y` <- +c(1, 3, 3, 3, 5) +`N` <- +5 +\end{verbatim} + +It is an error to supply a data value for a deterministic node. (See, however, +section \ref{section:obfun} on observable functions). + +\subsection{Node Array dimensions} + +\subsubsection*{Array declarations} + +\JAGS\ allows the option of declaring the dimensions of node arrays in +the model file. The declarations consist of the keyword \texttt{var} +(for variable) followed by a comma-separated list of array names, with +their dimensions in square brackets. The dimensions may be given in +terms of any expression of the data that returns a single integer +value. + +In the linear regression example, the model block could be preceded by +\begin{verbatim} +var x[N], Y[N], mu[N], alpha, beta, tau, sigma, x.bar; +\end{verbatim} + +\subsubsection*{Undeclared nodes} + +If a node array is not declared then JAGS has three methods of +determining its size. +\begin{enumerate} +\item {\bf Using the data.} The dimension of an undeclared node array + may be inferred if it is supplied in the data file. +\item {\bf Using the left hand side of the relations.} The maximal + index values on the left hand side of a relation are taken to be the + dimensions of the node array. For example, in this case: +\begin{verbatim} +for (i in 1:N) { + for (j in 1:M) { + Y[i,j] ~ dnorm(mu[i,j], tau) + } +} +\end{verbatim} +$Y$ would be inferred to be an $N \times M$ matrix. Using this method, +empty indices are not allowed on the left hand side of any relation. +\item {\bf Using the dimensions of the parents} If a whole node array + appears on the left hand side of a relation, then its dimensions can + be inferred from the dimensions of the nodes on the right hand side. + For example, if \verb+A+ is known to be an $N \times N$ matrix + and +\begin{verbatim} +B <- inverse(A) +\end{verbatim} +Then \verb+B+ is also an $N \times N$ matrix. +\end{enumerate} + +\subsubsection*{Querying array dimensions} + +The \JAGS\ compiler has two built-in functions for querying array +sizes. The \verb+length()+ function returns the number of elements in +a node array, and the \verb+dim()+ function returns a vector +containing the dimensions of an array. These two functions may be +used to simplify the data preparation. For example, if \verb+Y+ +represents a vector of observed values, then using the \verb+length()+ +function in a for loop: +\begin{verbatim} +for (i in 1:length(Y)) { + Y[i] ~ dnorm(mu[i], tau) +} +\end{verbatim} +avoids the need to put a separate data value \verb+N+ in the file +representing the length of \verb+Y+. + +For multi-dimensional arrays, the \verb+dim+ function serves a similar +purpose. The \verb+dim+ function returns a vector, which must be stored +in an array before its elements can be accessed. For this reason, calls +to the \verb+dim+ function must always be in a data block (see section +\ref{section:data:tranformations}). +\begin{verbatim} +data { + D <- dim(Z) +} +model { + for (i in 1:D[1]) { + for (j in 1:D[2]) { + Z[i,j] <- dnorm(alpha[i] + beta[j], tau) + } + } + ... +} +\end{verbatim} +Clearly, the \verb+length()+ and \verb+dim()+ functions can only +work if the size of the node array can be inferred, using one of the +three methods outlined above. + +Note: the \verb+length()+ and \verb+dim()+ functions are different +from all other functions in \JAGS: they do not act on nodes, but only +on node {\em arrays}. As a consequence, an expression such as +\verb+dim(a %*% b)+ is syntactically incorrect. + +\section{Compilation} + +When a model is compiled, a graph representing the model is created in +computer memory. Compilation can fail for a number of reasons: +\begin{enumerate} +\item The graph contains a directed cycle. These are forbidden +in \JAGS. +\item A top-level parameter is undefined. Any node that is used on +the right hand side of a relation, but is not defined on the left +hand side of any relation, is assumed to be a constant node. Its value +must be supplied in the data file. +\item The model uses a function or distribution that has not been +defined in any of the loaded modules. +\end{enumerate} +The number of parallel chains to be run by \JAGS\ is also defined at +compilation time. Each parallel chain should produce an independent +sequence of samples from the posterior distribution. By default, +\JAGS\ only runs a single chain. + +\section{Initialization} + +Before a model can be run, it must be initialized. There are three +steps in the initialization of a model: +\begin{enumerate} +\item The initial values of the model parameters are set. +\item A Random Number Generator (RNG) is chosen for each parallel chain, + and its initial value is set. +\item The Samplers are chosen for each parameter in the model. +\end{enumerate} + +\subsection{Parameter values} + +The user may supply an initial value file containing values for the +model parameters. The file may not contain values for logical or +constant nodes. The format is the same as the data file (see section +\ref{section:data}). + +If initial values are not supplied by the user, then each parameter +chooses its own initial value based on the values of its parents. The +initial value is chosen to be a ``typical value'' from the prior +distribution. The exact meaning of ``typical value'' depends on the +distribution of the stochastic node, but is usually the mean, median, +or mode. + +If you rely on automatic initial value generation and are running +multiple parallel chains, then the initial values will be the same in +all chains. You may not want this behaviour, especially if you are +using the Gelman and Rubin convergence diagnostic, which assumes that +the initial values are over-dispersed with respect to the posterior +distribution. In this case, you are advised to set the starting values +manually using the "parameters in" statement. + +\subsection{RNGs} +\label{section:rngs} + +Each chain in \JAGS\ has its own random number generator (RNG). RNGs +are more correctly referred to as {\em pseudo}-random number +generators. They generate a sequence of numbers that merely looks +random but is, in fact, entirely determined by the initial state. You +may optionally set the name of the RNG and its initial state in the +initial values file. + +The name of the RNG is set as follows. +\begin{verbatim} +.RNG.name <- "name" +\end{verbatim} +There are four RNGs supplied by the \texttt{base} module in \JAGS\ +with the following names: +\begin{verbatim} +"base::Wichmann-Hill" +"base::Marsaglia-Multicarry" +"base::Super-Duper" +"base::Mersenne-Twister" +\end{verbatim} + +There are two ways to set the starting state of the RNG. The simplest +is to supply an integer value to \texttt{.RNG.seed}, {\em e.g.} +\begin{verbatim} +".RNG.seed" <- 314159 +\end{verbatim} +The second is way to save the state of the RNG from one JAGS session +(see the ``PARAMETERS TO'' statement, section \ref{parameters:to}) and +use this as the initial state of a new chain. The state of any RNG in +JAGS can be saved and loaded as an integer vector with the name +\texttt{.RNG.state}. For example, +\begin{verbatim} +".RNG.state" <- as.integer(c(20899,10892,29018)) +\end{verbatim} +is a valid state for the Marsaglia-Multicarry generator. You cannot +supply an arbitrary integer to \texttt{.RNG.state}. Both the length of +the vector and the permitted values of its elements are determined by +the class of the RNG. The only safe way to use \texttt{.RNG.state} is +to re-use a previously saved state. + +If no RNG names are supplied, then RNGs will be chosen automatically +so that each chain has its own independent random number stream. The +exact behaviour depends on which modules are loaded. The \texttt{base} +module uses the four generators listed above for the first four +chains, then recycles them with different seeds for the next four +chains, and so on. + +By default, \JAGS\ bases the initial state on the time stamp. This +means that, when a model is re-run, it generates an independent set of +samples. If you want your model run to be reproducible, you must +explicitly set the \verb+.RNG.seed+ for each chain. + +\subsection{Samplers} + +A Sampler is an object that acts on a set of parameters and updates +them from one iteration to the next. During initialization of the +model, Samplers are chosen automatically for all parameters. + +The Model holds an internal list of {\em Sampler Factory} objects, +which inspect the graph, recognize sets of parameters that can be +updated with specific methods, and generate Sampler objects for +them. The list of Sampler Factories is traversed in order, starting with +sampling methods that are efficient, but limited to certain specific +model structures and ending with the most generic, possibly +inefficient, methods. If no suitable Sampler can be generated for one +of the model parameters, an error message is generated. + +The user has no direct control over the process of choosing +Samplers. However, you may indirectly control the process by loading a +module that defines a new Sampler Factory. The module will insert the +new Sampler Factory at the beginning of the list, where it will be +queried before all of the other Sampler Factories. You can also +optionally turn on and off sampler factories using the ``SET FACTORY'' +command. See \ref{set:factory}. + +A report on the samplers chosen by the model, and the stochastic nodes +they act on, can be generated using the ``SAMPLERS TO'' command. See +section \ref{samplers:to}. + +\section{Adaptation and burn-in} + +In theory, output from an MCMC sampler converges to the target +distribution ({\em i.e.} the posterior distribution of the model +parameters) in the limit as the number of iterations tends to +infinity. In practice, all MCMC runs are finite. By convention, the +MCMC output is divided into two parts: an initial ``burn-in'' period, +which is discarded, and the remainder of the run, in which the output +is considered to have converged (sufficiently close) to the target +distribution. Samples from the second part are used to create +approximate summary statistics for the target distribution. + +By default, \JAGS\ keeps only the current value of each node in the +model, unless a monitor has been defined for that node. The burn-in +period of a \JAGS\ run is therefore the interval between model +initialization and the creation of the first monitor. + +When a model is initialized, it is in {\em adaptive mode}, meaning +that the Samplers used by the model may modify their behaviour for +increased efficiency. Since this adaptation may depend on the entire +sample history, the sequence generated by an adapting sampler is no +longer a Markov chain, and is not guaranteed to converge to the target +distribution. Therefore, adaptive mode must be turned off at some +point during burn-in, and a sufficient number of iterations must take +place {\em after} the adaptive phase to ensure convergence. + +By default, adaptive mode is turned off half way through first update +of a \JAGS\ model, although the user may also control the length of +the adaptive phase directly. All samplers have a built in test to +determine whether they have converged to their optimal sampling +behaviour. If any sampler fails this validation test, a warning will +be printed. To ensure optimal sampling behaviour, the model should be +run again from scratch using a longer adaptation period. + +\section{Monitoring} +\label{section:monitoring} + +A {\em monitor} in \JAGS\ is an object that records sampled +values. The simplest monitor is a {\em trace monitor}, which stores +the sampled value of a node at each iteration. + +\JAGS\ cannot monitor a node unless it has been defined in the model +file. For vector- or array-valued nodes, this means that every +element must be defined. Here is an example of a simple for loop that +only defines elements $2$ to $N$ of \verb+theta+ + +\begin{verbatim} +for (i in 2:N) { + theta[i] ~ dnorm(0,1); +} +\end{verbatim} + +Unless \verb+theta[1]+ is defined somewhere else in the model file, +the multivariate node \verb+theta+ is undefined and therefore it +will not be possible to monitor \verb+theta+ as a whole. In such +cases you can request each element separately , e.g. \verb+theta[2]+, +\verb+theta[3]+, {\em etc.}, or request a subset that is fully defined, +e.g. \verb+theta[2:6]+. + +Monitors can be classified according to whether they pool values over +iterations and whether they pool values over parallel chains (The +standard trace monitor does neither). When monitor values are written +out to file using the CODA command, the output files created depend +on the pooling of the monitor, as shown in table \ref{table:coda}. By +default, all of these files have the prefix CODA, but this may be changed +to any other name using the ``stem'' option to the CODA command +(See \ref{coda}). + +\begin{table}[h] +\begin{tabular}{ccl} +\hline +Pool & Pool & Output files \\ +iterations & chains & \\ +\hline +no & no & CODAindex.txt, CODAchain1.txt, ... + CODAchainN.txt \\ +no & yes & CODAindex0.txt, CODAchain0.txt \\ +yes & no & CODAtable1.txt, ... CODAtableN.txt \\ +yes & yes & CODAtable0.txt \\ +\hline +\end{tabular} +\caption{Output files created by the CODA command depending on whether +a monitor pools its values over chains or over iterations \label{table:coda}} +\end{table} + +The standard CODA format for monitors that do not pool values over +iterations is to create an index file and one or more output files. +The index file is has three columns with, one each line, +\begin{enumerate} +\item A string giving the name of the (scalar) value being recorded +\item The first line in the output file(s) +\item The last line in the output file(s) +\end{enumerate} +The output file(s) contain two columns: +\begin{enumerate} +\item The iteration number +\item The value at that iteration +\end{enumerate} + +Some monitors pool values over iterations. For example a mean monitor +may record only the sample mean of a node, without keeping the +individual values from each iteration. Such monitors are written out +to a table file with two columns: +\begin{enumerate} +\item A string giving the name of the (scalar) value being recorded +\item The value (pooled over all iterations) +\end{enumerate} + +\chapter{Running \JAGS} + +\JAGS\ has a command line interface. To invoke jags interactively, +simply type \texttt{jags} at the shell prompt on Unix, or the Windows +command prompt on Windows. To invoke JAGS with a script file, type +\begin{verbatim} +jags <script file> +\end{verbatim} +A typical script file has the following commands: +\begin{verbatim} +model in "line.bug" # Read model file +data in "line-data.R" # Read data in from file +compile, nchains(2) # Compile a model with two parallel chains +parameters in "line-inits.R" # Read initial values from file +initialize # Initialize the model +update 1000 # Adaptation (if necessary) and burnin for 1000 iterations +monitor alpha # Set trace monitor for node alpha ... +monitor beta # ... and beta +monitor sigma # ... and sigma +update 10000 # Update model for 10000 iterations +coda * # All monitored values are written out to file +\end{verbatim} +More examples can be found in the file \verb+classic-bugs.tar.gz+ which is +available from the \JAGS\ web page. + +Output from \JAGS\ is printed to the standard output, even when a +script file is being used. The \JAGS\ interface is designed to be +forgiving. It will print a warning message if you make a mistake, but +otherwise try to keep going. This may create a cascade of error messages, +of which only the first is informative. + +\section{Scripting commands} +\label{section:scripting} + +\JAGS\ has a simple set of scripting commands with a syntax loosely +based on \textsf{Stata}. Commands are shown below preceded by a dot +(.). This is the \JAGS\ prompt. Do not type the dot in when you are +entering the commands. + +C-style block comments taking the form /* ... */ can be +embedded anywhere in the script file. Additionally, you may use +\R-style single-line comments starting with \#. + +If a scripting command takes a file name, then the name may be +optionally enclosed in quotes. Quotes are required when the file name +contains space, or any character which is not alphanumeric, or one of +the following: \verb+_+, \verb+-+, \verb+.+, \verb+/+, \verb+\+. + +In the descriptions below, angular brackets \verb+<>+, and the text +inside them, represents a parameter that should be replaced with the +correct value by you. Anything inside square brackets \verb+[]+ is +optional. Do not type the square brackets if you wish to use an +option. + +\subsection{MODEL IN} + +\begin{verbatim} +. model in <file> +\end{verbatim} +Checks the syntactic correctness of the model description in +\texttt{file} and reads it into memory. The next compilation +statement will compile this model. + +See also: MODEL CLEAR (\ref{model:clear}) + +\subsection{DATA IN} +\label{data:in} + +\begin{verbatim} +. data in <file> +\end{verbatim} +JAGS keeps an internal data table containing the values of observed +nodes inside each node array. The DATA IN statement reads data from a +file into this data table. + +Several data statements may be used to read in data from more than one +file. If two data files contain data for the same node array, the second +set of values will overwrite the first, and a warning will be printed. + +See also: DATA TO (\ref{data:to}). + +\subsection{COMPILE} + +\begin{verbatim} +. compile [, nchains(<n>)] +\end{verbatim} +Compiles the model using the information provided in the preceding +model and data statements. By default, a single Markov chain is +created for the model, but if the \texttt{nchains} option is given, +then \texttt{n} chains are created + +Following the compilation of the model, further DATA IN statements are +legal, but have no effect. A new model statement, on the other hand, +will replace the current model. + +\subsection{PARAMETERS IN} +\label{parameters:in} + +\begin{verbatim} +. parameters in <file> [, chain(<n>)] +\end{verbatim} +Reads the values in \texttt{file} and writes them to the corresponding +parameters in chain \texttt{n}. The file has the same format as the +one in the DATA IN statement. The \texttt{chain} option may be +omitted, in which case the parameter values in all chains are set to +the same value. + +The PARAMETERS IN statement may be used before a model has been +initialized. You may only supply the values of unobserved stochastic +nodes in the parameters file. Logical nodes and constant nodes are +forbidden. + +See also: PARAMETERS TO (\ref{parameters:to}) + +\subsection{INITIALIZE} + +\begin{verbatim} +. initialize +\end{verbatim} +Initializes the model using the previously supplied data and parameter +values supplied for each chain. + +\subsection{UPDATE} + +\begin{verbatim} +. update <n> [,by(<m>)] +\end{verbatim} +Updates the model by \texttt{n} iterations. + +The first UPDATE statement turns off adaptive mode for all samplers in +the model after \texttt{n/2} iterations. A warning is printed if +adaptation is incomplete. Incomplete adaptation means that the mixing +of the Markov chain is not optimal. It is still possible to continue +with a model that has not completely adapted, but it may be preferable +to run the model again with a longer adaptation phase, starting from +the MODEL IN statement. + +A progress bar is printed on the standard output consisting of 50 +asterisks. If the \texttt{by} option is supplied, a new asterisk is +printed every \texttt{m} iterations. If this entails more than 50 +asterisks, the progress bar will be wrapped over several lines. If +\texttt{m} is zero, the printing of the progress bar is suppressed. + +\subsection{ADAPT} + +\begin{verbatim} +. adapt <n> [,by(<m>)] +\end{verbatim} +Updates the model by \texttt{n} iterations and then turns of +adaptive mode. + +Use this instead of the first UPDATE statement if you want explicit +control over the length of the adaptive sampling phase. + +\subsection{MONITOR} +\label{section:monitor} + +In \JAGS, a monitor is an object that calculates summary statistics +from a model. The most commonly used monitor simply records the value +of a single node at each iteration. This is called a ``trace'' +monitor. +\begin{verbatim} +. monitor <varname> [, thin(n)] [type(<montype>)] +\end{verbatim} +The \texttt{thin} option sets the thinning interval of the monitor so +that it will only record every nth value. The \texttt{thin} option +selects the type of monitor to create. The default type is \texttt{trace}. + +More complex monitors can be defined that do additional calculations. +For example, the \verb+dic+ module defines a ``deviance'' monitor that +records the deviance of the model at each iteration, and a ``pD'' +monitor that calculates an estimate of the effective number of +parameters on the model \cite{spiegelhalter:etal:2002}. + +\begin{verbatim} +. monitor clear <varname> [type(<montype>)] +\end{verbatim} +Clears the monitor of the given type associated with variable +\texttt{<varname>}. + +\subsection{CODA} +\label{coda} + +\begin{verbatim} +. coda <varname> [, stem(<filename>)] +\end{verbatim} +If the named node has a trace monitor, this dumps the monitored values +of to files \texttt{CODAindex.txt}, \texttt{CODAindex1.out}, +\texttt{CODAindex2.txt}, \ldots in a form that can be read by the +\CODA\ package of \R. The stem option may be used to modify the +prefix from ``CODA'' to another string. The wild-card character ``*'' +may be used to dump all monitored nodes + +\subsection{EXIT} + +\begin{verbatim} +. exit +\end{verbatim} +Exits \JAGS. \JAGS\ will also exit when it reads an end-of-file character. + +\subsection{DATA TO} +\label{data:to} +\begin{verbatim} +. data to <filename> +\end{verbatim} +Writes the data ({\em i.e.} the values of the observed nodes) to a +file in the R \texttt{dump} format. The same file can be used in a +DATA IN statement for a subsequent model. + +See also: DATA IN (\ref{data:in}) + +\subsection{PARAMETERS TO} +\label{parameters:to} +\begin{verbatim} +. parameters to <file> [, chain(<n>)] +\end{verbatim} +Writes the current parameter values ({\em i.e.} the values of the +unobserved stochastic nodes) in chain \texttt{<n>} to a file in R dump +format. The name and current state of the RNG for chain \texttt{<n>} +is also dumped to the file. The same file can be used as input in a +PARAMETERS IN statement in a subsequent run. + +See also: PARAMETERS IN (\ref{parameters:in}) + +\subsection{SAMPLERS TO} +\label{samplers:to} +\begin{verbatim} +. samplers to <file> +\end{verbatim} +Writes out a summary of the samplers to the given file. The output appears +in three tab-separated columns, with one row for each sampled node +\begin{itemize} +\item The index number of the sampler (starting with 1). The index number +gives the order in which Samplers are updated at each iteration. +\item The name of the sampler, matching the index number +\item The name of the sampled node. +\end{itemize} +If a Sampler updates multiple nodes then it is represented by multiple rows +with the same index number. + +\subsection{LOAD} +\label{load} +\begin{verbatim} +. load <module> +\end{verbatim} +Loads a module into \JAGS\ (see chapter \ref{section:modules}). Loading +a module does not affect any previously initialized models, but will +affect the future behaviour of the compiler and model initialization. + +\subsection{UNLOAD} +\label{unload} +\begin{verbatim} +. unload <module> +\end{verbatim} +Unloads a module. Currently initialized models are unaffected, but +the functions, distribution, and factory objects in the model will not +be accessible to future models. + +\subsection{LIST MODULES} +\label{list:modules} +\begin{verbatim} +. list modules +\end{verbatim} +Prints a list of the currently loaded modules. + +\subsection{LIST FACTORIES} +\label{list:factories} +\begin{verbatim} +. list factories, type(<factype>) +\end{verbatim} +List the currently loaded factory objects and whether or not they are +active. The \verb+type+ option must be given, and the three possible values +of \verb+<factype>+ are \verb+sampler+, \verb+monitor+, and \verb+rng+. + +\subsection{SET FACTORY} +\label{set:factory} +\begin{verbatim} +. set factory "<facname>" <status>, type(<factype>) +\end{verbatim} +Sets the status of a factor object. The possible values of \verb+<status>+ +are \verb+on+ and \verb+off+. Possible factory names are given from the +LIST MODULES command. + +\subsection{MODEL CLEAR} +\label{model:clear} +\begin{verbatim} +. model clear +\end{verbatim} +Clears the current model. The data table (see section \ref{data:in}) +remains intact + +\subsection{Print Working Directory (PWD)} +\begin{verbatim} +. pwd +\end{verbatim} +Prints the name of the current working directory. This is where \JAGS\ +will look for files when the file name is given without a full path, +{\em e.g.} \verb+"mymodel.bug"+. + +\subsection{Change Directory (CD)} +\begin{verbatim} +. cd <dirname> +\end{verbatim} +Changes the working directory to \texttt{<dirname>} + +\subsection{Directory list (DIR)} +\begin{verbatim} +. dir +\end{verbatim} +Lists the files in the current working directory. + +\subsection{RUN} +\begin{verbatim} +. run <cmdfile> +\end{verbatim} +Opens the file \texttt{<cmdfile>} and reads further scripting commands +until the end of the file. Note that if the file contains an EXIT +statement, then the \JAGS\ session will terminate. + + +\section{Errors} + +There are two kinds of errors in \JAGS: runtime errors, which are due to +mistakes in the model specification, and logic errors which are internal +errors in the JAGS program. + +Logic errors are generally created in the lower-level parts of the \JAGS\ +library, where it is not possible to give an informative error message. +The upper layers of the \JAGS\ program are supposed to catch such errors +before they occur, and return a useful error message that will help you +diagnose the problem. Inevitably, some errors slip through. Hence, +if you get a logic error, there is probably an error in your input to +\JAGS, although it may not be obvious what it is. Please send a bug +report (see ``Feedback'' below) whenever you get a logic error. + +Error messages may also be generated when parsing files (model files, +data files, command files). The error messages generated in this case +are created automatically by the program \texttt{bison}. They +generally take the form ``syntax error, unexpected FOO, expecting BAR'' +and are not always abundantly clear. + +If a model compiles and initializes correctly, but an error occurs +during updating, then the current state of the model will be dumped +to a file named \verb+jags.dumpN.R+ where $N$ is the chain number. +You should then load the dumped data into R to inspect the state of +each chain when the error occurred. + +\chapter{Modules} +\label{section:modules} + +The \JAGS\ library is distributed along with certain dynamically +loadable modules that extend its functionality. A module can define +new objects of the following classes: +\begin{enumerate} +\item {\bf functions} and {\bf distributions}, the basic building +blocks of the BUGS language. +\item {\bf samplers}, the objects which update the parameters of the +model at each iteration, and {\bf sampler factories}, the objects that +create new samplers for specific model structures. If the module +defines a new distribution, then it will typically also define a new +sampler for that distribution. +\item {\bf monitors}, the objects that record sampled values for +later analysis, and {\bf monitor factories} that create them. +\item {\bf random number generators}, the objects that drive the +MCMC algorithm and {\bf RNG factories} that create them. +\end{enumerate} + +The \verb+base+ module and the \verb+bugs+ module are loaded automatically +at start time. Others may be loaded by the user. + +\section{The base module} + +The base module supply the base functionality for the \JAGS\ library +to function correctly. It is loaded first by default. + +\subsection{Base Samplers} + +The \verb+base+ module defines samplers that use highly generic update +methods. These sampling methods only require basic information about +the stochastic nodes they sample. Conversely, they may not be fully +efficient. + +Three samplers are currently defined: +\begin{enumerate} +\item The Finite sampler can sample a discrete-valued node with +fixed support of less than 20 possible values. The node must not +be bounded using the \verb+T(,)+ construct +\item The Real Slice Sampler can sample any scalar real-valued +stochastic node. +\item The Discrete Slice Sampler can sample any scalar +discrete-valued stochastic node. +\end{enumerate} + +\subsection{Base RNGs} + +The \verb+base+ module defines four RNGs, taken directly from \R, +with the following names: +\begin{enumerate} +\item \verb+"base::Wichmann-Hill"+ +\item \verb+"base::Marsaglia-Multicarry"+ +\item \verb+"base::Super-Duper"+ +\item \verb+"base::Mersenne-Twister"+ +\end{enumerate} + +A single RNG factory object is also defined by the \verb+base+ +module which will supply these RNGs for chains 1 to 4 respectively, if +``RNG.name'' is not specified in the initial values file. All chains +generated by the base RNG factory are initialized using the current +time stamp. + +If you have more than four parallel chains, then the base module will +recycle the same for RNGs, but using different seeds. If you want many +parallel chains then you may wish to load the \verb+lecuyer+ module. + +\subsection{Base Monitors} + +The \verb+base+ module defines the TraceMonitor class (type +``trace''). This is the monitor class that simply records the current +value of the node at each iteration. + +\section{The bugs module} + +The \verb+bugs+ module defines some of the functions and distributions +from \WinBUGS. These are described in more detail in sections +\ref{section:functions} and \ref{section:distributions}. The +\verb+bugs+ module also defines conjugate samplers for efficient Gibbs +sampling. + +\section{The mix module} + +The \verb+mix+ module defines a novel distribution +\verb+dnormmix(mu,tau,pi)+ representing a finite mixture of normal +distributions. In the parameterization of the \verb+dnormmix+ +distribution, $\mu$, $\tau$, and $\pi$ are vectors of the same length, +and the density of \verb+y ~ dnormmix(mu, tau, pi)+ is +\[ +f(y | \mu, \tau, \pi) = \sum_i \pi_i \tau_i^{\frac{1}{2}} \phi( \tau^{\frac{1}{2}}_i (y - \mu_i)) +\] +where $\phi()$ is the probability density function of a standard +normal distribution. + +The \verb+mix+ module also defines a sampler that is designed to act +on finite normal mixtures. It uses tempered transitions to jump +between distant modes of the multi-modal posterior distribution +generated by such models \cite{Neal94,Celeux99}. The tempered +transition method is computationally very expensive. If you want to +use the \verb+dnormmix+ distribution but do not care about label +switching, then you can disable the tempered transition sampler with +\begin{verbatim} +set factory "mix::TemperedMix" off, type(sampler) +\end{verbatim} + +\section{The dic module} + +The \verb+dic+ module defines new monitor classes for Bayesian model +criticism using deviance-based measures. + +\subsection{The deviance monitor} + +The deviance monitor records the deviance of the model ({\em i.e.} the +sum of the deviances of all the observed stochastic nodes) at each +iteration. The command +\begin{verbatim} +monitor deviance +\end{verbatim} +will create a deviance monitor {\em unless} you have defined a node +called ``deviance'' in your model. In this case, you will get a trace +monitor for your deviance node. + +\subsection{The \texttt{pD} monitor} + +The \verb+pD+ monitor is used to estimate the effective number of +parameters ($p_D$) of the model \cite{spiegelhalter:etal:2002}. It +requires at least two parallel chains in the model, but calculates +a single estimate of $p_D$ across all chains \cite{plummer:2002}. +A pD monitor can be created using the command: +\begin{verbatim} +monitor pD +\end{verbatim} +Like the deviance monitor, however, if you have defined a node called +``pD'' in your model then this will take precedence, and you will get +a trace monitor for your \verb+pD+ node. + +Since the $p_D$ monitor pools its value across all chains, its values +will be written out to the index file ``CODAindex0.txt'' and +output file ``CODAoutput0.txt'' when you use the CODA command. + +The effective number of parameters is the sum of separate contributions +from all observed stochastic nodes: $p_D = \sum_i p_{D_i}$. There is +also a monitor that stores the sample mean of $p_{D_i}$. These statistics +may be used as influence diagnostics \cite{spiegelhalter:etal:2002}. +The mean monitor for $p_{D_i}$ is created with: +\begin{verbatim} +monitor pD, type(mean) +\end{verbatim} +Its values can be written out to a file ``PDtable0.txt'' with +\begin{verbatim} +coda pD, type(mean) stem(PD) +\end{verbatim} + +\subsection{The \texttt{popt} monitor} + +The \verb+popt+ monitor works exactly like the mean monitor for $p_D$, +but records contributions to the optimism of the expected deviance +($p_{opt_i}$). The total optimism $p_{opt} = \sum_i p_{opt_i}$ can be +added to the mean deviance to give the penalized expected deviance +\cite{plummer:2008}. + +The mean monitor for $p_{opt_i}$ is created with +\begin{verbatim} +monitor popt, type(mean) +\end{verbatim} +Its values can be written out to a file ``POPTtable0.txt'' with +\begin{verbatim} +coda popt, type(mean) step(POPT) +\end{verbatim} +Under asymptotically favourable conditions in which $p_{D_i} \ll 1 +\forall i$, +\[ +p_{opt} \approx 2 p_D +\] +For generalized linear models, a better approximation is +\[ +p_{opt} \approx \sum_{i=1}^n \frac {p_{D_i}}{1 - p_{D_i}} +\] + +The \verb+popt+ monitor uses importance weights to estimate +$p_{opt}$. The resulting estimates may be numerically unstable when +$p_{D_i}$ is not small. This typically occurs in random-effects +models, so it is recommended to use caution with the \verb+popt+ +until I can find a better way of estimating $p_{opt_i}$. + +\section{The msm module} + +The \verb+msm+ module defines the matrix exponential function +\verb+mexp+ and the multi-state distribution \verb+dmstate+ which +describes the transitions between observed states in continuous-time +multi-state Markov transition models. + +\section{The glm module} + +The \verb+glm+ module implements samplers for efficient updating of +generalized linear mixed models. The fundamental idea is to do block +updating of the parameters in the linear predictor. The \verb+glm+ +module is built on top of the \textsf{Csparse} sparse matrix library +\cite{Davis2006} which allows updating of both fixed and random +effects in the same block. Currently, the methods only work on +parameters that have a normal prior distribution. + +Some of the samplers are based in the idea of introducing latent +normal variables that reduce the GLM to a linear model. This idea was +introduced by Albert and Chib \cite{AlbertChib93} for probit +regression with a binary outcome, and was later refined and extended +to logistic regression with binary outcomes by Holmes and Held +\cite{HolmesHeld06}. Another approach, auxiliary mixture sampling, +was developed by Fr{\"u}hwirth-Schnatter {\em et al} +\cite{Fruhwirth-Schnatter09} and is used for more general Poisson +regression and logistic regression models with binomial outcomes. +Gamerman \cite{Gamerman97} proposed a stochastic version of the iteratively +weighted least squares algorithm for GLMs, which is also implemented +in the \verb+glm+ module. However the IWLS sampler tends to break down +when there are many random effects in the model. It uses +Metropolis-Hastings updates, and the acceptance probability may be +very small under these circumstances. + +Block updating in GLMMs frees the user from the need to center +predictor variables, like this: +\begin{verbatim} +y[i] ~ dnorm(mu[i], tau) +mu[i] <- alpha + beta * (x[i] - mean(x)) +\end{verbatim} +The second line can simply be written +\begin{verbatim} +mu[i] <- alpha + beta * x[i] +\end{verbatim} +without affecting the mixing of the Markov chain. + +\chapter{Functions} +\label{section:functions} + +Functions allow deterministic nodes to be defined using the \verb+<-+ +(``gets'') operator. Most of the functions in \JAGS\ are scalar +functions taking scalar arguments. However, \JAGS\ also allows +arbitrary vector- and array-valued functions, such as the matrix +multiplication operator \verb+%*%+ and the transpose function +\verb+t()+ defined in the \verb+bugs+ module, and the matrix +exponential function \verb+mexp()+ defined in the \verb+msm+ +module. \JAGS\ also uses an enriched dialect of the BUGS language with +a number of operators that are used in the S language. + +Scalar functions taking scalar arguments are automatically vectorized. +They can also be called when the arguments are arrays with conforming +dimensions, or scalars. So, for example, the scalar $c$ can be added to +the matrix $A$ using +\begin{verbatim} +B <- A + c +\end{verbatim} +instead of the more verbose form +\begin{verbatim} +D <- dim(A) +for (i in 1:D[1]) + for (j in 1:D[2]) { + B[i,j] <- A[i,j] + c + } +} +\end{verbatim} + +\section{Base functions} +\label{section:functions:base} + +The functions defined by the \verb+base+ module all appear as infix or +prefix operators. The syntax of these operators is built into the +\JAGS\ parser. They are therefore considered part of the modelling +language. Table \ref{table:base:functions} lists them in reverse +order of precedence. + +\begin{table}[h] +\begin{center} +\begin{tabular}{lll} +\hline +Type & Usage & Description\\ +\hline +Logical & \verb+x || y+ & Or \\ +operators & \verb+x && y+ & And \\ + & \verb+!x+ & Not \\ +\hline +Comparison & \verb+x > y+ & Greater than\\ +operators & \verb+x >= y+ & Greater than or equal to \\ + & \verb+x < y+ & Less than \\ + & \verb+x <= y+ & Less than or equal to \\ + & \verb+x == y+ & Equal \\ +\hline +Arithmetic & \verb-x + y- & Addition \\ +operators & \verb+x - y+ & Subtraction\\ + & \verb+x * y+ & Multiplication \\ + & \verb+x / y+ & Division \\ + & \verb+x %special% y+ &User-defined operators\\ + & \verb+-x+ & Unary minus\\ +\hline +Power function & \verb+x^y+ & \\ +\hline +\end{tabular} +\caption{Base functions listed in reverse order of precedence + \label{table:base:functions}} +\end{center} +\end{table} + +Logical operators convert numerical arguments to logical values: zero +arguments are converted to FALSE and non-zero arguments to +TRUE. Logical and comparison operators return the value 1 if the +result is TRUE and 0 if the result is FALSE. Comparison operators are +non-associative: the expression \verb+x < y < z+, for example, is +syntactically incorrect. + +The \verb+%special%+ function is an exception in table +\ref{table:base:functions}. It is not a function defined by the +\verb+base+ module, but is a place-holder for any function +with a name starting and ending with the character ``\verb+%+'' Such +functions are automatically recognized as infix operators by the +\JAGS\ model parser, with precedence defined by table +\ref{table:base:functions}. + +\section{Functions in the bugs module} +\label{section:functions:bugs} + +\subsection{Scalar functions} + +Table \ref{table:bugs:scalar} lists the scalar-valued functions in the +\texttt{bugs} module that also have scalar arguments. These functions +are automatically vectorized when they are given vector, matrix, or +array arguments with conforming dimensions. + +Table \ref{table:bugs:link} lists the link functions in the +\texttt{bugs} module. These are smooth scalar-valued functions that +may be specified using an S-style replacement function notation. So, +for example, the log link +\begin{verbatim} +log(y) <- x +\end{verbatim} +is equivalent to the more direct use of its inverse, the exponential +function: +\begin{verbatim} +y <- exp(x) +\end{verbatim} +This usage comes from the use of link functions in generalized linear +models. + +Table \ref{table:bugs:dpq} shows functions to calculate the +probability density, probability function, and quantiles of some of +the distributions provided by the \texttt{bugs} module. These +functions are parameterized in the same way as the corresponding +distribution. For example, if $x$ has a normal distribution with mean +$\mu$ and precision $\tau$ +\begin{verbatim} +x ~ dnorm(mu, tau) +\end{verbatim} +Then the usage of the corresponding density, probability, and quantile +functions is: +\begin{verbatim} +density.x <- dnorm(x, mu, tau) # Density of normal distribution at x +prob.x <- pnorm(x, mu, tau) # P(X <= x) +quantile90.x <- qnorm(0.9, mu, tau) # 90th percentile +\end{verbatim} +For details of the parameterization of the other distributions, see +tables \ref{table:bugs:distributions:real} and +\ref{table:bugs:distributions:discrete}. + +\begin{table} +\begin{center} +\begin{tabular}{llll} +\hline +Usage & Description & Value & Restrictions on arguments \\ +\hline +\verb+abs(x)+ & Absolute value & Real & \\ +\verb+arccos(x)+ & Arc-cosine & Real & $-1 < x < 1$\\ +\verb+arccosh(x)+ & Hyperbolic arc-cosine & Real & $1 < x$ \\ +\verb+arcsin(x)+ & Arc-sine & Real & $-1 < x < 1$\\ +\verb+arcsinh(x)+ & Hyperbolic arc-sine & Real &\\ +\verb+arctan(x)+ & Arc-tangent & Real &\\ +\verb+arctanh(x)+ & Hyperbolic arc-tangent & Real & $-1 < x < 1$\\ +\verb+cos(x)+ & Cosine & Real & \\ +\verb+cosh(x)+ & Hyperbolic Cosine & Real & \\ +\verb+cloglog(x)+ & Complementary log log & Real & $0 < x < 1$ \\ +\verb+equals(x,y)+ & Test for equality & Logical & \\ +\verb+exp(x)+ & Exponential & Real & \\ +\verb+icloglog(x)+ & Inverse complementary & Real & \\ + & log log function & \\ +\verb+ilogit(x)+ & Inverse logit & Real & \\ +\verb+log(x)+ & Log function & Real & $x > 0$ \\ +\verb+logfact(x)+ & Log factorial & Real & $x > -1$ \\ +\verb+loggam(x)+ & Log gamma & Real & $x > 0$ \\ +\verb+logit(x)+ & Logit & Real & $0 < x < 1$ \\ +\verb+phi(x)+ & Standard normal cdf & Real & \\ +\verb+pow(x,z)+ & Power function & Real & If $x < 0$ then $z$ is integer \\ +\verb+probit(x)+ & Probit & Real & $0 < x < 1$ \\ +\verb+round(x)+ & Round to integer & Integer & \\ + & away from zero & & \\ +\verb+sin(x)+ & Sine & Real & \\ +\verb+sinh(x)+ & Hyperbolic Sine & Real & \\ +\verb+sqrt(x)+ & Square-root & Real & $x >= 0$ \\ +\verb+step(x)+ & Test for $x \geq 0$ & Logical & \\ +\verb+tan(x)+ & Tangent & Real & \\ +\verb+tanh(x)+ & Hyperbolic Tangent & Real & \\ +\verb+trunc(x)+ & Round to integer & Integer & \\ + & towards zero & \\ +\hline +\end{tabular} +\caption{Scalar functions in the \texttt{bugs} module \label{table:bugs:scalar}} +\end{center} +\end{table} + +\begin{table} +\begin{center} +\begin{tabular}{llll} +\hline +Distribution & Density & Distribution & Quantile \\ +\hline +Bernoulli & dbern & pbern & qbern \\ +Beta & dbeta & pbeta & qbeta \\ +Binomial & dbin & pbin & qbin \\ +Chi-square & dchisqr & pchisqr & qchisqr \\ +Double exponential & ddexp & pdexp & qdexp \\ +Exponential & dexp & pexp & qexp \\ +F & df & pf & qf \\ +Gamma & dgamma & pgamma & qgamma \\ +Generalized gamma & dgengamma & pgengamma & qgengamma \\ +Hypergeometric & dhyper & phyper & qhyper \\ +Log-normal & dlnorm & plnorm & qlnorm \\ +Negative binomial & dnegbin & pnegbin & qnegbin \\ +Normal & dnorm & pnorm & qnorm \\ +Pareto & dpar & ppar & qpar \\ +Poisson & dpois & ppois & qpois \\ +Student t & dt & pt & qt \\ +Weibull & dweib & pweib & qweib \\ +\hline +\end{tabular} +\caption{Wrappers for the d-p-q functions from the \texttt{Rmath} library + \label{table:bugs:dpq}} +\end{center} +\end{table} + +\begin{table} +\begin{center} +\begin{tabular}{llll} +\hline +Link function & Description & Range & Inverse \\ +\hline +\verb+cloglog(y) <- x+ & Complementary log log & $0 < y < 1$ & \verb+y <- icloglog(x)+ \\ +\verb+log(y) <- x+ & Log & $0 < y$ & \verb+y <- exp(x)+ \\ +\verb+logit(y) <- x+ & Logit & $0 < y < 1$ & \verb+y <- ilogit(x)+ \\ +\verb+probit(y) <- x+ & Probit & $0 < y < 1$ & \verb+y <- phi(x)+\\ +\hline +\end{tabular} +\caption{Link functions in the \texttt{bugs} module \label{table:bugs:link}} +\end{center} +\end{table} + +\section{Scalar-valued functions with vector arguments} + +Table \ref{table:bugs:scalar2} lists the scalar-valued functions in the +\texttt{bugs} module that take general arguments. Unless otherwise +stated in table \ref{table:bugs:scalar2}, the arguments to these functions +may be scalar, vector, or higher-dimensional arrays. + +The \verb+max()+ and \verb+min()+ functions work like the +corresponding \R\ functions. They take a variable number of arguments +and return the maximum/minimum element over all supplied +arguments. This usage is compatible with \WinBUGS, although more general. + +\begin{table} +\begin{tabular}{lll} +\hline +Function & Description & Restrictions \\ +\hline +\verb+inprod(x1,x2)+ & Inner product & Dimensions of $a$, $b$ conform \\ +\verb+interp.lin(e,v1,v2)+ & Linear Interpolation & $e$ scalar, \\ + & & $v1,v2$ conforming vectors \\ +\verb+logdet(a)+ & Log determinant & $a$ is a square matrix \\ +\verb+max(x1,x2,...)+ & Maximum element among all arguments & \\ +\verb+mean(x)+ & Mean of elements of $a$ & \\ +\verb+min(x1,x2,...)+ & Minimum element among all arguments & \\ +\verb+prod(x)+ & Product of elements of $a$ & \\ +\verb+sum(a)+ & Sum of elements of $a$& \\ +\verb+sd(a)+ & Standard deviation of elements of $a$ & \\ +\hline +\end{tabular} +\caption{Scalar-valued functions with general + arguments in the \texttt{bugs} module \label{table:bugs:scalar2}} +\end{table} + +\section{Vector- and array-valued functions} + +Table \ref{table:bugs:vector} lists vector- or matrix-valued functions +in the \texttt{bugs} module. + +The \texttt{sort} and \texttt{rank} functions behaves like their R +namesakes: \texttt{sort} accepts a vector and returns the same values +sorted in ascending order; \texttt{rank} returns a vector of ranks. +This is distinct from \WinBUGS, which has two scalar-valued functions +\verb+rank+ and \verb+ranked+. + +\begin{table} +\begin{center} +\begin{tabular}{lll} +\hline +Usage & Description & Restrictions \\ +\hline +\verb+inverse(a)+ & Matrix inverse & $a$ is a symmetric positive definite matrix \\ +\verb+mexp(a)+ & Matrix exponential & $a$ is a square matrix \\ +\verb+rank(v)+ & Ranks of elements of $v$ & $v$ is a vector \\ +\verb+sort(v)+ & Elements of $v$ in order & $v$ is a vector \\ +\verb+t(a)+ & Transpose & $a$ is a matrix \\ +\verb+a %*% b+ & Matrix multiplication & $a,b$ conforming vector or matrices\\ + +\hline +\end{tabular} +\caption{Vector- or matrix-valued functions in the \texttt{bugs} + module \label{table:bugs:vector}} +\end{center} +\end{table} + +\chapter{Distributions} +\label{section:distributions} + +Distributions are used to define stochastic nodes using the \verb+~+ +operator. The distributions defined in the bugs module are listed in +table \ref{table:bugs:distributions:real} (real-valued distributions), +\ref{table:bugs:distributions:discrete} (discrete-valued +distributions), and \ref{table:bugs:distributions:multi} +(multivariate distributions). + +Some distributions have restrictions on the valid parameter values, +and these are indicated in the tables. If a Distribution is +given invalid parameter values when evaluating the log-likelihood, it +returns $-\infty$. When a model is initialized, all stochastic nodes +are checked to ensure that the initial parameter values are valid for +their distribution. + + +\begin{table} + \begin{center} + \begin{tabular}{llcll} + \hline + Name & Usage & Density & Lower & Upper \\ + \hline + Beta & \verb+dbeta(a,b)+ & + \multirow{2}{*}{ + $\frac{\textstyle x^{a-1}(1-x)^{b-1}}{\textstyle \beta(a,b)}$ + } & $0$ & $1$ \\ + & $a > 0$, $b > 0$ \\ + Chi-square & \verb+dchisqr(k)+ & + \multirow{2}{*}{ + $\frac{\textstyle x^{\frac{k}{2} - 1} \exp(-x/2)} + {\textstyle 2^{\frac{k}{2}} \Gamma({\scriptstyle \frac{k}{2}})}$ + } & 0 & \\ + & $k > 0$ \\ + Double & \verb+ddexp(mu,tau)+ & + \multirow{2}{*}{$\tau \exp(-\tau | x-\mu |)/2$} & & \\ + exponential & $\tau > 0$ \\ + Exponential & \verb+dexp(lambda)+ & + \multirow{2}{*}{$\lambda \exp(-\lambda x)$} & 0 & \\ + & $\lambda > 0$ \\ + F & \verb+df(n,m)+ & + \multirow{2}{*}{ + $\textstyle \frac{\Gamma(\frac{n + m}{2})} + {\Gamma(\frac{n}{2}) \Gamma(\frac{m}{2})} + \left(\frac{n}{m} \right)^{\frac{n}{2}} x^{\frac{n}{2} - 1} + \left\{1 + \frac{nx}{m} \right\}^{-\frac{(n + m)}{2}}$} & 0 & \\ + & $n > 0$, $m > 0$ \\ + Gamma & \verb+dgamma(r, mu)+ & + \multirow{2}{*}{ + $\frac{\textstyle \mu^r x^{r - 1} \exp(-\mu x)} + {\textstyle \Gamma(r)}$} & 0 & \\ + & $\mu > 0$, $r > 0$ \\ + Generalized & \verb+dgen.gamma(r,mu,beta)+ & + \multirow{2}{*}{ + $\beta \mu^{\beta r} x^{\beta r - 1} \exp\{-(\mu x)^{\beta}\}$ + } & $0$ & \\ + gamma & $\mu >0$, $\beta > 0$, $r > 0$ \\ + Log-normal & \verb+dlnorm(mu,tau)+ & + \multirow{2}{*}{ + $\tau^{\frac{1}{2}} x^{-1} \exp \left\{-\tau (\log(x) - \mu)^2/2 \right\}$} & 0 \\ + ~ & $\tau > 0$ \\ + Normal & \verb+dnorm(mu,tau)+ & + \multirow{2}{*}{ + $\left(\frac{\tau}{2\pi}\right)^{\frac{1}{2}} \exp\{-(x - \mu)^2 \tau\}$} & & \\ + ~ & $\tau > 0$ \\ + Pareto & \verb+dpar(alpha, c)+ & + \multirow{2}{*}{ + $\alpha c^{\alpha} x^{-(\alpha + 1)}$ + } & $c$ & \\ + ~ & $\alpha > 0$, $c > 0$ \\ + Student t & \verb+dt(mu,tau,k)+ & + \multirow{2}{*}{ + $\textstyle \frac{\Gamma(\frac{k+1}{2})}{\Gamma(\frac{k}{2})} + \left(\frac{\tau}{k\pi} \right)^{\frac{1}{2}} + \left\{1 + \frac{\tau (x - \mu)^2}{k} \right\}^{-\frac{(k+1)}{2}}$} & & \\ + ~ & $\tau > 0$, $k > 0$ \\ + Uniform & \verb+dunif(a,b)+ & + \multirow{2}{*}{$\frac{\textstyle 1}{\textstyle b - a}$} & $a$ & $b$ \\ + ~ & $a < b$ \\ + Weibull & \verb+dweib(v, lambda)+ & + \multirow{2}{*}{$v \lambda x^{v - 1} \exp (- \lambda x^v)$} & 0 & \\ + ~ & $v > 0$, $\lambda > 0$ \\ + \hline + \end{tabular} + \caption{Univariate real-valued distributions in the \texttt{bugs} module + \label{table:bugs:distributions:real}} + \end{center} +\end{table} + +\begin{table} + \begin{center} + \begin{tabular}{llllll} + \hline + Name & Usage & Density & Lower & Upper \\ + \hline + Bernoulli & \verb+dbern(p)+ & + \multirow{2}{*}{$p^x (1 - p)^{1 -x}$} & + $0$ & $1$ \\ + ~ & $0 < p < 1$ \\ + Binomial & \verb+dbin(p,n)+ & + %\multirow{2}{*}{$\begin{array}{c} n \\ x \end{array} p^x (1-p)^{n-x}$} + \multirow{2}{*}{${n \choose x} p^x (1-p)^{n-x}$} + ~ & $0$ & $n$ \\ + ~ & $0 < p < 1$, $n \in \mathbb{N}^+$ \\ + Categorical & \verb+dcat(p)+ & \multirow{2}{*}{$\frac{\textstyle p_x}{\textstyle \sum_i p_i}$} & $1$ & $N$ \\ + ~ & $p \in (\mathbb{R}^+)^N$ \\ + Hypergeometric & \verb+dhyper(n1,n2,m1,psi)+ & + \multirow{2}{*}{ + $\textstyle {n_1 \choose x} {n_2 \choose m_1 - x} \psi^x$ + } & + $\text{max}(0,n_+ - m_1)$ & $\text{min}(n_1,m_1)$ \\ + ~ & $0 \leq n_i$, $0 < m_1 \leq n_+$ \\ + Negative & \verb+dnegbin(p, r)+ & + \multirow{2}{*}{${x + r -1 \choose x} p^r (1-p)^x$} & 0 & \\ + binomial & $0 < p < 1$, $r \in \mathbb{N}^+$ \\ + Poisson & \verb+dpois(lambda)+ & + \multirow{2}{*}{$\frac{\textstyle \exp(-\lambda) \lambda^x}{\textstyle x!}$} & 0 & \\ + ~ & $\lambda > 0$ \\ + \hline + \end{tabular} + \caption{Discrete univariate distributions in the \texttt{bugs} module + \label{table:bugs:distributions:discrete}} + \end{center} +\end{table} + + +\begin{table} + \begin{center} + \begin{tabular}{lll} + \hline + Name & Usage & Density \\ + \hline + Dirichlet & \verb+p ~ ddirch(alpha)+ & + \multirow{2}{*}{$\Gamma(\sum_i \alpha_i) \prod_j + \frac{\textstyle p_j^{\alpha_j - 1}}{\textstyle \Gamma(\alpha_j)}$} \\ + ~ & $\alpha_j \geq 0$ \\ + & \\ + Multivariate & \verb+x ~ dmnorm(mu,Omega)+ & + \multirow{2}{*}{ + $\left(\frac{|\Omega|}{2\pi}\right)^{\frac{1}{2}} exp\{-(x-\mu)^T \Omega (x-\mu) / 2\}$} \\ + normal & $\Omega$ positive definite \\ + Wishart & \verb+Omega ~ dwish(R,k)+ & + \multirow{2}{*}{ + $\frac{\textstyle |\Omega|^{(k-p-1)/2} |R|^{k/2} \exp\{-\text{Tr}(R\Omega/2)\}} + {\textstyle 2^{pk/2} \Gamma_p (k/2)}$ + } \\ + & $R$ pos. def. \\ + Multivariate & \verb+x ~ dmt(mu,Omega,k)+ & + \multirow{2}{*}{ + $\frac{\textstyle \Gamma \{(k+p)/2\}}{\textstyle \Gamma(k/2) (n\pi)^{p/2}} + |\Omega|^{1/2} + \left\{1 + \frac{1}{k} (x - \mu)^T \Omega (x - \mu) \right\}^{-\frac{(k+p)}{2}}$ } \\ + Student t & $\Omega$ pos. def. & \\ + Multinomial & \verb+x ~ dmulti(p, n)+ & + \multirow{2}{*}{$n! \prod_j + \frac{\textstyle p_j^{x_j}}{\textstyle x_j!}$} \\ + ~ & $\sum_i x_i = n$ \\ + & \\ + \hline + \end{tabular} + \caption{Multivariate distributions in the \texttt{bugs} module + \label{table:bugs:distributions:multi}} + \end{center} +\end{table} + +\chapter{Differences between \JAGS\ and \WinBUGS} + +Although \JAGS\ aims for the same functionality as \WinBUGS, there are +a number of important differences. + +\subsection{Data format} + +There is no need to transpose matrices and arrays when transferring +data between \R\ and \JAGS, since \JAGS\ stores the values of an array +in ``column major'' order, like \R\ and FORTRAN ({\em i.e.} filling +the left-hand index first). + +If you have an \textsf{S}-style data file for \WinBUGS\ and you wish +to convert it for \JAGS, then use the command \texttt{bugs2jags}, +which is supplied with the \CODA\ package. + +\subsection{Distributions} + +Structural zeros are allowed in the Dirichlet distribution. If +\begin{verbatim} +p ~ ddirch(alpha) +\end{verbatim} +and some of the elements of alpha are zero, then the corresponding +elements of p will be fixed to zero. + +The Multinomial (\verb+dmulti+) and Categorical (\verb+dcat+) +distributions, which take a vector of probabilities as a parameter, +may use unnormalized probabilities. The probability vector is +normalized internally so that +\[ +p_i \rightarrow \frac{p_i}{\sum_j p_j} +\] + +\subsection{Observable Functions} +\label{section:obfun} + +Logical nodes in the \BUGS\ language are a convenient way of +describing the relationships between observables (constant and +stochastic nodes), but are not themselves observable. You cannot +supply data values for a logical node. + +This restriction can occasionally be inconvenient, as there are +important cases where the data are a deterministic function of +unobserved variables. Two important examples are +\begin{enumerate} +\item Censored data, which commonly occurs in survival analysis. In +the most general case, we know that unobserved failure time $T$ +lies in the interval $(L,U]$. +\item Aggregate data when we observe the sum of two or more +unobserved variables. +\end{enumerate} +\JAGS\ contains two novel distributions to handle these situations. +\begin{enumerate} +\item The \texttt{dinterval} distribution represents interval-censored +data. It has two parameters: $t$ the original continuous variable, and +$c[]$, a vector of cut points of length $M$, say. If \texttt{X $\sim$ +dinterval(t, c)} then + +\begin{tabular}{lll} +$X=0$ & if & $t \leq c[1]$\\ +$X=m$ & if & $c[m] < t \leq c[m+1]$ for $1 \leq m < M$\\ +$X = M$ & if & $c[M] < t$. +\end{tabular} + +\item The \texttt{dsum} distribution represents the sum of two or more +variables. It takes a variable number of parameters. If \texttt{Y $\sim$ +dsum(x1,x2,x3)} then $Y=x1+x2+x3$. +\end{enumerate} +These distributions exist to give a likelihood to data that is, in fact, +a deterministic function of the parameters. The relation +\begin{verbatim} +Y ~ dsum(x1, x2) +\end{verbatim} +is logically equivalent to +\begin{verbatim} +Y <- x1 + x2 +\end{verbatim} +But the latter form does not create a contribution to the likelihood, +and does not allow you to define $Y$ as data. The likelihood function +is trivial: it is 1 if the parameters are consistent with the data and +0 otherwise. The \texttt{dsum} distribution also requires a special +sampler, which can currently only handle the case where the parameters +of \texttt{dsum} are unobserved stochastic nodes, and where the +parameters are either all discrete-valued or all continuous-valued. A node +cannot be subject to more than one \texttt{dsum} constraint. + +\subsection{Data transformations} +\label{section:data:tranformations} + +\JAGS\ allows data transformations, but the syntax is different from +\BUGS. \BUGS\ allows you to put a stochastic node twice on the left +hand side of a relation, as in this example taken from the manual +\begin{verbatim} + for (i in 1:N) { + z[i] <- sqrt(y[i]) + z[i] ~ dnorm(mu, tau) + } +\end{verbatim} +This is forbidden in \JAGS. You must put data transformations in a +separate block of relations preceded by the keyword \texttt{data}: +\begin{verbatim} +data { + for (i in 1:N) { + z[i] <- sqrt(y[i]) + } +} +model { + for (i in 1:N) { + z[i] ~ dnorm(mu, tau) + } + ... +} +\end{verbatim} +This syntax preserves the declarative nature of the \BUGS\ language. +In effect, the data block defines a distinct model, which describes +how the data is generated. Each node in this model is forward-sampled +once, and then the node values are read back into the data table. The +data block is not limited to logical relations, but may also include +stochastic relations. You may therefore use it in simulations, +generating data from a stochastic model that is different from the one +used to analyse the data in the \texttt{model} statement. + +This example shows a simple location-scale problem in which the ``true'' +values of the parameters \texttt{mu} and \texttt{tau} are generated +from a given prior in the \texttt{data} block, and the generated +data is analyzed in the \texttt{model} block. +\begin{verbatim} +data { + for (i in 1:N) { + y[i] ~ dnorm(mu.true, tau.true) + } + mu.true ~ dnorm(0,1); + tau.true ~ dgamma(1,3); +} +model { + for (i in 1:N) { + y[i] ~ dnorm(mu, tau) + } + mu ~ dnorm(0, 1.0E-3) + tau ~ dgamma(1.0E-3, 1.0E-3) +} +\end{verbatim} +Beware, however, that every node in the \texttt{data} statement will +be considered as data in the subsequent \texttt{model} statement. This +example, although superficially similar, has a quite different interpretation. +\begin{verbatim} +data { + for (i in 1:N) { + y[i] ~ dnorm(mu, tau) + } + mu ~ dnorm(0,1); + tau ~ dgamma(1,3); +} +model { + for (i in 1:N) { + y[i] ~ dnorm(mu, tau) + } + mu ~ dnorm(0, 1.0E-3) + tau ~ dgamma(1.0E-3, 1.0E-3) +} +\end{verbatim} +Since the names \texttt{mu} and \texttt{tau} are used in both +\texttt{data} and \texttt{model} blocks, these nodes will be +considered as {\em observed} in the model and their values will be +fixed at those values generated in the \texttt{data} block. + +\subsection{Directed cycles} + +Directed cycles are forbidden in \JAGS. There are two important +instances where directed cycles are used in \BUGS. +\begin{itemize} +\item Defining autoregressive priors +\item Defining ordered priors +\end{itemize} +For the first case, the \texttt{GeoBUGS} extension to \WinBUGS\ provides +some convenient ways of defining autoregressive priors. These should be +available in a future version of \JAGS. + +\subsection{Censoring, truncation and prior ordering} +\label{section:censoring} + +These are three, closely related issues that are all handled using +the \texttt{I(,)} construct in \BUGS. + +Censoring occurs when a variable $X$ is not observed directly, +but is observed only to lie in the range $(L,U]$. Censoring is +an {\em a posteriori} restriction of the data, and is represented +in WinBUGS by the \texttt{I(,)} construct, {\em e.g.} +\begin{verbatim} +X ~ dnorm(theta, tau) I(L,U) +\end{verbatim} +where $L$ and $U$ are constant nodes. + +Truncation occurs when a variable is known {\em a priori} to lie in +a certain range. Although \BUGS\ has no construct for representing +truncated variables, it turns out that there is no difference between +censoring and truncation for top-level parameters ({\em i.e.} variables +with no unobserved parents). Hence, for example, this +\begin{verbatim} +theta ~ dnorm(0, 1.0E-3) I(0, ) +\end{verbatim} +is a perfectly valid way to describe a parameter $\theta$ with a +half-normal prior distribution. + +Prior ordering occurs when a vector of nodes is known {\em a priori} +to be strictly increasing or decreasing. It can be represented in +WinBUGS with symmetric $I(,)$ constructs, {\em e.g.} +\begin{verbatim} +X[1] ~ dnorm(0, 1.0E-3) I(,X[2]) +X[2] ~ dnorm(0, 1.0E-3) I(X[1],) +\end{verbatim} +ensures that $X[1] \leq X[2]$. + +\JAGS\ makes an attempt to separate these three concepts. + +Censoring is handled in \JAGS\ using the new distribution +\texttt{dinterval} (section \ref{section:obfun}). This can be +illustrated with a survival analysis example. A right-censored +survival time $t_i$ with a Weibull distribution is described in +\WinBUGS\ as follows: +\begin{verbatim} +t[i] ~ dweib(r, mu[i]) I(c[i], ) +\end{verbatim} +where $t_i$ is unobserved if $t_i > c_i$. In \JAGS\ this becomes +\begin{verbatim} +is.censored[i] ~ dinterval(t[i], c[i]) +t[i] ~ dweib(r, mu[i]) +\end{verbatim} +where \verb+is.censored[i]+ is an indicator variable that takes the +value 1 if $t_i$ is censored and 0 otherwise. See the MICE and KIDNEY +examples in the ``classic bugs'' set of examples. + +Truncation is represented in \JAGS\ using the \texttt{T(,)} construct, +which has the same syntax as the \texttt{I(,)} construct in \WinBUGS, +but has a different interpretation. If +\begin{verbatim} +X ~ dfoo(theta) T(L,U) +\end{verbatim} +then {\em a priori} $X$ is known to lie between $L$ and $U$. This +generates a likelihood +\[ +\frac{p(x \mid \theta)}{P(L \leq X \leq U \mid \theta)} +\] +if $L \leq X \leq U$ and zero otherwise, where $p(x \mid \theta)$ is +the density of $X$ given $\theta$ according to the distribution +\texttt{foo}. Note that calculation of the denominator may be +computationally expensive. + +Prior ordering of top-level parameters in the model can be achieved +using the \texttt{sort} function, which sorts a vector in ascending +order. + +Symmetric truncation relations like this +\begin{verbatim} +alpha[1] ~ dnorm(0, 1.0E-3) I(,alpha[2]) +alpha[2] ~ dnorm(0, 1.0E-3) I(alpha[1],alpha[3]) +alpha[3] ~ dnorm(0, 1.0E-3) I(alpha[2],) +\end{verbatim} +Should be replaced by this +\begin{verbatim} +for (i in 1:3) { + alpha0[i] ~ dnorm(0, 1.0E-3) +} +alpha[1:3] <- sort(alpha0) +\end{verbatim} + +\chapter{Feedback} + +Please send feedback to \url{martyn_plummer@users.sourceforge.net}. +I am particularly interested in the following problems: + +\begin{itemize} +\item Crashes, including both segmentation faults and uncaught exceptions. +\item Incomprehensible error messages +\item Models that should compile, but don't +\item Output that cannot be validated against \WinBUGS +\item Documentation erors +\end{itemize} + +If you want to send a bug report, it must be reproducible. Send the +model file, the data file, the initial value file and a script file +that will reproduce the problem. Describe what you think should +happen, and what did happen. + +\chapter{Acknowledgments} + +Many thanks to the \BUGS\ development team, without whom \JAGS\ would +not exist. Thanks also to Simon Frost for pioneering \JAGS\ on +Windows and Bill Northcott for getting \JAGS\ on Mac OS X to +work. Kostas Oikonomou found many bugs while getting \JAGS\ to work on +Solaris using Sun development tools and libraries. Bettina Gruen, +Chris Jackson, Greg Ridgeway and Geoff Evans also provided useful +feedback. Special thanks to Jean-Baptiste Denis who has been very +diligent in providing feedback on JAGS. + +Testing of \JAGS\ on IRIX 6.5 was carried out on Helix Systems at the +National Institutes of Health, Bethesda, MD (\url{http://helix.nih.gov}). + +\bibliographystyle{plain} +\bibliography{jags_user_manual} + +\end{document} + + +