Hierarchical Modelling - Practical 3

Hierarchical modelling using JAGS

We will continue with our hemlock light and growth example. Recall from the slides that we are modelling the relationship between plant growth rate and light as a non-linear curve defined by the Michaelis-Menten function - our process model. The equation here also indexes [s] for site which we come to later in the practical:

\[a[s] * (x[i,s]-c) / ((a[s]/b)+(x[i,s]-c))\]

# hemlock <- read.csv("Hemlock-light-data-simple.csv")
hemlock <- read.csv(url("https://raw.githubusercontent.com/NERC-CEH/beem_data/main/Hemlock-light-data-simple.csv"))
str(hemlock)

'data.frame':   77 obs. of  2 variables:
 $ Light               : num  98.5 98.5 98.5 23.5 11.8 11.8 19.4 19.4 98.7 98.7 ...
 $ Observed.growth.rate: num  18.2 26 25.7 6.2 10 11.2 7.3 15.7 44.2 39.7 ...

head(hemlock)

  Light Observed.growth.rate
1  98.5                 18.2
2  98.5                 26.0
3  98.5                 25.7
4  23.5                  6.2
5  11.8                 10.0
6  11.8                 11.2

Let’s read in the data, and take a quick look at how it’s structured and also plot it. We take a rudimentary approach to modelling the curve first, using non-linear least squares regression (a non-linear version of the familiar OLS).

par(mfrow=c(1,1))
plot(hemlock$Light, hemlock$Observed.growth.rate)

# Quick and dirty non-linear least squares fit to visualise
attach(hemlock)
start<-list(a=40, b=2, c=6)
m1<-nls(Observed.growth.rate~a*(Light-c)/((a/b)+(Light-c)), start=start)
summary(m1)

Formula: Observed.growth.rate ~ a * (Light - c)/((a/b) + (Light - c))

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
a   38.493      4.364   8.820 3.64e-13 ***
b    1.735      0.498   3.485 0.000831 ***
c    4.740      2.067   2.293 0.024664 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.409 on 74 degrees of freedom

Number of iterations to convergence: 6 
Achieved convergence tolerance: 5.505e-06

a<-coef(m1)["a"]
b<-coef(m1)["b"]
c<-coef(m1)["c"]
plot(Light,Observed.growth.rate)
x<-Light
curve(a*(x-c)/((a/b)+(x-c)), add=TRUE, col="blue")

detach(hemlock)

Now lets move ahead with the Bayesian analysis for a single non-linear fit using the software JAGS and the R package rjags, which ‘speaks to’ JAGS via our model file. First we need to set up the conditions for the model - how we want to run the MCMC chains, the data etc. Then we call to JAGS using our model file (open it and take a closer look).

Something worth noting is that we have tended to refer to our variance (sigma) in our model specifications, but JAGS uses precision, the inverse of variance (which typically takes the notation tau).

jags_model <- "
# JAGS model

# y is the observed growth rate
# x is the measurement of light, L

# a is alpha, the max. growth rate at infinite light
# c is light level at which growth is zero
# b is rate at which curve tails off (gamma)

########################################################################################
model{
  ### likelihood
  # Data model
  for (i in 1:n)
  {
    y[i] ~ dnorm(mu[i], tau.o) #Note JAGS uses precision, not variance
  }  
  
  # process model
  for (i in 1:n)
  {
    mu[i] ~ dnorm(mu2[i], tau.p) #Note JAGS uses precision, not variance
    
    mu2[i] <- a * (x[i]-c) / ((a/b)+(x[i]-c))
    
  }
  
  # priors
  a ~ dgamma(0.01,0.01)
  c ~ dunif(-10,10)
  b ~ dgamma(0.01,0.01)
  
  sigma.o ~ dnorm(5, 1/(0.5*0.5)) ## assume prior knowledge of observation error (5 with SD of 0.5)
  sigma.p ~ dunif(0, 50)
  
  # derived quantities: convert precisions to standard deviations
  tau.o <- 1/(sigma.o * sigma.o)
  tau.p <- 1/(sigma.p * sigma.p)
  
} # end of model
#########################################################################################
"
writeLines(jags_model, con="jags_model.txt")

library(rjags)

Loading required package: coda

Linked to JAGS 4.3.1

Loaded modules: basemod,bugs

# Initialize the MCMC chains. This needs to be a list, so we define a function that calls the list. Anything that is to be ESTIMATED in the model needs to be initialised.
inits <- function() list(a=runif(1,35,40), b=runif(1,1,3), c=runif(1,-5,5), 
                         sigma.o=1, sigma.p=1) #small tau makes variance big
inits

function() list(a=runif(1,35,40), b=runif(1,1,3), c=runif(1,-5,5), 
                         sigma.o=1, sigma.p=1)

# Specify data to be used by your JAGS program (giving dataframe and variables); must be a list. 
data=list(
    n=nrow(hemlock),
    x=as.numeric(hemlock$Light),
    y=as.numeric(hemlock$Observed.growth.rate )
    )
data

$n
[1] 77

$x
 [1] 98.5 98.5 98.5 23.5 11.8 11.8 19.4 19.4 98.7 98.7 98.7 10.9  6.5  6.5 11.1
[16] 50.4 57.5 56.3 53.6 53.5 57.0 67.8 70.9 75.1 22.6 64.3 12.8 29.0 39.6 29.6
[31] 28.7 71.1 98.0 98.0 98.0  5.6  6.9 39.1 28.8 38.2 43.4 38.2 56.4 53.5 29.4
[46] 39.1 18.7 18.0 25.4 57.1 25.9 15.5 14.6 42.9 37.1 66.4 39.3 14.9 14.9 27.4
[61] 19.7 22.9 98.3 98.3 98.3  7.7  7.7  8.8  8.8  8.8 41.3 43.4 40.5 80.5 82.8
[76] 83.5 85.9

$y
 [1] 18.2 26.0 25.7  6.2 10.0 11.2  7.3 15.7 44.2 39.7 30.0  4.8  4.2  5.3  4.7
[16] 28.2 22.7 36.3 37.3 35.3 37.3 33.8 32.0 34.2 12.5 26.2 11.5  7.5 25.0 14.3
[31] 17.2 28.3 12.7 30.3 25.0  2.7  3.5 19.2 30.7 32.2 46.7 24.8 31.7 26.8 29.7
[46] 30.8 28.7 20.2  4.7 25.5  6.8 10.5  4.2 26.3 31.3 13.5 14.3 16.8  8.2 11.7
[61] 23.2 21.2 37.0 34.2 29.2  5.3  6.0  6.5  4.0  6.7 16.8 17.0 23.2 31.7 31.2
[76] 23.5 27.3

# Parameters monitored (i.e., for which estimates are saved)
# params <- c("a","b","c","sigma.o","sigma.p")

