3

I'm working on fitting a multi-level logistic regression model with group level predictors. I am using JAGS via R. I am getting different behaviors when I fit the model with the runjags versus the R2Jags packages.

I've tried to write a reproducible example that shows the issue. Below, I simulate data from a binomial model, index the data to 8 plots and 2 blocks, and then fit a multi-level logistic regression to recover the success probabilities (b1 and b2) in the code below. Scroll to the bottom to see the summaries of the two fits.

My question is:

  1. Why are the posteriors from these two fits different? I am using the same data, a single model specification, and setting the random number generator before each. Why does the mean of the posteriors differ, and why are the Rhat values so different? 
# -------------------------------------------------------------------
# Loading required packages
# -------------------------------------------------------------------
library(rjags) 
library(R2jags)
library(MCMCvis)

Package version information:

jags.version()
[1] ‘4.3.0’

R2jags_0.5-7   MCMCvis_0.13.5 rjags_4-10

# -------------------------------------------------------------------
# Simulate data
# -------------------------------------------------------------------
set.seed(10)

N.plots = 8
N.blocks = 2
trials=400

n = rep(100,trials)
N=length(n)
plotReps=N/N.plots
blockReps=N/N.blocks

# Block 1
b1<-rep(c(.25,.75,.9,.1),each=plotReps)-.05
# Block 2
b2<-rep(c(.25,.75,.9,.1),each=plotReps)+.05

y = rbinom(trials, 100, p = c(b1,b2))

# vectors indexing plots and blocks
plot = rep(1:8,each=plotReps)
block = rep(1:2,each=blockReps)

# pass data to list for JAGS
data = list(
  y = y,
  n = n,
  N = length(n),
  plot = plot,
  block= block,
  N.plots = N.plots,
  N.blocks = N.blocks
)

# -------------------------------------------------------------------
# Code for JAGS model
# -------------------------------------------------------------------

modelString <- "model { 
  ## Priors

  # hyperpriors
  mu.alpha ~ dnorm(0, 0.0001)

  sigma.plot ~ dunif(0,100) 
  tau.plot <- 1 / sigma.plot^2

  sigma.block ~ dunif(0,100) 
  tau.block <- 1 / sigma.block^2

  # priors 
  for(i in 1:N.plots){     
    eps.plot[i]~dnorm(0,tau.plot)
  }

  for(i in 1:N.blocks){
    eps.block[i]~dnorm(0,tau.block)
  }

  # Likelihood
  for(i in 1:N){
    logit(p[i]) <- mu.alpha + eps.plot[plot[i]] + eps.block[block[i]]
    y[i] ~ dbin(p[i], n[i])

  }
}"

# -------------------------------------------------------------------
# Initial values
# -------------------------------------------------------------------
# set inits for rjags
inits = list(list(mu.alpha = 0,sigma.plot=2,sigma.block=2),
             list(mu.alpha = 0,sigma.plot=2,sigma.block=2),
             list(mu.alpha = 0,sigma.plot=2,sigma.block=2)) 

# set inits function for R2jags
initsFun<-function(){list(
  mu.alpha=0,
  sigma.plot=2,
  sigma.block=2
)}

# -------------------------------------------------------------------
# Set JAGS parameters and random seed
# -------------------------------------------------------------------
# scalars that specify the 
# number of iterations in the chain for adaptation
# number of iterations for burn-in
# number of samples in the final chain
n.adapt = 500
n.update = 5000
n.iterations = 1000
n.thin = 1
parsToMonitor = c("mu.alpha","sigma.plot","sigma.block","eps.plot","eps.block")

# -------------------------------------------------------------------
# Call to JAGS via rjags
# -------------------------------------------------------------------
set.seed(2)
# tuning (n.adapt)
jm = jags.model(textConnection(modelString), data = data, inits = inits,
                n.chains = length(inits), n.adapt = n.adapt)

# burn-in (n.update)
update(jm, n.iterations = n.update)

# chain (n.iter)
samples.rjags = coda.samples(jm, variable.names = c(parsToMonitor), n.iter = n.iterations, thin = n.thin)

# -------------------------------------------------------------------
# Call to JAGS via R2jags
# -------------------------------------------------------------------
set.seed(2)
samples.R2jags <-jags(data=data,inits=initsFun,parameters.to.save=parsToMonitor,model.file=textConnection(modelString),
                      n.thin=n.thin,n.chains=length(inits),n.burnin=n.adapt,n.iter=n.iterations,DIC=T)

