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)

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

 

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