# call to JAGS
# Set run conditions: number of iterations for adaptation & runs, number of chains, etc.
n.adapt=500 #plays with algorithm for first 500 samples
n.update = 1000 #burnin
n.iter = 5000 #no. of iterations
jm = jags.model(file = "jags_model.txt",data=data,inits,n.chains=3,n.adapt=n.adapt)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 77
   Unobserved stochastic nodes: 82
   Total graph size: 484

Initializing model

# the jags file name needs to be in the same working directory, or set the path name
# Burnin the chain
update(jm, n.iter=n.update)

The next block of code creates the chains and stores them as an MCMC list. We give the model object name created from when we used the jags.model function, and then variable.names tells JAGS which variables to “monitor” (i.e. the variables for which we want posterior distributions). Then we specify how many iterations we want and to keep - n.thin = 10 means we keep every 10th iteration (which we may choose to do to reduce autocorrelation in our chain).

From summary(objectname) we can get summary statistics for the MCMC chains. The full output of the coda object is a list of matrices, where each matrix is the output of the chains for each parameter (the columns). We can also use the object for plotting (see further below).

# generate coda object for exploring parameters and deviance.
zm <- coda.samples(jm,variable.names=c("a", "b", "c","sigma.o","sigma.p"),n.iter=n.iter, n.thin=10)
summary(zm)

Iterations = 1501:6500
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
a       40.033 4.8197 0.039353        0.21033
b        1.657 0.4970 0.004058        0.02417
c        3.488 2.5880 0.021131        0.09995
sigma.o  5.050 0.5188 0.004236        0.01375
sigma.p  5.434 1.0624 0.008675        0.03912

2. Quantiles for each variable:

           2.5%    25%    50%    75%  97.5%
a       32.4016 36.687 39.508 42.672 51.504
b        0.9156  1.329  1.589  1.906  2.736
c       -2.9772  2.150  3.922  5.364  7.155
sigma.o  4.0340  4.701  5.051  5.400  6.075
sigma.p  3.2996  4.796  5.474  6.137  7.386

We also create a jags object, which in brief, are more handy for posterior distributions in matrix form, e.g. posterior distributions of predictions (mu).

# generate JAGS object
zj <- jags.samples(jm,variable.names=c("a", "b","c", "mu2"), n.iter=n.iter, n.thin=10)
zj

$a
mcarray:
[1] 40.42775

Marginalizing over: iteration(5000),chain(3) 

$b
mcarray:
[1] 1.636964

Marginalizing over: iteration(5000),chain(3) 

$c
mcarray:
[1] 3.438509

Marginalizing over: iteration(5000),chain(3) 

$mu2
mcarray:
 [1] 31.252374 31.252374 31.252374 17.407380  9.595709  9.595709 15.206086
 [8] 15.206086 31.266739 31.266739 31.266739  8.738852  3.611863  3.611863
[15]  8.933883 25.605556 26.853602 26.658128 26.196286 26.178564 26.772855
[22] 28.322637 28.702968 29.180810 16.960254 27.861541 10.488908 19.788620
[29] 23.169854 20.017178 19.672290 28.726665 31.216260 31.216260 31.216260
[36]  2.310076  4.157379 23.036601 19.711221 22.791232 24.116874 22.791232
[43] 26.674639 26.178564 19.941591 23.036601 14.783639 14.345482 18.293875
[50] 26.789084 18.515174 12.636980 11.959759 23.998504 22.481290 28.142435
[57] 23.090160 12.189518 12.189518 19.151645 15.382548 17.111358 31.237963
[64] 31.237963 31.237963  5.193592  5.193592  6.511420  6.511420  6.511420
[71] 23.607366 24.116874 23.404403 29.739596 29.960728 30.026186 30.244399

Marginalizing over: iteration(5000),chain(3) 

hat <- summary(zj$mu,quantile,c(0.025,.5,.975))$stat # gives you the median, upper and lower quaniles for mu
hat

          [,1]     [,2]     [,3]     [,4]      [,5]      [,6]     [,7]     [,8]
2.5%  28.23156 28.23156 28.23156 14.87633  6.818390  6.818390 12.68836 12.68836
50%   31.21294 31.21294 31.21294 17.42503  9.573266  9.573266 15.19961 15.19961
97.5% 34.58827 34.58827 34.58827 19.89686 12.406633 12.406633 17.78859 17.78859
          [,9]    [,10]    [,11]     [,12]     [,13]     [,14]     [,15]
2.5%  28.24160 28.24160 28.24160  5.812118 -1.469351 -1.469351  6.047795
50%   31.22733 31.22733 31.22733  8.719984  3.644365  3.644365  8.915143
97.5% 34.61263 34.61263 34.61263 11.705875  8.336314  8.336314 11.873749
         [,16]    [,17]    [,18]    [,19]    [,20]    [,21]    [,22]    [,23]
2.5%  23.52044 24.74626 24.58066 24.13528 24.11346 24.67574 26.05322 26.35741
50%   25.61757 26.84556 26.65139 26.19765 26.17978 26.76605 28.32132 28.70399
97.5% 27.64412 28.90911 28.70141 28.23284 28.21355 28.82034 30.58803 31.03822
         [,24]    [,25]    [,26]     [,27]    [,28]    [,29]    [,30]    [,31]
2.5%  26.74622 14.42191 25.64784  7.827993 17.28696 20.88262 17.52849 17.17100
50%   29.18640 16.97236 27.85897 10.469083 19.81339 23.20662 20.04318 19.69856
97.5% 31.62641 19.47348 30.05707 13.183900 22.15291 25.29021 22.36855 22.04168
         [,32]    [,33]    [,34]    [,35]     [,36]      [,37]    [,38]
2.5%  26.37978 28.20400 28.20400 28.20400 -3.752676 -0.5975828 20.73794
50%   28.72810 31.17766 31.17766 31.17766  2.401274  4.1681636 23.07240
97.5% 31.06641 34.52930 34.52930 34.52930  7.655669  8.6369946 25.16510
         [,39]    [,40]    [,41]    [,42]    [,43]    [,44]    [,45]    [,46]
2.5%  17.20871 20.47008 21.92290 20.47008 24.59556 24.11346 17.44629 20.73794
50%   19.73717 22.82680 24.14928 22.82680 26.66799 26.17978 19.96784 23.07240
97.5% 22.07720 24.93781 26.19549 24.93781 28.71958 28.21355 22.29860 25.16510
         [,47]    [,48]    [,49]    [,50]    [,51]    [,52]     [,53]    [,54]