# -------------------------------------------------------------------
# Summarize posteriors using MCMCvis
# -------------------------------------------------------------------
sum.rjags <- MCMCvis::MCMCsummary(samples.rjags,params=c("mu.alpha","eps.plot","sigma.plot","sigma.block","eps.block"))
sum.rjags

sum.R2jags2 <- MCMCvis::MCMCsummary(samples.R2jags,params=c("mu.alpha","eps.plot","sigma.plot","sigma.block","eps.block"))
sum.R2jags2

Here is the output from an rjags fit:

                     mean         sd         2.5%         50%       97.5% Rhat n.eff
mu.alpha      0.07858079 21.2186737 -48.99286669 -0.04046538 45.16440893 1.11  4063
eps.plot[1]  -1.77570813  0.8605892  -3.45736942 -1.77762035 -0.02258692 1.00  2857
eps.plot[2]  -0.37359614  0.8614370  -2.07913650 -0.37581522  1.36611635 1.00  2846
eps.plot[3]   0.43387001  0.8612820  -1.24273657  0.42332033  2.20253810 1.00  2833
eps.plot[4]   1.31279883  0.8615840  -0.38750596  1.31179143  3.06307745 1.00  2673
eps.plot[5]  -1.34317034  0.8749558  -3.06843578 -1.34747145  0.44451006 1.00  2664
eps.plot[6]  -0.40064738  0.8749104  -2.13233876 -0.41530587  1.37910977 1.00  2677
eps.plot[7]   0.36515253  0.8738092  -1.35364716  0.35784379  2.15597251 1.00  2692
eps.plot[8]   1.71826293  0.8765952  -0.01057452  1.70627507  3.50314147 1.00  2650
sigma.plot    1.67540914  0.6244529   0.88895789  1.53080631  3.27418094 1.01   741
sigma.block  19.54287007 26.1348353   0.14556791  6.68959552 93.21927035 1.22    94
eps.block[1] -0.55924545 21.2126905 -46.34099332 -0.24261169 48.81435107 1.11  4009
eps.block[2]  0.35658731 21.2177540 -44.65998407  0.25801739 49.31921639 1.11  4457

and here is the output from an R2jags fit:

                   mean         sd         2.5%         50%       97.5% Rhat n.eff
mu.alpha     -0.09358847 19.9972601 -45.81215297 -0.03905447 47.32288503 1.04  1785
eps.plot[1]  -1.70448172  0.8954054  -3.41749845 -1.70817566  0.08187877 1.00  1141
eps.plot[2]  -0.30070570  0.8940527  -2.01982416 -0.30458798  1.46954632 1.00  1125
eps.plot[3]   0.50295713  0.8932038  -1.20985348  0.50458106  2.29271214 1.01  1156
eps.plot[4]   1.37862742  0.8950657  -0.34965321  1.37627777  3.19545411 1.01  1142
eps.plot[5]  -1.40421696  0.8496819  -3.10743244 -1.41880218  0.25843323 1.01  1400
eps.plot[6]  -0.45810643  0.8504694  -2.16755579 -0.47087931  1.20827684 1.01  1406
eps.plot[7]   0.30319019  0.8492508  -1.39045509  0.28668886  1.96325582 1.01  1500
eps.plot[8]   1.65474420  0.8500635  -0.03632306  1.63399429  3.29585024 1.01  1395
sigma.plot    1.66375532  0.6681285   0.88231891  1.49564854  3.45544415 1.04   304
sigma.block  20.64694333 23.0418085   0.41071589 11.10308188 85.56459886 1.09    78
eps.block[1] -0.45810120 19.9981027 -46.85060339 -0.33090743 46.27709625 1.04  1795
eps.block[2]  0.58896195 19.9552211 -46.39310677  0.28183123 46.57874408 1.04  1769

Here are trace plots for mu.alpha from the 2 fits. First, from the rjags fit:

Trace plot for mu.alpha from rjags fit

Second, from the R2jags fit:

Trace plot for mu.alpha from R2Jags fit

