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.16 Random-intercept–random-slopes Poisson data in R

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

set.seed = 1056                                              # set seed to replicate example


N <- 5000                                                      # 10 groups, each with 500 observations
NGroups <- 10

x1 <- runif(N)
x2 <- ifelse( x1<=0.500, 0, NA)
x2 <- ifelse(x1> 0.500, 1, x2)

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

a <- rnorm(NGroups, mean = 0, sd = 0.1)
b <- rnorm(NGroups, mean = 0, sd = 0.35)

eta <- 2 + 4 * x1 - 7 * x2 + a[Groups] + b[Groups]*x1
mu <- exp(eta)
y <- rpois(N, lambda = mu)

pric <- data.frame(
  y = y,
  x1 = x1,
  x2 = x2,
  Groups = Groups)

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

 

Code 8.17 Random-intercept–random-slopes Poisson model in R using JAGS

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

library(R2jags)

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

re <- as.numeric(pric$Groups)


NGroups <- length(unique(pric$Groups))

model.data <- list(Y = pric$y,
                              X = X,
                              N = nrow(pric),
                              b0 = rep(0, K),
                              B0 = diag(0.0001, K),
                              re = re,
                              a0 = rep(0, NGroups),
                              A0 = diag(1, NGroups))

sink("RICGLMM.txt")

cat("
model{
    #Priors
    beta ~ dmnorm(b0[], B0[,])
    a ~ dmnorm(a0[], tau.ri * A0[,])
    b ~ dmnorm(a0[], tau.rs * A0[,])

    tau.ri ~ dgamma( 0.01, 0.01 )
    tau.rs ~ dgamma( 0.01, 0.01 )
    sigma.ri <- pow(tau.ri,-0.5)
    sigma.rs <- pow(tau.rs,-0.5)

    # 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]] + b[re[i]] * X[i,2]
        }
    }
    "
,fill = TRUE)

sink()

# Initial values
inits <- function () {
    list(
        beta = rnorm(K, 0, 0.01),
        a = rnorm(NGroups, 0, 0.1),
        b = rnorm(NGroups, 0, 0.1))}

# Identify parameters
params <- c("beta", "sigma.ri", "sigma.rs","a","b")

# Run MCMC
PRIRS <- jags(data = model.data,
              inits = inits,
              parameters = params,
              model.file = "RICGLMM.txt",
              n.thin = 10,
              n.chains = 3,
              n.burnin = 3000,
              n.iter = 4000)

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

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

 

Output on screen:

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

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

    n.sims = 300 iterations saved

 

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

a[1]                   0.164        0.053            0.066                0.266          1.000      300

a[2]                   0.016        0.054            -0.101                0.116          1.006      300

a[3]                   0.075        0.052            -0.032                0.174          1.005      300

a[4]                  -0.135        0.054            -0.244              -0.028          1.002      300

a[5]                   0.019        0.050            -0.074                0.121          1.007      300

a[6]                 -0.166         0.053            -0.276              -0.074          0.998       300

a[7]                 -0.028         0.052            -0.129                0.071          1.007      300

a[8]                  0.138         0.050             0.036                  0.231        1.022       300

a[9]                 -0.135        0.055            -0.238              -0.028          1.002       300

a[10]                0.067        0.056            -0.039               0.177           0.999       300

b[1]                -0.270         0.172             -0.625               0.043          1.004       300

b[2]                -0.497         0.181             -0.850              -0.139          1.008       190

b[3]                   0.468       0.164              0.157                0.814          1.009       300

b[4]                   0.586      0.174              0.248                0.935          1.009       210

b[5]                   0.030      0.172             -0.315                0.338         1.000        300

b[6]                   0.611      0.166              0.313                0.969          1.023        270

b[7]                 -0.205      0.171              -0.520               0.150           1.007       200

b[8]                 -0.426      0.174              -0.773              -0.103          1.004       300

b[9]                 -0.426      0.172              -0.777              -0.115          0.999       300

b[10]                0.180      0.170              -0.187                0.524          1.003       300

beta[1]             2.062      0.045                1.973                2.155          1.002       300

beta[2]             3.899      0.154                3.612                4.203          1.005       290

beta[3]            -6.993      0.048               -7.084              -6.904          1.057          45

sigma.ri            0.139      0.038                0.089               0.225          1.001        300

sigma.rs           0.483       0.138                0.289               0.811          1.010        190

deviance  17137.256        6.483       17126.668       17150.339           1.008        210

 

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.0 and DIC = 17158.2

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