model {
for (i in 1:N) {
y[i] ~ dbin(theta, n[i])
prop[i] <- y[i]/n[i] # (extra 0.00001 avoids numerical errors if prop[i] = 0 or 1)
Ds[i] <- 2*n[i]*(prop[i]*log((prop[i]+0.00001)/theta)
+ (1-prop[i])*log((1-prop[i]+0.00001)/(1-theta)))
# sign of deviance residual
sign[i] <- 2*step(prop[i] - theta) - 1
dev.res[i] <- sign[i]*sqrt(Ds[i])
}
dev.sat <- sum(Ds[])
theta ~ dunif(0, 1)
}
Data:
list(N=12,
y=c(4,3,2,2,3,4,2,5,3,3,6,3),
n=c(14,19,32,12,16,41,24,48,20,18,58,30))
node mean sd MC error 2.5% median 97.5% start sample
dev.res[1] 1.644 0.237 0.002556 1.187 1.642 2.111 1001 10000
dev.res[2] 0.4679 0.2504 0.002702 -0.03402 0.4663 0.9584 1001 10000
dev.res[3] -1.116 0.2869 0.003094 -1.676 -1.116 -0.5581 1001 10000
dev.res[4] 0.4574 0.2008 0.002165 0.07027 0.4561 0.8507 1001 10000
dev.res[5] 0.7572 0.2363 0.002547 0.2995 0.7556 1.221 1001 10000
dev.res[6] -0.4866 0.3436 0.003706 -1.155 -0.4876 0.1902 1001 10000
dev.res[7] -0.6045 0.257 0.002771 -1.106 -0.6049 -0.106 1001 10000
dev.res[8] -0.3783 0.3763 0.004056 -1.108 -0.3805 0.3585 1001 10000
dev.res[9] 0.3791 0.2552 0.002746 -0.1211 0.3777 0.8781 1001 10000
dev.res[10] 0.5601 0.2459 0.002651 0.08607 0.5586 1.042 1001 10000
dev.res[11] -0.4335 0.4129 0.004445 -1.235 -0.4357 0.3759 1001 10000
dev.res[12] -0.3727 0.2953 0.003191 -0.9462 -0.3738 0.2074 1001 10000
dev.sat 7.487 1.419 0.01223 6.47 6.941 11.55 1001 10000
theta 0.1223 0.01811 1.95E-4 0.08918 0.1216 0.1596 1001 10000