model{
  for (i in 1:N){
    y[i] ~ dbin(pi[i], n[i])
    pi[i] <- exp(b0+epsilon[i])/(1+exp(b0+epsilon[i]))
    epsilon[i] ~ dnorm(gamma, tau)
    }

gamma <- 0
b0 ~ dnorm(2, 0.0001)
tau ~ dgamma(0.01, 0.01)

}