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

N <- 2000                                                 # 10 groups, each with 200 observations
NGroups <- 10

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

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


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

eta <- 1 + 0.2 * x1 - 0.75 * x2 + a[Groups]
mu <- exp(eta)
y <- rpois(N, lambda = mu)

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

Code 8.14 Bayesian random intercept Poisson model in R using JAGS

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

library(R2jags)

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

re <- as.numeric(poir$Groups)
Nre <- length(unique(poir$Groups))

model.data <- list(
  Y = poir$y,                                                  # response
  X = X,                                                         # covariates
  N = nrow(poir),                                          # sample size
  re = poir$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] ~ dpois(mu[i])
        log(mu[i])<- 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))}

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

# Run MCMC
PRI <- 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(PRI, 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.816       0.214           0.356          1.232       1.002        300

a[2]                    -0.558       0.223          -1.060         -0.129      1.004        300

a[3]                    -0.339       0.223          -0.776          0.097      1.013        210

a[4]                      0.041      0.216          -0.444          0.476       1.003       300

a[5]                      0.531      0.218            0.062          0.956      1.009       290

a[6]                     0.073       0.218           -0.398          0.483      1.003       300

a[7]                    -0.708       0.226          -1.208         -0.305      1.007       210

a[8]                     0.266       0.216          -0.198           0.705      1.001       300

a[9]                    -0.419       0.227          -0.945           0.047      1.003       300

a[10]                   0.163       0.219          -0.294           0.595       1.003      300

beta[1]                0.883       0.218           0.460            1.329      1.002      300

beta[2]                0.231       0.054           0.125            0.328      1.014      120

beta[3]               -0.769       0.055          -0.885           -0.665     1.017        96

sigma.re              0.576       0.176           0.357            0.979      0.998      300

tau.re                  3.717       1.854            1.044            7.868      0.998      300

deviance       6604.581       4.876      6597.501      6616.895      1.028      100

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 = 11.7 and DIC = 6616.3

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