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://people.stat.sc.edu/hansont/stat740/jags_user_manual.pdf. The goal today is to install JAGS and learn a little about the basic functionality by executing code from R.

1. Installing JAGS

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.3265919      0.0646055      0.0005275      0.0005275
##
## 2. Quantiles for each variable:
##
##   2.5%    25%    50%    75%  97.5%
## 0.2082 0.2813 0.3240 0.3693 0.4594``````
Q3. What is the coda.samples function doing in this case?
``plot(codaSamples)``