model {
for (i in 2:N) {
y[i] ~ dbin(theta[i], n[i])
logit(theta[i]) <- logit.theta[i]
logit.theta[i] ~ dnorm(mu, inv.omega.squared)
}
inv.omega.squared <- 1/pow(omega, 2)
omega ~ dunif(0, 10)
mu ~ dunif(-10, 10)
# Mixed predictions of centre 1:
logit.theta1.cv ~ dnorm(mu, inv.omega.squared)
# generate replicate log-odds:
logit(theta1.cv) <- logit.theta1.cv
# generate replicate deaths:
y1.cv ~ dbin(theta1.cv, n[1])
# use mid p-value:
P.mixed <- step(y1.cv - y[1] - 0.00001)
+ 0.5*equals(y1.cv, y[1])
}

Data:
list(N = 12, y=c(41,25,24,23,25,42,24,53,26,25,58,31),
n=c(143,187,323,122,164,405,239,482,195,177,581,301))

Inits:
list(mu = -2, omega = 0.1)

   node   mean   sd   MC error   2.5%   median   97.5%   start   sample
   P.mixed   0.00165   0.03965   4.267E-4   0.0   0.0   0.0   1001   10000
   omega   0.1996   0.1132   0.004138   0.01784   0.1857   0.4557   1001   10000
   y1.cv   16.77   5.38   0.06203   7.0   16.0   28.0   1001   10000