top of page

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

library(MASS)

â€‹

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 <- rnegbin(mu, theta=2)

â€‹

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

â€‹

â€‹

â€‹

Code 8.22 Bayesian random intercept negative binomial in R using JAGS

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

library(R2jags)

â€‹

X <- model.matrix(~ x1 + x2, data = nbri)
K <- ncol(X)
Nre <- length(unique(nbri\$Groups))

â€‹

model.data <- list(
Y = nbri\$y,                                                 # response
X = X,                                                        # covariates
K = K,                                                        # num. betas
N = nrow(nbri),                                         # sample size
re = nbri\$Groups,                                      # random effects
b0 = rep(0,K),
B0 = diag(0.0001, K),
a0 = rep(0,Nre),
A0 = diag(Nre))

â€‹

sink("GLMM_NB.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)

â€‹

# Prior for alpha
numS ~ dnorm(0, 0.0016)
denomS ~ dnorm(0, 1)
alpha <- abs(numS / denomS)

â€‹

# Likelihood
for (i in 1:N) {
Y[i] ~ dnegbin(p[i], 1/alpha)
p[i] <- 1 /( 1 + alpha * mu[i])

log(mu[i]) <- eta[i]
eta[i] <- inprod(beta[], X[i,]) + a[re[i]]
}
}
"
,fill = TRUE)

â€‹

sink()

â€‹

# Define initial values
inits <- function () {
list(beta = rnorm(K, 0, 0.1),
a = rnorm(Nre, 0, 1),
num = rnorm(1, 0, 25),
denom = rnorm(1, 0, 1),
numS = rnorm(1, 0, 25) ,
denomS = rnorm(1, 0, 1))}

â€‹

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

â€‹

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

â€‹

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

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

Anchor 1

Output on screen:

â€‹

Inference for Bugs model at "GLMM_NB.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%         Rha        t n.eff

a[1]                -0.402          0.555        -1.112         0.516         7.815               3

a[2]                 0.358          0.545        -0.348         1.275         7.805               3

a[3]                 0.186          0.557        -0.527         1.124         7.588               3

a[4]                 0.805          0.537         0.088         1.679         7.675               3

a[5]                 0.692          0.548         -0.011         1.573         8.947              3

a[6]                -1.160          0.536         -1.777         -0.191       6.966              3

a[7]                 0.673          0.538         -0.026         1.576         8.766              3

a[8]                 0.586          0.545         -0.068        1.478         9.084               3

a[9]                 0.976          0.551         0.300         1.912         5.406               3

a[10]               0.505          0.541         -0.189         1.415         7.851              3

alpha               0.462         0.027           0.409         0.515        1.013           170

beta[1]            0.768          0.544         -0.135         1.444         8.943              3

beta[2]            0.172          0.071          0.028         0.299         1.006          300

beta[3]           -0.791          0.072         -0.926         -0.639       1.050            76

sigma.re          0.942          0.394         0.455         1.778         2.317               4

tau.re               1.771          1.286         0.316         4.838         2.317               4

deviance    7888.346          5.039   7880.194   7899.858         1.005           220

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 = 12.7 and DIC = 7901.0

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

bottom of page