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

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.11 Bayesian negative binomial in R using JAGS

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

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

â€‹

NB.data <- list(
Y = negbml\$nby,
N = nrow(negbml))

â€‹

library(R2jags)

# Attach(negbml)
X <- model.matrix(~ x1 + x2)
K <- ncol(X)

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

sink("NBGLM.txt")

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

â€‹

# Prior for dispersion
theta ~ dunif(0.001, 5)

â€‹

# Likelihood function
for (i in 1:N){
Y[i] ~ dpois(g[i])
g[i] ~ dgamma(theta, rateParm[i])
rateParm[i] <- theta / mu[i]
log(mu[i]) <- eta[i]
eta[i] <- inprod(beta[], X[i,])
}
}
"
,fill = TRUE)

â€‹

sink()

â€‹

inits <- function () {
list(
beta = rnorm(K, 0, 0.1),                                    # regression parameters
theta = runif(0.00, 5)                                         # dispersion
)}

â€‹

params <- c("theta""beta")

â€‹

NB2 <- jags(data = model.data,
inits = inits,
parameters = params,
model = "NBGLM.txt",
n.thin = 1,
n.chains = 3,
n.burnin = 3000,
n.iter = 5000)

â€‹

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

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

Anchor 1

Output on screen:

â€‹

Inference for Bugs model at "NBGLM.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

beta[1]              0.989          0.096             0.806                1.180        1.005          600

beta[2]             2.043           0.081             1.885                2.204        1.001        4700

beta[3]            -1.621           0.144           -1.898               -1.333        1.006          540

theta                 0.295           0.011            0.275                 0.316        1.001        6000

deviance     6609.402        62.560        6485.851         6732.318         1.001        6000

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 = 1957.5 and DIC = 8566.9

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

â€‹

bottom of page