Logistic regression in JAGS

Question

I'm new to Bayesian analysis. I have a hierarchical model with a binary response variable. There is only one predictor (categorical), which has 3 levels: HLL, LHL and LLL. I prepared my data file by dummy-coding all these levels. My model specification is as follows:

cat("

model{
  for(i in 1:Ny){
    y[i] ~ dbern(mu[s[i]])
  }
for(j in 1:Ns){
mu[j] <- ilogit(b0[j] + b1*HLL[j] + b2*LHL[j] + b3*LLL[j])

b0[j] ~ dnorm(mu0, sigma0)
b1[j] ~ dnorm(mu1, sigma1)
b2[j] ~ dnorm(mu2, sigma2)
b3[j] ~ dnorm(mu3, sigma3)

  }

mu0 ~ dnorm(0, 0.001)
sigma0 ~ dunif(0, 100)

mu1 ~ dnorm(0, 0.001)
sigma1 ~ dunif(0, 100)

mu2 ~ dnorm(0, 0.001)
sigma2 ~ dunif(0, 100)

mu3 ~ dnorm(0, 0.001)
sigma3 ~ dunif(0, 100)
}
", fill = TRUE, file = "generalModel.txt")

Error

Basically, I want to get the estimates for HLL and LHL (using LLL as ref. level). This model doesn't run, and I'm not sure why. Here's the error message:

Calling 3 simulations using the parallel method...
Following the progress of chain 1 (the program will wait for all chains to finish
before continuing):
Welcome to JAGS 4.2.0 on Sun Jul 10 00:10:00 2016
JAGS is free software and comes with ABSOLUTELY NO WARRANTY
Loading module: basemod: ok
Loading module: bugs: ok
. . Reading data file data.txt
. Compiling model graph
   Resolving undeclared variables
   Allocating nodes
RUNTIME ERROR:
Invalid vector argument to ilogit
Deleting model
. Reading parameter file inits1.txt
Can't set RNG name. No model!
Can't set initial values. No model!
. Can't initialize. No model!
. Adaptation skipped: model is not in adaptive mode.
. Updating 1000
-------------------------------------------------| 1000
Can't update. No model!

. Can't set monitor. No model!
. Can't set monitor. No model!
. Can't set monitor. No model!
. Can't set monitor. No model!
. Updating 10000
-------------------------------------------------| 10000
Can't update. No model!

. No model
. Can't dump CODA output. No model!
. Can't dump samplers. No model!
. Updating 0
Can't update. No model!
Can't update. No model!
. Deleting model
. 
All chains have finished
Note: the model did not require adaptation
Error in runjags.readin(directory = startinfo$directory, silent.jags = silent.jags,  : 
  All the simulations appear to have crashed - check the model output in failed.jags() for clues
In addition: Warning messages:
1: No initial values were provided - JAGS will use the same initial values for all chains 
2: You attempted to start parallel chains without setting different PRNG for each chain, which is not recommended.  Different .RNG.name values have been added to each set of initial values. 
Note: Either one or more simulation(s) failed, or there was an error in processing
the results.  You may be able to retrieve any successful simulations using:
results.jags("/private/var/folders/nv/gznh75k93cv1wp35q1hvkkg00000gn/T/RtmpYRkQYd/runjagsfiles7c8d79109b99",
recover.chains=TRUE)
See the help file for that function for possible options.
To remove failed simulation folders use cleanup.jags() - this will be run
automatically when the runjags package is unloaded

Intercept-only model works fine

I successfully ran intercept-only models that are identical to the one above. In this case, I've run one model, say, using only LLL and another model using only HLL. Then I plotted the difference of both posteriors, and the result was pretty much in line with the estimate for HLL in a glmer() model where LLL was the ref. level.

cat("

model{
  for(i in 1:Ny){
    y[i] ~ dbern(mu[s[i]])
  }
  for(j in 1:Ns){
    mu[j] <- ilogit(b0[j])
    b0[j] ~ dnorm(mu0, sigma)
  }

  mu0 ~ dnorm(0, 0.001)
  sigma ~ dunif(0, 100)
}
", fill = TRUE, file = "modelA.txt")

Any ideas? Thanks!

Answer 1

Also (in extension of the previous answer), be careful with the parameterization of standard deviation in JAGS. dnorm is parameterized by precision, not standard deviation. It's much clearer (IMO) to write it like this:

b ~ dnorm( mu , 1/sigma^2 ) # Notice 1/sigma^2, not sigma

sigma ~ dunif( 0 , 100 )

Usually the uniform prior is put on sigma, not on precision.

[Wanted to add this as a comment instead of an answer, but the system wouldn't let me.]

Answer 2

The important part of the output is this error:

RUNTIME ERROR:
Invalid vector argument to ilogit

This comes from the following part of your model:

mu[j] <- ilogit(b0[j] + b1*HLL[j] + b2*LHL[j] + b3*LLL[j])

The parameters b1, b2 and b3 are indexed by j within the loop below this, so b1*HLL[j] results in a vector (of length Ns), which ilogit can't handle. It is possible that you meant to index b1 etc by j within this line i.e.:

mu[j] <- ilogit(b0 + b1[j]*HLL[j] + b2[j]*LHL[j] + b3[j]*LLL[j])

Note that I have dropped the index on b0 to have a fixed intercept, otherwise the model might not converge as quickly, but you can always add this back in later once you have the rest of the model working. There may be other differences between this and what you intended, but without more information on what you are trying to do it is hard for me to say. For example, it is possible that you don't want to index b0, b1, b2 and b3 by j, and rather you want a single b0, b1, b2 and b3 defined outside the loop - this wouldn't be hierarchical though, which you say you want.

A couple of other points:

1) I would be very careful with priors of the following type:

sigma0 ~ dunif(0, 100)

This may be quite an informative prior for sigma0, and will tend to pull the parameter towards higher values.

2) Make sure you load the GLM module - this will allow block updating of the parameters which should improve convergence

3) This is a more minor point but is is more customary to write:

logit(mu[j]) <- ....

Than:

mu[j] <- ilogit(...)

Aside from the benefit of making this clearer that you are writing a GLM, it would also change the error message JAGS gives you to something that might give you more clues as to the cause of coding errors.