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