gregor-fausto
  • 435
  • 2
  • 9
  • 2
    Have you tried increasing the number of reps in your posterior? Have you checked that they both converged? – Dason Feb 19 '20 at 17:03
  • 2
    It looks like your mu.alpha parameter is very poorly constrained by the data. You end up with the central 95% of its distribution ranging from [-49.3, 51.7] in the rjags case, and [-36.3, 46.1] in the R2jags case. The two different fits look like they are converging on similar posterior distributions. I think any difference is due to sampling variation in the MCMC. Setting the same random seed for each one will not guarantee the same behavior since the two packages may set up the call to JAGS in different ways. – qdread Feb 19 '20 at 17:04
  • I just added trace plots for mu.alpha, which I think suggest that they have not converged (?). – gregor-fausto Feb 19 '20 at 17:10
  • Yes, those trace plots certainly suggest that they haven't converged – qdread Feb 19 '20 at 17:11
  • I upped the number of reps to 10000, and the summaries of the posterior are more similar between the two. However, the trace plot for mu.alpha is still sampling a lot of parameter space, so perhaps that's the issue that mu.alpha is poorly constrained by the data. I was trying to understand what a random-effects model for a multilevel logistic should look like with JAGS; I was expecting this to be similar to glmer(cbind(y,n) ~ 1 + (1|plot) + (1|block),family="binomial"). – gregor-fausto Feb 19 '20 at 17:19
  • 1
    Thanks to both of you. I think my issue here is actually unrelated to rjags or R2jags... As @qdread suggested, my mu.alpha is poorly constrained by the data. I think this is because I only have 2 blocks. When I run this without `block` in the model, `mu.alpha` converges. – gregor-fausto Feb 19 '20 at 17:32

2 Answers2

0

While part of the issue is related to a lack of convergence for mu.alpha, another issue is how both packages determine the number of samples to collect from the posterior distribution. Additionally, the update call after jags.model should be:

update(jm, n.iter = n.update)

instead of

update(jm, n.iterations = n.update)

For rjags you can pretty easily specify the number of adaptation steps, update steps, and iteration steps. Looking at samples.rjags it is quite clear that each chain has a posterior of length n.iterations, for a total of (in this example) 3000 samples (n.iterations * n.chains). Conversely, R2jags::jags will sample the posterior a number of times equal to the n.iter argument minus the n.burnin argument. So, as you have specified this you have 1) not included the n.update steps into R2jags::jags and 2) only sampled the posterior a total of 1500 times (each chain only keeps 500 samples) compared to 3000 times from rjags.

If you wanted do a similar burn-in and sample the same number of times you could instead run:

samples.R2jags <-jags(
  data=data,
  inits=inits,
  parameters.to.save=parsToMonitor,
  model.file=textConnection(modelString),
  n.thin=n.thin,
  n.chains=length(inits),
  n.burnin=n.adapt + n.update ,
  n.iter=n.iterations +n.update + n.adapt,
  DIC=T
)

Finally, R2jags loads the glm module by default while rjags does not. This could potentially cause some discrepancies as the samplers that get used are likely different (at least in this case because you are fitting a glm). You could load the glm module with a rjags::load.module('glm') call before calling jags.model.

And while this is not related to the question, per se, I would avoid you i in your for loops within the model for each loop (use a different letter if the number of iterations varies among loops):

modelString <- "model { 
  ## Priors

  # hyperpriors
  mu.alpha ~ dnorm(0, 0.0001)

  sigma.plot ~ dunif(0,100) 
  tau.plot <- 1 / sigma.plot^2

  sigma.block ~ dunif(0,100) 
  tau.block <- 1 / sigma.block^2

  # priors 
  for(i in 1:N.plots){     
    eps.plot[i]~dnorm(0,tau.plot)
  }

  for(j in 1:N.blocks){
    eps.block[j]~dnorm(0,tau.block)
  }

  # Likelihood
  for(k in 1:N){
    logit(p[k]) <- mu.alpha + eps.plot[plot[k]] + eps.block[block[k]]
    y[k] ~ dbin(p[k], n[k])

  }
}"
mfidino
  • 3,030
  • 1
  • 9
  • 13
0

I'm pretty sure the reason your posteriors are different is because Jags doesn't care about a seed set in R code.

However! While the set.seed() does nothing directly for Jags and is useless when calling Jags directly via rjags, it does get propagated when you use R2Jags.

Let's compare:

  • rjags is a lower level interface. If you don't provide a choice of random generator and seeds in your inits Jags will base their initialisation is based on the current timestamp.
  • R2Jags wraps around rjags functions. The jags() (R2Jags) function calls jags.model() (rjags). If you check out the R-code of the jags() function you will see that it generates a seed for each chain using the runif() function in R. As the Jags seeds are relying on the output of the runif() function in R, setting a seed in R will ensure you get the same seeds for Jags every time you run.
Philip De Vries
  • 417
  • 4
  • 13