model {
for (i in 1: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.02225   0.1457   0.001438   0.0   0.0   0.0   1001   10000
   omega   0.4248   0.1297   0.001985   0.23   0.4059   0.7387   1001   10000
   y1.cv   19.35   8.917   0.08326   6.0   18.0   40.0   1001   10000