Bayesian random intercept binary 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

# Data from Code 8.5
N <- 4000 #20 groups, each with 200 observations
NGroups <- 20

x1 <- runif(N)
x2 <- runif(N)

Groups <- rep(1:20, each = 200)

a <- rnorm(NGroups, mean = 0, sd = 0.5)

eta <- 1 + 0.2 * x1 - 0.75 * x2 + a[Groups]
mu <- 1/(1+exp(-eta))

y <- rbinom(N, prob=mu, size=1)

logitr <- data.frame(
y = y,
x1 = x1,
x2 = x2,
Groups = Groups,
RE = a[Groups]
)

Code 8.8 Bayesian random intercept binary logistic model in R using JAGS

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

library(R2jags)

X <- model.matrix(~ x1 + x2, data = logitr)
K <- ncol(X)
re <- as.numeric(logitr$Groups)
Nre <- length(unique(logitr$Groups))

model.data <- list(
Y = logitr$y, # Response
X = X, # Covariates
K = K, # Num. betas
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.re * A0[,])
num ~ dnorm(0, 0.0016)
denom ~ dnorm(0, 1)

sigma.re <- abs(num / denom)
tau.re <- 1 / (sigma.re * sigma.re)

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

sink()

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

params <- c("beta", "a", "sigma.re", "tau.re")

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

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

Anchor 1

GET SOURCE

Output:

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

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

n.sims = 600 iterations saved

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

a[1] 0.204 0.188 -0.155 0.547 1.329 10

a[2] -0.088 0.207 -0.524 0.280 1.203 14

a[3] -0.044 0.184 -0.375 0.295 1.149 18

a[4] -0.144 0.181 -0.484 0.214 1.169 17

a[5] -0.754 0.180 -1.133 -0.375 1.253 13

a[6] 0.277 0.198 -0.083 0.658 1.207 15

a[7] 0.525 0.229 0.013 0.914 1.284 11

a[8] 0.112 0.184 -0.242 0.487 1.223 13

a[9] -0.544 0.167 -0.865 -0.221 1.236 13

a[10] -0.031 0.188 -0.420 0.329 1.302 11

a[11] 0.082 0.199 -0.316 0.443 1.129 20

a[12] -0.971 0.188 -1.314 -0.586 1.152 19

a[13] 0.703 0.193 0.379 1.125 1.230 13

a[14] 0.037 0.190 -0.319 0.414 1.107 25

a[15] 0.180 0.198 -0.189 0.568 1.164 17

a[16] 0.705 0.195 0.364 1.079 1.204 14

a[17] 0.578 0.199 0.192 0.943 1.235 13

a[18] 0.039 0.170 -0.294 0.388 1.213 14

a[19] -0.499 0.186 -0.869 -0.145 1.201 14

a[20] 0.070 0.176 -0.281 0.388 1.413 8

beta[1] 0.856 0.158 0.561 1.170 1.308 10

beta[2] 0.022 0.120 -0.214 0.236 1.004 410

beta[3] -0.582 0.119 -0.809 -0.358 1.008 290

sigma.re 0.515 0.103 0.360 0.759 1.010 180

tau.re 4.206 1.577 1.734 7.726 1.010 180

deviance 5032.369 6.691 5021.427 5046.705 1.002 600

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 = 22.4 and DIC = 5054.8

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

HSI

HSI