Using JAGS for MCMC sampling

Chapter 3.3: Introduction to JAGS

In this example, we use JAGS to conduct simple linear regression. Before executing this code, but sure to install JAGS and the R package rjags.

The response is the mass of a T. Rex and the covariate is the age. The model is \[ mass_i\sim\mbox{Normal}(\beta_1+\beta_2age_i,\sigma^2).\] The priors are \(\beta_1,\beta_2\sim\mbox{Normal}(0,1000)\) and \(\sigma^2\sim\mbox{InvGamma}(0.1,0.1)\).

Load T-Rex data

 library(rjags)

 mass <- c(29.9, 1761, 1807, 2984, 3230, 5040, 5654)
 age  <- c(2, 15, 14, 16, 18, 22, 28)
 n    <- length(age)

 # JAGS require all the data to be packaged as a list
 data <- list(mass=mass,age=age,n=n) 

(1) Define the model as a string

 model_string <- textConnection("model{

   # Likelihood (dnorm uses a precision, not variance)
   for(i in 1:n){
     mass[i] ~ dnorm(beta1 + beta2*age[i],tau) #tau = 1/sigma^2
   }

   # Priors
   tau   ~  dgamma(0.1, 0.1)
   sigma <- 1/sqrt(tau)
   beta1 ~  dnorm(0, 0.001)
   beta2 ~  dnorm(0, 0.001)

 }")

(2) Load the data and compile the MCMC code

 inits <- list(beta1=rnorm(1),beta2=rnorm(1),tau=10)
 model <- jags.model(model_string,data = data, inits=inits, n.chains=2,quiet=TRUE)

(3) Burn-in for 10000 samples

 update(model, 10000, progress.bar="none")

(4) Generate 20000 post-burn-in samples and retain the parameters named in params

 params  <- c("beta1","beta2","sigma")
 samples <- coda.samples(model, 
            variable.names=params, 
            n.iter=20000, progress.bar="none")

(5) Summarize the output

 summary(samples)
## 
## Iterations = 10001:30000
## Thinning interval = 1 
## Number of chains = 2 
## Sample size per chain = 20000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean      SD Naive SE Time-series SE
## beta1    2.477   31.47   0.1573         0.1573
## beta2   51.979   39.04   0.1952         0.3849
## sigma 2809.370 1169.48   5.8474        10.0124
## 
## 2. Quantiles for each variable:
## 
##          2.5%     25%      50%     75%   97.5%
## beta1  -59.55  -18.58    2.723   23.68   63.58
## beta2  -20.88   25.25   50.475   77.65  132.84
## sigma 1095.13 2015.22 2620.528 3384.47 5647.02
 plot(samples)

plot of chunk p5