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

set.seed(18472)

nobs <- 750 

​

x1_2 <- rbinom(nobs,size=1,prob=0.7)
x2 <- rnorm(nobs,0,1)
xb <- 1 - 1.5*x1_2 - 3.5*x2

​

exb <- exp(xb)
py <- rpois(nobs, exb)
pois <- data.frame(py, x1_2, x2)


poi <- glm(py ~ x1_2 + x2, family=poisson, data=pois)

​

​

 

Code 6.4 Bayesian Poisson model in R using JAGS

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

require(R2jags)

​

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

​

model.data <- list(Y = pois$py,
                              X = X,
                              K = K,                                            # number of betas
                              N = nrow(pois))                             # sample size

​

sink("Poi.txt")

​

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

​

    for (i in 1:N) {
        Y[i] ~ dpois(mu[i])
        log(mu[i]) <- inprod(beta[], X[i,])
        }
    }"
,fill = TRUE)

​

sink()

 

inits <- function () {list(beta = rnorm(K, 0, 0.1))}

 

params <- c("beta")

 

POI <- jags(data = model.data,
                    inits = inits,
                    parameters = params,
                    model = "Poi.txt",
                    n.thin = 1,
                    n.chains = 3,
                    n.burnin = 4000,
                    n.iter = 5000)

​

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

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

Anchor 1

Output on screen:

​

Inference for Bugs model at "Poi.txt", fit using jags, 3 chains,

    each with 5000 iterations (first 4000 discarded)

    n.sims = 3000 iterations saved

 

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

beta[1]               1.005         0.028              0.983             1.026       1.001       3000

beta[2]             -1.507          0.008            -1.516            -1.497       1.002       3000

beta[3]             -3.500          0.012            -3.509            -3.492       1.002       3000

deviance      2443.973      170.754       2433.825       2444.234       1.006       3000

 

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)

bottom of page