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

(c) 2017,  Joseph M. Hilbe, Rafael S. de Souza and Emille E. O. Ishida 

 

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)

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

 

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).