2.5%  12.26051 11.82537 15.76761 24.68812 15.99589 10.08411  9.375828 21.79058
50%   14.77598 14.33734 18.31328 26.78169 18.53402 12.61378 11.936981 24.02861
97.5% 17.37232 16.95348 20.74220 28.83910 20.95606 15.24440 14.570216 26.07298
         [,55]    [,56]    [,57]     [,58]     [,59]    [,60]    [,61]    [,62]
2.5%  20.13035 25.88891 20.79389  9.604614  9.604614 16.63635 12.86041 14.57136
50%   22.51292 28.14144 23.12544 12.164334 12.164334 19.17319 15.37792 17.12541
97.5% 24.64966 30.37341 25.21313 14.802457 14.802457 21.55790 17.95711 19.61894
         [,63]    [,64]    [,65]    [,66]    [,67]     [,68]     [,69]
2.5%  28.22086 28.22086 28.22086 1.041210 1.041210  2.974819  2.974819
50%   31.19793 31.19793 31.19793 5.192152 5.192152  6.490925  6.490925
97.5% 34.56454 34.56454 34.56454 9.257478 9.257478 10.100777 10.100777
          [,70]    [,71]    [,72]    [,73]    [,74]    [,75]    [,76]    [,77]
2.5%   2.974819 21.35670 21.92290 21.13591 27.17002 27.33254 27.38925 27.54159
50%    6.490925 23.64043 24.14928 23.43748 29.73819 29.95441 30.02138 30.24058
97.5% 10.100777 25.70859 26.19549 25.51186 32.35998 32.66454 32.74783 33.06291

# install.packages("lattice")
library(lattice) # this library gives us some additional plotting options
# Plot parameters
# You can set ask to TRUE
plot(zm, ask = FALSE)

plot(zm[, c("a", "b")], ask=FALSE) # subset of parameters

xyplot(zm,ask = FALSE)

densityplot(zm, ask = FALSE)

From these objects we can now plot our predictions along with the data.

# Plot predictions
bu <- summary(zj$mu2, quantile, c(.025,.5,.976))$stat
pl <- cbind(hemlock$Light,t(bu))
dat <- pl[order(pl[,1]),]

plot(hemlock$Light, hemlock$Observed.growth.rate, xlab="Light", ylab="Growth Rate", pch=16, col="blue")
lines(dat[,1], dat[,2],lty="dashed")
lines(dat[,1], dat[,3])
lines(dat[,1], dat[,4],lty="dashed")

We have not covered diagnostics here, but the code below gives some examples of what is possible using the ‘coda’ object you created.

# Convergence diagnostics
rejectionRate(zm) # sampling conjugate

      a       b       c sigma.o sigma.p 
      0       0       0       0       0 

gelman.diag(zm) # var in chains, stable=1

Potential scale reduction factors:

        Point est. Upper C.I.
a             1.01       1.02
b             1.00       1.01
c             1.01       1.03
sigma.o       1.01       1.02
sigma.p       1.01       1.01

Multivariate psrf

1.01

heidel.diag(zm) # requires convergence

[[1]]
                                      
        Stationarity start     p-value
        test         iteration        
a       passed       1         0.9440 
b       passed       1         0.6242 
c       passed       1         0.7882 
sigma.o passed       1         0.0686 
sigma.p passed       1         0.0600 
                                 
        Halfwidth Mean  Halfwidth
        test                     
a       passed    40.26 0.6720   
b       passed     1.62 0.0616   
c       passed     3.35 0.3265   
sigma.o passed     5.11 0.0523   
sigma.p passed     5.30 0.1541   

[[2]]
                                      
        Stationarity start     p-value
        test         iteration        
a       passed       1         0.675  
b       passed       1         0.301  
c       passed       1         0.362  
sigma.o passed       1         0.496  
sigma.p passed       1         0.511  
                                 
        Halfwidth Mean  Halfwidth
        test                     
a       passed    39.91 0.6422   
b       passed     1.66 0.0628   
c       passed     3.65 0.2998   
sigma.o passed     5.00 0.0407   
sigma.p passed     5.53 0.1048   

[[3]]
                                      
        Stationarity start     p-value
        test         iteration        
a       passed         1       0.0574 
b       passed       501       0.4722 
c       passed         1       0.2157 
sigma.o passed         1       0.1021 
sigma.p passed         1       0.1034 
                                 
        Halfwidth Mean  Halfwidth
        test                     
a       passed    39.93 0.8159   
b       passed     1.66 0.0774   
c       failed     3.47 0.3860   
sigma.o passed     5.04 0.0462   
sigma.p passed     5.47 0.1349   

raftery.diag(zm) # how many iter you need for convergence

