model {

  ## PRIORS
  for (m in 1:5){            
    beta[m] ~ dnorm(0, 0.01)           # Linear effects
  }
  alpha ~ dnorm(0, .01)
  sigma ~ dunif(0, 5)
  tau <- 1/(sigma*sigma)
  psi ~ dunif(0, 1)
  for(j in 1:nnests){
     a[j] ~ dnorm(0, tau)
   }

   for(i in 1:N){
       SibNeg[i] ~ dpois(mu[i])
       mu[i] <- lambda[i]*z[i]+0.00001 ## hack required for Rjags -- otherwise 'incompatible' 
       z[i] ~ dbern(psi)
       log(lambda[i]) <- offset[i] + alpha + inprod(X[i,],beta) + a[nest[i]]
   }
}