model {
 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)

sigI ~ dunif(0, 5)
tauI <- 1/(sigI*sigI)

for(j in 1:nnests){
a[j] ~ dnorm(0, tau)
}

for(i in 1:N){

b[i] ~ dnorm(0, tauI)  ## fit lognormal-Poisson???

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]) <- log(offset[i]) + alpha + beta[1]*SexParent[i] + beta[2]*FoodTreatment[i] +
                    beta[3]*ArrivalTime[i] + beta[4]*SexParent[i]*FoodTreatment[i] +
                    beta[5]*SexParent[i]*ArrivalTime[i] + a[nest[i]] + b[i]

}

}