Random intercept binomial logistic model in R using JAGS

From: Bayesian Models for Astrophysical Data, Cambridge Univ. Press

you are kindly asked to include the complete citation if you used this material in a publication

Code 8.9 Random intercept binomial logistic data in R

============================================

y <- c(6,11,9,13,17,21,8,10,15,19,7,12,8,5,13,17,5,12,9,10)
m <- c(45,54,39,47,29,44,36,57,62,55,66,48,49,39,28,35,39,43,50,36)
x1 <- c(1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0)
x2 <- c(1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0)
Groups <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20)
logitr <- data.frame(y,m,x1,x2,Groups)

Code 8.10 Random intercept binomial logistic model in using JAGS

===================================================

library(R2jags)

X <- model.matrix(~ x1 + x2, data = logitr)
K <- ncol(X)

re <- length(unique(logitr$Groups))
Nre <- length(unique(Groups))

model.data <- list(
Y = logitr$y, # response
X = X, # covariates
m = m, # binomial denominator
N = nrow(logitr), # sample size
re = logitr$Groups, # random effects
b0 = rep(0,K),
B0 = diag(0.0001, K),
a0 = rep(0,Nre),
A0 = diag(Nre))

sink("GLMM.txt")

cat("
model{
# Diffuse normal priors for regression parameters
beta ~ dmnorm(b0[], B0[,])

# Priors for random effect group
a ~ dmnorm(a0, tau * A0[,])
num ~ dnorm(0, 0.0016)
denom ~ dnorm(0, 1)
sigma <- abs(num / denom)
tau <- 1 / (sigma * sigma)

# Likelihood function
for (i in 1:N){
Y[i] ~ dbin(p[i], m[i])
logit(p[i]) <- eta[i]
eta[i] <- inprod(beta[], X[i,]) + a[re[i]]
}
}",fill = TRUE)

sink()

inits <- function () {
list(
beta = rnorm(K, 0, 0.1),
a = rnorm(Nre, 0, 0.1),
num = rnorm(1, 0, 25),
denom = rnorm(1, 0, 1))}

params <- c("beta", "a", "sigma")

LOGIT0 <- jags(data = model.data,
inits = inits,
parameters = params,
model.file = "GLMM.txt",
n.thin = 10,
n.chains = 3,
n.burnin = 4000,
n.iter = 5000)

print(LOGIT0, intervals=c(0.025, 0.975), digits=3)

===================================================

Anchor 1

GET SOURCE

Output on screen:

Inference for Bugs model at "GLMM.txt", fit using jags,

3 chains, each with 5000 iterations (first 4000 discarded), n.thin = 10

n.sims = 300 iterations saved

mu.vect sd.vect 2.5% 97.5% Rhat n.eff

a[1] -0.581 0.445 -1.401 0.262 1.017 97

a[2] -0.432 0.374 -1.075 0.341 1.002 300

a[3] -0.007 0.401 -0.806 0.788 1.007 300

a[4] -0.137 0.363 -0.890 0.634 1.003 300

a[5] 1.126 0.421 0.337 1.956 0.999 300

a[6] 0.623 0.371 -0.019 1.429 0.999 300

a[7] -0.062 0.415 -0.862 0.677 1.007 190

a[8] -0.568 0.374 -1.356 0.071 0.996 300

a[9] -0.015 0.355 -0.696 0.652 1.005 260

a[10] 0.172 0.375 -0.483 0.931 0.998 300

a[11] -0.608 0.427 -1.438 0.165 1.018 96

a[12] -0.024 0.425 -0.860 0.741 0.998 300

a[13] -0.191 0.407 -0.949 0.555 0.997 300

a[14] -0.559 0.441 -1.502 0.266 1.009 150

a[15] 0.884 0.418 0.100 1.738 1.028 69

a[16] 0.817 0.403 0.107 1.591 1.003 300

a[17] -0.429 0.482 -1.372 0.472 1.004 290

a[18] 0.066 0.380 -0.670 0.817 1.004 280

a[19] -0.141 0.399 -0.919 0.666 1.027 68

a[20] 0.086 0.392 -0.615 0.911 0.998 300

beta[1] -1.103 0.322 -1.721 -0.430 1.006 300

beta[2] 0.249 0.331 -0.413 0.876 1.003 300

beta[3] -0.295 0.376 -1.008 0.446 1.011 140

sigma 0.676 0.177 0.394 1.041 1.002 300

deviance 97.575 6.516 86.888 111.331 1.005 300

For each parameter, n.eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)

pD = 21.3 and DIC = 118.9

DIC is an estimate of expected predictive error (lower deviance is better).

HSI

HSI