• JAGS
  • background
  • model specification syntax
  • workflow
  • tips & tricks

Bayes rule for data analysis:

\[\underbrace{P(\theta \, | \, D)}_{posterior} \propto \underbrace{P(\theta)}_{prior} \times \underbrace{P(D \, | \, \theta)}_{likelihood}\]

normalizing constant:

\[ \int P(\theta') \times P(D \mid \theta') \, \text{d}\theta' = P(D) \]

easy to solve only if:

• \(\theta\) is a single discrete variable with reasonably sized domain
• \(P(\theta)\) is conjugate prior for the likelihood function \(P(D \mid \theta)\)
• we are very lucky

Markov Chain Monte Carlo

• sequence of samples from \(X \sim P\) s.t. every sample depends on its predecessor:
  1. the stationary distribution of the chain is \(P\)
  2. ideally, only few steps are necessary to get representative samples of \(P\)
• Metropolis Hastings
  • versatile but often inefficient
  • depends on proposal distribution
• Gibbs sampling
  • fast and efficient, but not universally applicable
  • depends on availability of conditional posterior

assessing quality of sample chains

• convergence / representativeness
  • trace plots
  • R hat
• efficiency
  • autocorrelation
  • effective sample size

BUGS project (1989 - present ???)

• Bayesian inference Using Gibbs Sampling
• developed by UK-based biostatisticians

WinBUGS (1997 - 2007)

• Windows based; with GUI
• programmed in component pascal

OpenBUGS (2005 - present)

• Windows & Linux; MacOS through Wine
• component pascal

JAGS (2007 - present)

• runs on Windows, Linux and MacOS
• no GUI
• written in C++

• declarative language
  • model as a directed acyclic graph
  • nodes are variables, edges are deterministic or probabilistic dependencies
• automatically selects adequate sampler
• call from & process output in R via packages:
  • R2Jags
  • rjags
  • runjags

require('rjags')
modelString = "
model{
  theta ~ dbeta(1,1)
  k ~ dbinom(theta, N)
}"
# prepare for JAGS
dataList = list(k = 14, N = 20)
# set up and run model
jagsModel = jags.model(file = textConnection(modelString), 
                       data = dataList,
                       n.chains = 2)
update(jagsModel, n.iter = 5000)
codaSamples = coda.samples(jagsModel, 
                           variable.names = c("theta"),
                           n.iter = 5000)

require('ggmcmc')
require('dplyr')
ms = ggs(codaSamples)
ms %>% group_by(Parameter) %>% 
   summarise(mean = mean(value),
             HDIlow  = HPDinterval(as.mcmc(value))[1],
             HDIhigh = HPDinterval(as.mcmc(value))[2])

## Source: local data frame [1 x 4]
## 
##   Parameter      mean    HDIlow   HDIhigh
##      (fctr)     (dbl)     (dbl)     (dbl)
## 1     theta 0.6823977 0.4998315 0.8735755

  tracePlot = ggs_traceplot(ms)  
  densPlot = ggs_density(ms) + 
    stat_function(fun = function(x) dbeta(x, 15,7), color = "black")

  require('gridExtra')
  grid.arrange(tracePlot, densPlot, ncol = 2) 

• use JAGS version 4
  • faster
  • allows for use of "="
  • added samplers and distributions (not yet all documented!?)
• syntax is a mix of BUGS and R
  • careful with parameterization of probability distributions !!!
• declarative model specification
  • directed acyclic graphs / Bayes nets
    • nodes are variables, vertices deterministic or probabilistic dependencies
  • order invariant
  • no variable reassignment, no control flow statements etc.
    • exception: "ifelse" and "step" commands for boolean switching
    • "for"-loops for vector construction
• command-line interface (documented in JAGS manual)
  • we will not use that here, but communicate via R

• you can write this into a separtate file called "myModel.jags.R"
  • this is untidy file naming but gives you default R syntax highlighting

# declare size of variables if needed
var ... ; 
# do some data massaging (usually done in R)
data{
  ...
}
# specify model
model{
  ...
}

var myData[5,100], myMu[5], mySigma[5]; 
data{
  N = sum(counts) # counts is input to JAGS from R
  indVector = c(3,2,4) # works like in R
  specialConditions = counts[indVector] # like R; new in JAGS 4
}
model{
  for (i in 1:dim(myData)[1]){
    for (j in 1:dim(myData)[2]){
      myData[i,j] ~ dnorm(myMu[i], mySigma[i])
    }
  }
  for (i in 1:5){
    tmp[i] ~ dbeta(1,1)
    myMu[i] = 100 + tmp[i] - 0.5
    mySigma[i] ~ dunif(0, 1000)
  }
}

• semicolon after declaration of variable dimensions
• declare dimensions to avoid confusing the JAGS compiler
• JAGS error messages are occassionally underinformative
• all model specification lines are probabilistic ("~") or deterministic ("=" or "<-")
  • "=" only works for JAGS 4+
• all probabilistic dependencies must be samples from a distribution known to JAGS!
  • this is therefore all ruled out:

  x ~ dbeta(1,1) - 0.5

  myCostumDistribution = ... # something fancy
  x ~ myCustomDistribution(0,1)

Your turn: what's this model good for?

• obs is a vector of 100 observations (real numbers)

model{
  mu ~ dunif(-2,2)
  var ~ dunif(0, 100)
  tau = 1/var
  for (i in 1:100){
    obs[i] ~ dnorm(mu,tau)
  }
}

NB: JAGS' precision \(\tau = 1/\sigma^2\), not standard deviation \(\sigma\) in dnorm

model{
  mu ~ dunif(-2,2)
  var ~ dunif(0, 100)
  tau = 1/var
  for (i in 1:100){
    obs[i] ~ dnorm(mu,tau)
  }
}

\[\mu \sim \text{Unif}(-2,2)\] \[\sigma^2 \sim \text{Unif}(0,100)\] \[obs_i \sim \text{Norm}(\mu, \sigma)\]

set.seed(1789)
fakeData = rnorm(200, mean = 0, sd = 1)
samplesMH = MH(f, 60000,2,10000) # this is the output of the program of Ex 4 in HW 2
                           # inferring mean and standard deviation for fakeData

modelString = "
model{
  mu ~ dunif(-2,2)
  variance ~ dunif(0,100)
  sigma = sqrt(variance)
  for (i in 1:length(obs)){
    obs[i] ~ dnorm(mu,1/variance)
  }
}"
jagsModel = jags.model(file = textConnection(modelString), 
                       data = list(obs = fakeData),
                       n.chains = 2)
update(jagsModel, n.iter = 5000)
samplesJAGS = coda.samples(jagsModel, 
                           variable.names = c("mu", "sigma"),
                           n.iter = 5000)

grid.arrange(ggs_density(ggs(samplesMH)) + theme(plot.background=element_blank()),
         ggs_density(ggs(samplesJAGS)) + theme(plot.background=element_blank()))

grid.arrange(ggs_traceplot(ggs(samplesMH)) + theme(plot.background=element_blank()),
      ggs_traceplot(ggs(samplesJAGS)) + theme(plot.background=element_blank()))

# function from Kruschke, defined in 'DBDA2E-utilities.R' 
DbdaAcfPlot(samplesMH) + theme(plot.background=element_blank())

## NULL

# function from Kruschke, defined in 'DBDA2E-utilities.R'
DbdaAcfPlot(samplesJAGS) + theme(plot.background=element_blank())

## NULL

develop step by step, & monitor each new intermediate variable

modelString = "
model{
  mu ~ dnorm(0,1)
}"
jagsModel = jags.model(file = textConnection(modelString), 
                       data = list(obs = fakeData),
                       n.chains = 2, n.adapt = 10)
update(jagsModel, n.iter = 10)
codaSamples = coda.samples(jagsModel, variable.names = c("mu"), n.iter = 5000)
ggs_density(ggs(codaSamples)) + theme(plot.background=element_blank())

model{
  mu ~ dunif(-2,2)
  variance ~ dunif(0,100)
  sigma = sqrt(variance)
  for (i in 1:length(obs)){
    obs[i] ~ dnorm(mu,1/variance) # supplied data
    obsReplicate[i] ~ dnorm(mu, 1/variance) # virtual data
  }
}

• boolean operations x <= y, x != y, x || y as usual (see manual)
• functions ifelse, equals and step for boolean conditioning

model{
  flag ~ dbern(0.5)
  parameter1 = ifelse(flag, 0, -100)
  parameter2 = ifelse(flag, 1, 100)
}

What is this model trying to achieve? And, why does it not work?

model{
  flag ~ dbern(0.5)
  SOMEPDF = ifelse(flag, dnorm, dunif)
  parameter1 = ifelse(flag, 0, -100)
  parameter2 = ifelse(flag, 1, 100)
  for(i in 1:length(obs)){
    obs[i] ~ SOMEPDF(parameter1, parameter2)
  }
}

data{
  for (i in 1:length(obs)){
    ones[i] = 1 # create a vector of ones
  }
}
model{
  flag ~ dbern(0.5)
  parameter1 = ifelse(flag, 0, -100)
  parameter2 = ifelse(flag, 1, 100)
  for(i in 1:length(obs)){
    theta[i] = ifelse(flag, 
                      dnorm(obs[i],parameter1, parameter2), 
                      dunif(obs[i],parameter1, parameter2))
    ones[i] ~ dbern(theta[i])
  }
}

• read Kruschke Chapter 9 in preparation
• optional reading:
  • Lee (2011)
  • Lee (2008)
• start on homework 3