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 6.11

library(MASS)

â€‹

set.seed(141)
nobs <- 2500

â€‹

x1 <- rbinom(nobs,size=1, prob=0.6)
x2 <- runif(nobs)
xb <- 1 + 2.0*x1 - 1.5*x2

a <- 3.3

theta <- 0.303                                                                 # 1/a

exb <- exp(xb)

nby <- rnegbin(n=nobs, mu=exb, theta=theta)
negbml <- data.frame(nby, x1, x2)

â€‹

â€‹

â€‹

Code 6.14 Negative binomial with zero trick using JAGS directly

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

library(R2jags)

â€‹

X <- model.matrix(~ x1 + x2, data=negbml)

K <- ncol(X)
N <- nrow(negbml)

â€‹

model.data <- list(
Y = negbml
\$nby,
N =N,
K =K,
X =X,
Zeros = rep(
0, N)
)

â€‹

sink("NB0.txt")

cat("
model{
# Priors regression parameters
for (i in 1:K) { beta[i] ~ dnorm(0, 0.0001)}

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

C <- 10000
for (i in 1:N) {
# Log-likelihood function using zero trick:
Zeros[i] ~ dpois(Zeros.mean[i])
Zeros.mean[i] <- -L[i] + C

l1[i] <- 1/alpha * log(u[i])
l2[i] <- Y[i] * log(1 - u[i])
l3[i] <- loggam(Y[i] + 1/alpha)
l4[i] <- loggam(1/alpha)
l5[i] <- loggam(Y[i] + 1)
L[i] <- l1[i] + l2[i] + l3[i] - l4[i] - l5[i]
u[i] <- (1/alpha) / (1/alpha + mu[i])
log(mu[i]) <- max(-20, min(20, eta[i]))
eta[i] <- inprod(X[i,], beta[])
}
}
"
,fill = TRUE)

â€‹

sink()

inits1 <- function () {
list(beta = rnorm(K,
0, 0.1),
numS = rnorm(
1, 0, 25) ,
denomS = rnorm(
1, 0, 1)
) }

params1 <- c("beta", "alpha")

NB01 <- jags(data = model.data,
inits = inits1,
parameters = params1,
model =
"NB0.txt",
n.thin =
1,
n.chains =
3,
n.burnin = 3
000,
n.iter =
5000)

â€‹

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

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

Anchor 1

Output on screen:

â€‹

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

3 chains, each with 5000 iterations (first 3000 discarded)

n.sims = 6000 iterations saved

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

alpha                            3.393      0.123                     3.161                       3.644       1.001       6000

beta[1]                         0.979      0.090                     0.798                       1.154       1.001       5300

beta[2]                         2.044      0.082                     1.887                       2.207       1.002       1200

beta[3]                        -1.601      0.139                    -1.871                     -1.316       1.001       6000

deviance        50011548.685      2.906        50011545.105        50011555.966        1.000             1

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 = 4.2 and DIC = 50011552.9

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

bottom of page