JAGS stands for Just Another Gibbs Sampler and is a program that can build a large number of Bayesian models.

In order to following along during this lecture, please download and install JAGS from https://sourceforge.net/projects/mcmc-jags/. This lecture is based on version 4.3.0 of JAGS.

The JAGS user manual, available from https://sourceforge.net/projects/mcmc-jags/files/Manuals/4.x/jags_user_manual.pdf/download, is a good to read as it provides many details on how to use JAGS.

We will initially use the rjags package, which is described in detail in the JAGS user manual. You can install and/or load it from the following syntax chunk:

if (!require(rjags)) install.packages("rjags")
## Loading required package: rjags
## Loading required package: coda
## Linked to JAGS 4.3.0
## Loaded modules: basemod,bugs
library(rjags)

Additionally, we will use the R2jags package which can be installed using the following chunk of syntax. The R2jags package does a lot of the labor intensive work for us, making it much easier (and faster) to run JAGS.

if (!require(R2jags)) install.packages("R2jags")
## Loading required package: R2jags
## 
## Attaching package: 'R2jags'
## The following object is masked from 'package:coda':
## 
##     traceplot
library(R2jags)

Finally, we will use the coda package to examine the JAGS model output for convergence and to summarize the posterior distribution. You can install and load coda from the chunk below:

if (!require(coda)) install.packages("coda")
library(coda)

To demonstrate JAGS in our introduction, we will use the syntax provided on p. 39 of Levy and Mislevy (2017) that estimates the probability parameter \(\theta\) for a 7-success in 10-trail set of data, using a binomial distribution. The text file model01.jags, which is in the current directory of the repo, has a shortened version of this syntax. Note, you can avoid having to go outside of RStudio by enclosing the model in quotes and saving it as a variable. To pass it to the jags.model() function, put the name of the variable in the textConnection function.

To run JAGS several steps must be undertaken: 1. The “data” (i.e., everything JAGS depends on and won’t sample) must be listed 2. The model must be compiled 3. The model must be “adapted” (tuned) 4. The samples must be drawn from the posterior distribution

# parts of data needing to be passed to JAGS:
J = 10
y = 7
alpha = 1
beta = 1

model01.text = "

model{
  # PRIOR DISTRIBUTION ######################################################
  theta ~ dbeta(alpha, beta)
  
  # P(X|THETA) - COND. DIST. OF DATA (AKA LIKELIHOOD) #######################
  y ~ dbin(theta, J)

}
"

# compile the JAGS model and run first part of adaptation phase (you may have to change the path to the model file)
model01 = jags.model(file = textConnection(model01.text), 
                     data = list(J = J, y = y, alpha = alpha, beta = beta), 
                     inits = list(theta=runif(1)),
                     n.chains = 5, 
                     n.adapt = 1000)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 1
##    Unobserved stochastic nodes: 1
##    Total graph size: 5
## 
## Initializing model
# not needed but included for didactic purposes (shows the algorithm that will be used)
list.samplers(model01)
## $`bugs::BinomSlicer`
## [1] "theta"
# draw samples from the posterior of the model
model01.samples = coda.samples(model = model01, 
                               variable.names = "theta", 
                               n.iter = 1000)

# samples are drawn as a coda mcmc.list object, so you can use some generic function with them:
summary(model01.samples)
## 
## Iterations = 1001:2000
## Thinning interval = 1 
## Number of chains = 5 
## Sample size per chain = 1000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean             SD       Naive SE Time-series SE 
##       0.667190       0.129928       0.001837       0.002319 
## 
## 2. Quantiles for each variable:
## 
##   2.5%    25%    50%    75%  97.5% 
## 0.3920 0.5799 0.6802 0.7627 0.8868
# plotting the chain:
plot(model01.samples)