[[1]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95 
                                               
         Burn-in  Total Lower bound  Dependence
         (M)      (N)   (Nmin)       factor (I)
 a       18       20676 3746          5.52     
 b       33       32151 3746          8.58     
 c       22       24880 3746          6.64     
 sigma.o 9        9308  3746          2.48     
 sigma.p 67       66809 3746         17.80     


[[2]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95 
                                               
         Burn-in  Total Lower bound  Dependence
         (M)      (N)   (Nmin)       factor (I)
 a       13       14172 3746         3.78      
 b       32       35108 3746         9.37      
 c       27       29031 3746         7.75      
 sigma.o 12       12728 3746         3.40      
 sigma.p 18       18442 3746         4.92      


[[3]]

Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95 
                                               
         Burn-in  Total Lower bound  Dependence
         (M)      (N)   (Nmin)       factor (I)
 a       20       21639 3746          5.78     
 b       60       56485 3746         15.10     
 c       30       34364 3746          9.17     
 sigma.o 11       12009 3746          3.21     
 sigma.p 50       54338 3746         14.50     

load.module("dic")

module dic loaded

dic.ex <- dic.samples(jm,n.iter, type="pD")

# autocorrelation with chains
autocorr.plot(zm)

crosscorr(zm)

                  a            b            c      sigma.o     sigma.p
a        1.00000000 -0.773464160 -0.543704600 -0.018331895  0.02256673
b       -0.77346416  1.000000000  0.631290122 -0.006723272  0.04116641
c       -0.54370460  0.631290122  1.000000000 -0.005288184 -0.02045251
sigma.o -0.01833190 -0.006723272 -0.005288184  1.000000000 -0.53399238
sigma.p  0.02256673  0.041166407 -0.020452513 -0.533992376  1.00000000

crosscorr.plot(zm)

Multi-site extension to the model

Now we will explore adding another level to the model - the nature of the relationship between light and plant growth rate can vary according to site, of which we have three in our new dataset. We will model that the maximum growth rate (the asymptote or parameter a) can vary according to site such that this parameter is a random variable with a mean and variance. Lets explore our data briefly again.

jags_model <- "
# JAGS model
########################################################################################
model{
  ### likelihood
  # Data model
    for(s in 1:3){
      for(i in 1:n){
        y[i,s] ~ dnorm(mu[i,s], tau.o) #Note JAGS uses precision, not variance
      }  
    }
  # process model
    for(s in 1:3){
      for(i in 1:n){
        mu[i,s] ~ dnorm(mu2[i,s], tau.p) #Note JAGS uses precision, not variance
    
        mu2[i,s] <- a[s] * (x[i,s]-c) / ((a[s]/b)+(x[i,s]-c))
        }
    }
  
  # priors
  for(s in 1:3){
    a[s] ~ dnorm(mu.a, tau.a)
  }
  mu.a ~ dunif(0,100)
  c ~ dunif(-10,10)
  b ~ dgamma(0.01,0.01)
  
  sigma.o ~ dnorm(5, 1/(0.5*0.5)) ## assume prior knowledge of observation error (5 with SD of 0.5)
  sigma.p ~ dunif(0, 50)
  sigma.a ~ dunif(0, 50)
  
  # derived quantities: convert precisions to standard deviations
  tau.o <- 1/(sigma.o * sigma.o)
  tau.p <- 1/(sigma.p * sigma.p)
  tau.a <- 1/(sigma.a * sigma.a)
  
} # end of model
#########################################################################################
"
writeLines(jags_model, con="jags_model.txt")

# hemlock2 <- read.csv("Hemlock-light-data-hierarchical_2.csv")
# hemlock <- read.csv(url("https://raw.githubusercontent.com/NERC-CEH/beem_data/main/Hemlock-light-data-simple.csv"))
hemlock2 <- read.csv(url("https://raw.githubusercontent.com/NERC-CEH/beem_data/main/Hemlock-light-data-hierarchical_2.csv"))
str(hemlock2) # 231 observations

'data.frame':   231 obs. of  3 variables:
 $ Light               : num  5.6 6.5 6.5 6.9 7.7 7.7 8.8 8.8 8.8 10.9 ...
 $ Observed.growth.rate: num  2.7 4.2 5.3 3.5 5.3 6 6.5 4 6.7 4.8 ...
 $ Site                : int  1 1 1 1 1 1 1 1 1 1 ...

head(hemlock2)

  Light Observed.growth.rate Site
1   5.6                  2.7    1
2   6.5                  4.2    1
3   6.5                  5.3    1
4   6.9                  3.5    1
5   7.7                  5.3    1
6   7.7                  6.0    1

par(mfrow=c(1,1))
colour <- c("red", "blue", "green")
col.list <- rep(0, length(hemlock2$Site))
col.list[hemlock2$Site=="1"] <- 1 
col.list[hemlock2$Site=="2"] <- 2
col.list[hemlock2$Site=="3"] <- 3
plot(hemlock2$Light, hemlock2$Observed.growth.rate, col=c(colour[col.list]))

Again, we need to set up the conditions for the model - how we want to run the MCMC chains, the data etc. Notice how we have to adapt our initial starting values for the new parameteres in the model - open up the new multi-level file.

Note: when you have two product symbols in your conditional distribution with different indices, you have nested for loops.

# Initialize the MCMC chains. This needs to be a list.
inits2 <- function() list(a=rep(30,3), b=runif(1,1,3), c=runif(1,-5,5), 
                         mu.a=30, sigma.o=1, sigma.p=1, sigma.a=1) #small tau makes variance big
inits2

function() list(a=rep(30,3), b=runif(1,1,3), c=runif(1,-5,5), 
                         mu.a=30, sigma.o=1, sigma.p=1, sigma.a=1)

## create dataframe for y data (observed growth rate)

y.dat <- array(NA, c(77,3))
y.dat[1:77,1] <- hemlock2$Observed.growth.rate[which(hemlock2$Site==1)]
y.dat[1:77,2] <- hemlock2$Observed.growth.rate[which(hemlock2$Site==2)]
y.dat[1:77,3] <- hemlock2$Observed.growth.rate[which(hemlock2$Site==3)]
head(y.dat)

     [,1] [,2] [,3]
[1,]  2.7  4.7  2.6
[2,]  4.2  2.7  5.5
[3,]  5.3  5.9  5.2
[4,]  3.5  2.3  4.3
[5,]  5.3  4.7  4.2
[6,]  6.0  7.0  6.6

## create dataframe for x data (Light)
x.dat <- array(NA, c(77,3))
x.dat[1:77,1] <- hemlock2$Light[which(hemlock2$Site==1)]
x.dat[1:77,2] <- hemlock2$Light[which(hemlock2$Site==2)]
x.dat[1:77,3] <- hemlock2$Light[which(hemlock2$Site==3)]
head(x.dat)

     [,1] [,2] [,3]
[1,]  5.6  3.9  5.5
[2,]  6.5  4.8  5.7
[3,]  6.5  6.7  5.9
[4,]  6.9  7.1  6.5
[5,]  7.7  7.6  7.8
[6,]  7.7  8.3  7.8

# Specify data; must be a list. 
data=list(
  n=77,
  x=x.dat,
  y=y.dat
)
data

$n
[1] 77

$x
      [,1] [,2]  [,3]
 [1,]  5.6  3.9   5.5
 [2,]  6.5  4.8   5.7
 [3,]  6.5  6.7   5.9
 [4,]  6.9  7.1   6.5
 [5,]  7.7  7.6   7.8
 [6,]  7.7  8.3   7.8
 [7,]  8.8  8.7   8.0
 [8,]  8.8  9.2   8.3
 [9,]  8.8  9.3   9.2
[10,] 10.9  9.3  10.7
[11,] 11.1  9.4  11.1
[12,] 11.8 10.4  11.6
[13,] 11.8 11.5  12.2
[14,] 12.8 11.5  13.5
[15,] 14.6 13.1  13.8
[16,] 14.9 14.0  15.1
[17,] 14.9 14.9  15.6
[18,] 15.5 16.0  17.4
[19,] 18.0 17.5  17.5
[20,] 18.7 18.1  17.9
[21,] 19.4 18.4  18.6
[22,] 19.4 20.0  18.9
[23,] 19.7 21.0  20.6
[24,] 22.6 21.1  21.0
[25,] 22.9 21.5  22.4
[26,] 23.5 22.9  23.9
[27,] 25.4 24.0  24.3
[28,] 25.9 24.3  25.1
[29,] 27.4 27.1  28.3
[30,] 28.7 27.3  28.6
[31,] 28.8 28.9  28.9
[32,] 29.0 29.0  29.1
[33,] 29.4 30.2  30.5
[34,] 29.6 31.2  30.7
[35,] 37.1 36.2  36.4
[36,] 38.2 38.2  37.4
[37,] 38.2 38.4  37.8
[38,] 39.1 38.8  38.2
[39,] 39.1 39.7  39.8
[40,] 39.3 40.1  39.8
[41,] 39.6 40.5  40.1
[42,] 40.5 40.7  40.8
[43,] 41.3 41.0  41.4
[44,] 42.9 41.8  41.8
[45,] 43.4 42.8  42.7
[46,] 43.4 44.2  43.7
[47,] 50.4 49.4  51.5
[48,] 53.5 52.0  51.8
[49,] 53.5 53.6  52.9
[50,] 53.6 54.6  55.1
[51,] 56.3 55.9  55.5
[52,] 56.4 57.3  56.8
[53,] 57.0 57.9  57.2
[54,] 57.1 58.0  57.5
[55,] 57.5 58.7  58.3
[56,] 64.3 64.9  62.5
[57,] 66.4 66.2  65.1
[58,] 67.8 66.6  68.5
[59,] 70.9 69.1  69.7
[60,] 71.1 72.2  70.2
[61,] 75.1 76.1  74.3
[62,] 80.5 80.8  81.2
[63,] 82.8 82.0  82.2
[64,] 83.5 84.1  84.7
[65,] 85.9 85.0  86.6
[66,] 98.0 96.8  96.2
[67,] 98.0 97.0  96.8
[68,] 98.0 97.1  97.2
[69,] 98.3 97.5  97.3
[70,] 98.3 98.2  97.6
[71,] 98.3 98.5  98.2
[72,] 98.5 99.4  98.5
[73,] 98.5 99.4  98.8
[74,] 98.5 99.5  98.9
[75,] 98.7 99.6  99.7
[76,] 98.7 99.7 100.0
[77,] 98.7 99.9 100.3

$y
      [,1]      [,2]     [,3]
 [1,]  2.7  4.700000  2.60000
 [2,]  4.2  2.700000  5.50000
 [3,]  5.3  5.900000  5.20000
 [4,]  3.5  2.300000  4.30000
 [5,]  5.3  4.700000  4.20000
 [6,]  6.0  7.000000  6.60000
 [7,]  6.5  7.000000  4.30000
 [8,]  4.0  2.700000  3.00000
 [9,]  6.7  5.000000  6.50000
[10,]  4.8  5.800000  5.30000
[11,]  4.7  3.200000  3.00000
[12,] 10.0  8.200000 10.60000
[13,] 11.2  9.300000 12.30000
[14,] 11.5  9.700000  9.90000
[15,]  4.2  3.200000  1.40000
[16,] 16.8 15.200000 17.60000
[17,]  8.2  6.700000  7.70000
[18,] 10.5 11.500000  8.60000
[19,] 20.2 19.900000 20.20000
[20,] 28.7 27.700000 28.50000
[21,]  7.3  7.800000  6.50000
[22,] 15.7 15.300000 16.90000
[23,] 23.2 22.800000 22.80000
[24,] 12.5 13.000000 14.80000
[25,] 21.2 22.600000 22.20000
[26,]  6.2  7.700000  9.60000
[27,]  4.7  3.900000  5.50000
[28,]  6.8  5.400000 10.60000
[29,] 11.7 12.300000  8.70000
[30,] 17.2 15.500000 18.10000
[31,] 30.7 32.500000 33.70000
[32,]  7.5  6.700000  5.00000
[33,] 29.7 31.300000 29.70000
[34,] 14.3 15.700000 12.50000
[35,] 31.3 47.136898 63.89346
[36,] 32.2 28.957552 53.82420
[37,] 24.8 33.877043 85.72795
[38,] 19.2 32.241437 69.38528
[39,] 30.8 23.643437 35.94980
[40,] 14.3 39.105755 40.44215
[41,] 25.0 43.813701 37.19800
[42,] 23.2 32.523262 46.49454
[43,] 16.8 21.814450 49.02929
[44,] 26.3 30.259437 50.90613
[45,] 46.7 18.943678 66.51076
[46,] 17.0 45.630525 52.30740
[47,] 28.2 20.099826 90.76589
[48,] 35.3 18.316243 87.56340
[49,] 26.8 15.446236 39.93826
[50,] 37.3  8.010457 50.31908
[51,] 36.3 31.162236 67.13457
[52,] 31.7 31.934466 59.57232
[53,] 37.3 16.407239 50.69410
[54,] 25.5 19.764969 55.26655
[55,] 22.7 41.048627 56.47245
[56,] 26.2 34.287661 67.75295
[57,] 13.5 31.308202 76.56862
[58,] 33.8 16.114231 40.12719
[59,] 32.0 32.443813 48.78845
[60,] 28.3 32.601453 41.05356
[61,] 34.2 37.722595 50.32160
[62,] 31.7 28.820051 57.95963
[63,] 31.2 17.373145 56.57387
[64,] 23.5 31.298986 54.23894
[65,] 27.3 28.693491 67.87955
[66,] 12.7 39.714315 26.83006
[67,] 30.3 39.255526 39.53883
[68,] 25.0 36.226181 52.91146
[69,] 37.0 15.376946 60.52105
[70,] 34.2 44.137594 56.51756
[71,] 29.2 29.533540 71.09812
[72,] 18.2 37.031727 71.52053
[73,] 26.0 21.899740 63.90977
[74,] 25.7 21.608192 69.75787
[75,] 44.2 24.773272 43.06696
[76,] 39.7 45.428707 31.40532
[77,] 30.0 21.766299 57.56964

# Parameters monitored (i.e., for which estimates are saved)
# params <- c("a","b","c","sigma.o","sigma.p", "sigma.a", "mu.a")

Then we call to JAGS using our model file and the rjags function jags.model.

### fit the model in JAGS
n.adapt = 500
n.update = 1000
n.iter = 5000
jm2 = jags.model(file = "jags_model.txt", data = data, inits2, n.chains = 3, n.adapt=n.adapt)

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 231
   Unobserved stochastic nodes: 240
   Total graph size: 1528

Initializing model

# Burn in the chain
update(jm2, n.iter=n.update)

# Generate coda object
zm2 <- coda.samples(jm2,variable.names=c("a","b","c","sigma.o","sigma.p","sigma.a", "mu.a"), n.iter=n.iter,thin=1)
summary(zm)

Iterations = 1501:6500
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
a       40.033 4.8197 0.039353        0.21033
b        1.657 0.4970 0.004058        0.02417
c        3.488 2.5880 0.021131        0.09995
sigma.o  5.050 0.5188 0.004236        0.01375
sigma.p  5.434 1.0624 0.008675        0.03912

2. Quantiles for each variable:

           2.5%    25%    50%    75%  97.5%
a       32.4016 36.687 39.508 42.672 51.504
b        0.9156  1.329  1.589  1.906  2.736
c       -2.9772  2.150  3.922  5.364  7.155
sigma.o  4.0340  4.701  5.051  5.400  6.075
sigma.p  3.2996  4.796  5.474  6.137  7.386

# Generate jags object
zj2 <- jags.samples(jm2,variable.names=c("a", "b","c", "mu2"), n.iter=n.iter, n.thin=10)
zj2

$a
mcarray:
[1] 36.08776 37.44723 88.92779

Marginalizing over: iteration(5000),chain(3) 

$b
mcarray:
[1] 2.211446

Marginalizing over: iteration(5000),chain(3) 

$c
mcarray:
[1] 5.785925

Marginalizing over: iteration(5000),chain(3) 

$mu2
mcarray:
            [,1]      [,2]        [,3]
 [1,] -0.6934448 -5.117440 -0.83128808
 [2,]  1.2850200 -2.645600 -0.38056993
 [3,]  1.2850200  1.701050  0.06545186
 [4,]  2.0950898  2.494504  1.37609112
 [5,]  3.6034248  3.436498  4.08230590
 [6,]  3.6034248  4.670523  4.08230590
 [7,]  5.4651579  5.335114  4.48330139
 [8,]  5.4651579  6.127841  5.07747127
 [9,]  5.4651579  6.281581  6.80907352
[10,]  8.4722087  6.281581  9.53625827
[11,]  8.7273492  6.433778 10.23218047
[12,]  9.5835708  7.875833 11.08456351
[13,]  9.5835708  9.311388 12.08260457
[14,] 10.7159355  9.311388 14.15706250
[15,] 12.5237726 11.162465 14.61951198
[16,] 12.7998976 12.098578 16.55729007
[17,] 12.7998976 12.969200 17.27521864
[18,] 13.3328650 13.954114 19.74299414
[19,] 15.3105391 15.174367 19.87499217
[20,] 15.8033581 15.627221 20.39786090
[21,] 16.2730532 15.846720 21.29361991
[22,] 16.2730532 16.945738 21.67021675
[23,] 16.4676764 17.576957 23.72573724
[24,] 18.1671637 17.637924 24.19085660
[25,] 18.3261492 17.878053 25.76670925
[26,] 18.6357232 18.674076 27.37048038
[27,] 19.5479765 19.255075 27.78421833
[28,] 19.7722726 19.407282 28.59491176
[29,] 20.4098084 20.713029 31.62996935
[30,] 20.9230665 20.799043 31.89865381
[31,] 20.9611491 21.456349 32.16478705
[32,] 21.0367392 21.495698 32.34081003
[33,] 21.1856591 21.953119 33.54256275
[34,] 21.2590115 22.314600 33.71001980
[35,] 23.5674193 23.894203 38.08522466
[36,] 23.8466090 24.436966 38.78134476
[37,] 23.8466090 24.488848 39.05441350
[38,] 24.0658982 24.591377 39.32448525
[39,] 24.0658982 24.816209 40.37577965
[40,] 24.1135614 24.913618 40.37577965
[41,] 24.1843484 25.009529 40.56789370
[42,] 24.3917529 25.056934 41.01022895
[43,] 24.5701112 25.127364 41.38290518
[44,] 24.9109258 25.311306 41.62811342
[45,] 25.0133151 25.533618 42.17058244
[46,] 25.0133151 25.831459 42.75874994
[47,] 26.2716216 26.818019 46.87956258
[48,] 26.7416377 27.250826 47.02327563
[49,] 26.7416377 27.500027 47.54184527
[50,] 26.7560491 27.649667 48.54107617
[51,] 27.1289844 27.837579 48.71754755
[52,] 27.1422237 28.032027 49.28048757
[53,] 27.2208427 28.112967 49.45051544
[54,] 27.2338113 28.126322 49.57707429
[55,] 27.2853078 28.218742 49.91059822
[56,] 28.0771618 28.963609 51.57211471
[57,] 28.2936330 29.104764 52.53107575
[58,] 28.4315497 29.147254 53.71273317
[59,] 28.7201913 29.403294 54.11155586
[60,] 28.7380664 29.699614 54.27506068
[61,] 29.0781678 30.042717 55.55967683
[62,] 29.4902416 30.417810 57.51808626
[63,] 29.6513808 30.507557 57.78289985
[64,] 29.6988702 30.659213 58.42570196
[65,] 29.8564949 30.722195 58.89667158
[66,] 30.5457234 31.451463 61.06911064
[67,] 30.5457234 31.462461 61.19445404
[68,] 30.5457234 31.467945 61.27738887
[69,] 30.5608780 31.489781 61.29804479
[70,] 30.5608780 31.527612 61.35982702
[71,] 30.5608780 31.543679 61.48256276
[72,] 30.5709351 31.591359 61.54351998
[73,] 30.5709351 31.591359 61.60420587
[74,] 30.5709351 31.596609 61.62437453
[75,] 30.5809555 31.601850 61.78465372
[76,] 30.5809555 31.607081 61.84427232
[77,] 30.5809555 31.617516 61.90362852

Marginalizing over: iteration(5000),chain(3) 

# summary(zj2$mu.a, quantile, c(0.025, 0.5, 0.975))$stat

# Plot parameters
# library(lattice) required for these plots
xyplot(zm2,ask = TRUE)

densityplot(zm2,ask = TRUE)

summary(zm2)

Iterations = 1501:6500
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
a[1]    36.025  3.0226 0.024679       0.050147
a[2]    37.496  3.0180 0.024642       0.049200
a[3]    88.507  7.3375 0.059910       0.188950
b        2.227  0.2594 0.002118       0.007085
c        5.821  0.9628 0.007861       0.018183
mu.a    53.640 18.0636 0.147489       0.217640
sigma.a 32.684  9.3573 0.076402       0.114217
sigma.o  5.037  0.5017 0.004096       0.015027
sigma.p  9.938  0.6476 0.005287       0.012503

2. Quantiles for each variable:

          2.5%    25%    50%    75%   97.5%
a[1]    30.431 33.953 35.910 37.946  42.322
a[2]    31.864 35.458 37.390 39.428  43.626
a[3]    75.860 83.351 88.009 92.969 104.407
b        1.751  2.047  2.214  2.391   2.769
c        3.700  5.238  5.905  6.493   7.462
mu.a    16.332 42.133 53.638 65.726  90.008
sigma.a 15.813 25.257 32.535 40.259  48.842
sigma.o  4.034  4.709  5.034  5.372   6.034
sigma.p  8.699  9.506  9.922 10.363  11.234

summary(zm2)$stat

             Mean         SD    Naive SE Time-series SE
a[1]    36.024512  3.0225682 0.024679166     0.05014748
a[2]    37.495944  3.0180138 0.024641980     0.04920039
a[3]    88.506537  7.3374549 0.059910069     0.18894984
b        2.226947  0.2594296 0.002118234     0.00708510
c        5.821200  0.9628198 0.007861391     0.01818302
mu.a    53.640447 18.0636442 0.147489038     0.21764002
sigma.a 32.684102  9.3573275 0.076402259     0.11421733
sigma.o  5.037443  0.5017150 0.004096486     0.01502715
sigma.p  9.937884  0.6475765 0.005287440     0.01250329

# Plot predictions
bu2 <- summary(zj2$mu2,quantile,c(.025,.5,.976))$stat
bu2

, , 1

            [,1]      [,2]      [,3]      [,4]      [,5]      [,6]     [,7]
2.5%  -5.2719791 -2.489593 -2.489593 -1.396183 0.5511916 0.5511916 2.859725
50%   -0.6228312  1.329798  1.329798  2.133757 3.6260280 3.6260280 5.469650
97.6%  3.4742256  4.877585  4.877585  5.470465 6.6075444 6.6075444 8.055556
          [,8]     [,9]     [,10]     [,11]     [,12]     [,13]     [,14]
2.5%  2.859725 2.859725  6.383341  6.675920  7.633815  7.633815  8.868874
50%   5.469650 5.469650  8.471259  8.728315  9.582602  9.582602 10.717759
97.6% 8.055556 8.055556 10.524780 10.744521 11.508493 11.508493 12.556054
         [,15]    [,16]    [,17]    [,18]    [,19]    [,20]    [,21]    [,22]
2.5%  10.80438 11.09946 11.09946 11.65422 13.64714 14.13410 14.59691 14.59691
50%   12.51809 12.79448 12.79448 13.32502 15.30447 15.79651 16.26504 16.26504
97.6% 14.28565 14.55561 14.55561 15.07728 17.05383 17.54878 18.03052 18.03052
         [,23]    [,24]    [,25]    [,26]    [,27]    [,28]    [,29]    [,30]
2.5%  14.78429 16.39204 16.54166 16.82660 17.66493 17.86388 18.43610 18.89076
50%   16.46136 18.15963 18.31777 18.62450 19.53191 19.75888 20.39808 20.91454
97.6% 18.23421 19.99748 20.16173 20.48989 21.48034 21.71924 22.40281 22.94557
         [,31]    [,32]    [,33]    [,34]    [,35]    [,36]    [,37]    [,38]
2.5%  18.92702 18.99561 19.12054 19.18510 21.20268 21.43460 21.43460 21.61806
50%   20.95274 21.02858 21.17845 21.25345 23.58478 23.86507 23.86507 24.08297
97.6% 22.98649 23.07022 23.23253 23.31226 25.92345 26.25004 26.25004 26.50783
         [,39]    [,40]    [,41]    [,42]    [,43]    [,44]    [,45]    [,46]
2.5%  21.61806 21.65786 21.71951 21.89593 22.04698 22.32561 22.40703 22.40703
50%   24.08297 24.13038 24.20292 24.41088 24.59031 24.93405 25.03593 25.03593
97.6% 26.50783 26.56239 26.64069 26.86431 27.06352 27.46354 27.58048 27.58048
         [,47]    [,48]    [,49]    [,50]    [,51]    [,52]    [,53]    [,54]
2.5%  23.41544 23.78905 23.78905 23.79943 24.08747 24.09820 24.15912 24.17023
50%   26.29806 26.77159 26.77159 26.78678 27.16362 27.17625 27.25465 27.26830
97.6% 29.07034 29.64403 29.64403 29.65943 30.11835 30.13337 30.22951 30.24511
         [,55]    [,56]    [,57]    [,58]    [,59]    [,60]    [,61]    [,62]
2.5%  24.20752 24.79730 24.97057 25.08430 25.30450 25.31959 25.57627 25.88614
50%   27.31979 28.11212 28.32941 28.46397 28.75149 28.76994 29.10374 29.51378
97.6% 30.31020 31.29275 31.56033 31.74134 32.10458 32.12843 32.56141 33.10538
         [,63]    [,64]    [,65]    [,66]    [,67]    [,68]    [,69]    [,70]
2.5%  26.01224 26.04982 26.16623 26.65807 26.65807 26.65807 26.66852 26.66852
50%   29.67586 29.72369 29.87549 30.56117 30.56117 30.56117 30.57663 30.57663
97.6% 33.31457 33.37647 33.58568 34.49683 34.49683 34.49683 34.51610 34.51610
         [,71]    [,72]    [,73]    [,74]    [,75]    [,76]    [,77]
2.5%  26.66852 26.67547 26.67547 26.67547 26.68240 26.68240 26.68240
50%   30.57663 30.58611 30.58611 30.58611 30.59633 30.59633 30.59633
97.6% 34.51610 34.53067 34.53067 34.53067 34.54414 34.54414 34.54414

, , 2

             [,1]      [,2]      [,3]       [,4]      [,5]     [,6]     [,7]
2.5%  -11.7469924 -8.052169 -1.942947 -0.8795375 0.3191769 1.852940 2.675408
50%    -4.9277810 -2.522880  1.737959  2.5236462 3.4621173 4.679768 5.342896
97.6%   0.6040506  2.154169  5.195910  5.7809243 6.5057354 7.456692 7.965973
          [,8]     [,9]    [,10]    [,11]     [,12]    [,13]    [,14]     [,15]
2.5%  3.622138 3.797076 3.797076 3.979608  5.664618  7.28750  7.28750  9.317829
50%   6.131888 6.285499 6.285499 6.437439  7.879624  9.31002  9.31002 11.154527
97.6% 8.612872 8.739334 8.739334 8.865400 10.075580 11.33747 11.33747 13.051463
         [,16]    [,17]    [,18]    [,19]    [,20]    [,21]    [,22]    [,23]
2.5%  10.32750 11.24581 12.24894 13.46411 13.92013 14.13181 15.20719 15.81106
50%   12.08320 12.95789 13.93981 15.16122 15.61572 15.83269 16.92532 17.55658
97.6% 13.94709 14.79230 15.76481 16.98084 17.44188 17.66757 18.80737 19.44824
         [,24]    [,25]    [,26]    [,27]    [,28]    [,29]    [,30]    [,31]
2.5%  15.87131 16.10207 16.86802 17.41013 17.54972 18.75156 18.83186 19.43019
50%   17.61793 17.85785 18.66010 19.24360 19.39676 20.70359 20.79076 21.44983
97.6% 19.51394 19.76779 20.59548 21.20729 21.36636 22.77330 22.87019 23.56703
         [,32]    [,33]    [,34]    [,35]    [,36]    [,37]    [,38]    [,39]
2.5%  19.46311 19.87217 20.19790 21.57810 22.03642 22.08121 22.16656 22.35338
50%   21.48928 21.94562 22.30785 23.89884 24.44568 24.49920 24.60095 24.82870
97.6% 23.60932 24.10286 24.50325 26.25128 26.87903 26.94154 27.06010 27.31354
         [,40]    [,41]    [,42]    [,43]    [,44]    [,45]    [,46]    [,47]
2.5%  22.43393 22.51557 22.55792 22.62078 22.77474 22.96444 23.20806 24.00745
50%   24.92452 25.01890 25.06594 25.13614 25.31906 25.54023 25.83788 26.82678
97.6% 27.42809 27.54242 27.59844 27.67812 27.89688 28.14930 28.48620 29.64140
         [,48]    [,49]    [,50]    [,51]    [,52]    [,53]    [,54]    [,55]
2.5%  24.35722 24.55341 24.67441 24.82844 24.98475 25.04744 25.05778 25.13258
50%   27.26268 27.51371 27.65964 27.84426 28.04050 28.12171 28.13474 28.22425
97.6% 30.17033 30.47223 30.65111 30.88723 31.13636 31.23664 31.25283 31.36301
         [,56]    [,57]    [,58]    [,59]    [,60]    [,61]    [,62]    [,63]
2.5%  25.72222 25.83262 25.86285 26.05552 26.28608 26.53847 26.81805 26.88763
50%   28.96192 29.10351 29.14705 29.39725 29.69334 30.03864 30.41197 30.50289
97.6% 32.29575 32.47236 32.52405 32.84577 33.22867 33.66469 34.12701 34.24539
         [,64]    [,65]    [,66]    [,67]    [,68]    [,69]    [,70]    [,71]
2.5%  26.99859 27.04917 27.58633 27.59536 27.59987 27.61692 27.64486 27.65635
50%   30.64974 30.71431 31.43688 31.44770 31.45285 31.47362 31.51090 31.52585
97.6% 34.44590 34.52287 35.49319 35.50705 35.51396 35.54289 35.59286 35.61423
         [,72]    [,73]    [,74]    [,75]    [,76]    [,77]
2.5%  27.69046 27.69046 27.69421 27.69797 27.70209 27.70957
50%   31.57242 31.57242 31.57749 31.58268 31.58791 31.59861
97.6% 35.68000 35.68000 35.68805 35.69441 35.70161 35.71596

, , 3

            [,1]      [,2]        [,3]      [,4]      [,5]      [,6]     [,7]
2.5%  -5.1036933 -4.543644 -3.98465603 -2.388722 0.7926723 0.7926723 1.245938
50%   -0.8422903 -0.390623  0.05932395  1.359577 4.0604418 4.0604418 4.457971
97.6%  3.5162909  3.880787  4.24662466  5.309897 7.5788644 7.5788644 7.935692
          [,8]     [,9]     [,10]    [,11]     [,12]     [,13]    [,14]
2.5%  1.913717 3.837045  6.772447  7.49474  8.389297  9.416832 11.53063
50%   5.048992 6.778182  9.513635 10.20186 11.057502 12.056800 14.13039
97.6% 8.448801 9.937466 12.415038 13.07181 13.882303 14.849237 16.88366
         [,15]    [,16]    [,17]    [,18]    [,19]    [,20]    [,21]    [,22]
2.5%  11.99633 13.94584 14.66768 17.10197 17.23352 17.74591 18.62548 18.99061
50%   14.59093 16.53993 17.25463 19.70538 19.84051 20.36201 21.26094 21.63851
97.6% 17.34424 19.28982 20.02910 22.55380 22.68749 23.21618 24.13426 24.52434
         [,23]    [,24]    [,25]    [,26]    [,27]    [,28]    [,29]    [,30]
2.5%  20.99244 21.44565 22.97951 24.52985 24.93612 25.72751 28.69809 28.95901
50%   23.69200 24.15861 25.73535 27.35025 27.76305 28.57034 31.61112 31.88112
97.6% 26.63476 27.11451 28.72344 30.34724 30.76422 31.57957 34.62353 34.89868
         [,31]    [,32]    [,33]    [,34]    [,35]    [,36]    [,37]    [,38]
2.5%  29.22086 29.39387 30.58497 30.75292 35.12182 35.81337 36.08450 36.35774
50%   32.14738 32.32081 33.52378 33.69351 38.08469 38.77826 39.05132 39.32200
97.6% 35.16541 35.34418 36.53739 36.70755 41.03808 41.73181 42.00288 42.26954
         [,39]    [,40]    [,41]    [,42]    [,43]    [,44]    [,45]    [,46]
2.5%  37.42234 37.42234 37.61580 38.06563 38.43329 38.68211 39.22793 39.81588
50%   40.37761 40.37761 40.57240 41.01512 41.38875 41.63345 42.17617 42.76797
97.6% 43.30988 43.30988 43.49953 43.93793 44.30850 44.55921 45.10716 45.69274
         [,47]    [,48]    [,49]    [,50]    [,51]    [,52]    [,53]    [,54]
2.5%  43.90840 44.05598 44.55594 45.52796 45.69439 46.22578 46.39482 46.52445
50%   46.87810 47.02396 47.54412 48.53925 48.71832 49.28256 49.45235 49.57956
97.6% 49.86494 50.01193 50.53659 51.56275 51.74722 52.32464 52.49897 52.62738
         [,55]    [,56]    [,57]    [,58]    [,59]    [,60]    [,61]    [,62]
2.5%  46.84152 48.44528 49.34118 50.43271 50.79631 50.94361 52.09416 53.78793
50%   49.91099 51.57000 52.53278 53.72067 54.12310 54.28491 55.57590 57.52924
97.6% 52.98521 54.74312 55.77797 57.06231 57.49233 57.67658 59.12748 61.36555
         [,63]    [,64]    [,65]    [,66]    [,67]    [,68]    [,69]    [,70]
2.5%  54.00437 54.53196 54.92727 56.72750 56.82524 56.89321 56.90966 56.96216
50%   57.79390 58.43670 58.90750 61.05876 61.18234 61.26394 61.28366 61.34531
97.6% 61.67499 62.40588 62.97725 65.55003 65.70580 65.79359 65.81460 65.88762
         [,71]    [,72]    [,73]    [,74]    [,75]    [,76]    [,77]
2.5%  57.07131 57.11665 57.16716 57.18284 57.31621 57.36721 57.41352
50%   61.46699 61.52925 61.58979 61.60936 61.76787 61.82907 61.89130
97.6% 66.03228 66.10135 66.16870 66.19318 66.38841 66.46161 66.53934

Additional resources

We have not used this package here, but do take a look at the functionality offered by the package ‘jagshelper’ - https://cran.r-project.org/web/packages/jagshelper/vignettes/jagshelper-vignette.html

install.packages("jagshelper")
library(jagshelper)

Plummer (2017) JAGS user manual - https://people.stat.sc.edu/hansont/stat740/jags_user_manual.pdf