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

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

 

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