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.13 Negative binomial: direct parameterization using JAGS and dnegbin

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

library(R2jags)

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

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

# Model
sink("NBDGLM.txt")

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

    # Prior for alpha
    alpha ~ dunif(0.001, 5)

     # Likelihood function
    for (i in 1:N){
        Y[i] ~ dnegbin(p[i], 1/alpha)                                               # for indirect, (p[i], alpha)
        p[i] <- 1/(1 + alpha*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
    alpha = runif(0.00, 5)
  )
}

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

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

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

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

 

Output on screen:

Inference for Bugs model at "NBDGLM.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.398         0.124              3.159                3.644         1.001      3000

beta[1]                0.994         0.090               0.821               1.163         1.010        290

beta[2]                2.039         0.083               1.875               2.200         1.005        440

beta[3]               -1.622         0.139             -1.893              -1.356         1.007        440

deviance      11548.717        2.900       11545.148        11555.896         1.007        670

 

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 = 11552.9

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