Lab Overview

JAGS (Just Another Gibbs Sampler) is a system that builds MCMC samplers. For a comprehensive overview of JAGS, see the user manual http://www.uvm.edu/~bbeckage/Teaching/DataAnalysis/Manuals/manual.jags.pdf. The goal today is to install JAGS and learn a little about the basic functionality by executing code from R.

1. Installing JAGS

  1. Download and install JAGS: https://sourceforge.net/projects/mcmc-jags/files/
  2. In R, install the rjags package.

2. JAGS Overview

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:

  1. Definition of the model
  2. Compilation
  3. Initialization
  4. Adaptation and burn-in
  5. Monitoring

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. We will use the coda package to assess MCMC samples.

3. Tutorial

3a Beta-Binomial Model

## JAGS lab - STAT 532
require(rjags)               # Must have previously installed package rjags.
## Loading required package: rjags
## Loading required package: coda
## Linked to JAGS 4.3.0
## Loaded modules: basemod,bugs
# Simulate the data from binomial distribution:
set.seed(09282017)
n <- 50
theta.true <- .3
y.sim <- rbinom(n, 1, theta.true)

  
dataList = list(    # Put the information into a list.
  y = y.sim ,
  Ntotal = n 
)

# Define the model as a string object:
modelString = "
model {
  for ( i in 1:Ntotal ) {
    y[i] ~ dbern( theta )
  }
  theta ~ dbeta( 1 , 1 )
}
" 
writeLines( modelString , con="TEMPmodel.txt" ) #creates a text file object

# Initialize the chains based on MLE of data.
# Option: Use function that generates random values for each chain:
initsList = function() {
  resampledY = sample( y.sim , replace=TRUE )
  thetaInit = sum(resampledY)/length(resampledY)
  thetaInit = 0.001+0.998*thetaInit # keep away from 0,1
  return( list( theta=thetaInit ) )
}

initsList = function() {
  return( list( theta=runif(1) ) )
}
initsList()
## $theta
## [1] 0.7358637

Q1. Summarize what the initsList function is doing.

# Run the chains:
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList , inits=initsList , n.chains=3 , n.adapt=500 )
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 50
##    Unobserved stochastic nodes: 1
##    Total graph size: 53
## 
## Initializing model

Q2 Describe the five arguments in the jags.model() function.

codaSamples = coda.samples( jagsModel , variable.names=c("theta") ,
                            n.iter=5000)

summary(codaSamples)
## 
## Iterations = 1:5000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 5000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean             SD       Naive SE Time-series SE 
##      0.3262505      0.0646601      0.0005279      0.0005209 
## 
## 2. Quantiles for each variable:
## 
##   2.5%    25%    50%    75%  97.5% 
## 0.2046 0.2812 0.3241 0.3694 0.4577

Q3 What is the coda.samples function doing in this case?

plot(codaSamples)