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 5.33

y <- c(6,11,9,13,17,21,8,10,15,19,7,12)
m <- c(45,54,39,47,29,44,36,57,62,55,66,48)
x1 <- c(1,1,1,1,1,1,0,0,0,0,0,0)
x2 <- c(1,1,0,0,1,1,0,0,1,1,0,0)
x3 <- c(1,0,1,0,1,0,1,0,1,0,1,0)

bindata <- data.frame(y, m, x1, x2, x3)

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

K <- ncol(X)

 

Code 5.35 - Beta–binomial model with inverse link in R using JAGS

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

library(R2jags)

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

glogit.data <- list(Y = bindata$y,
                             N = nrow(bindata),
                             X = X,
                             K = K,
                             m = m)

sink("BBI.txt")

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

    # Prior for sigma
    sigma ~ dgamma(0.01,0.01)

    for (i in 1:N){
        Y[i] ~ dbin(p[i],m[i])
        p[i] ~ dbeta(shape1[i],shape2[i])

        shape1[i]<-sigma*pi[i]
        shape2[i]<-sigma*(1-pi[i])

        logit(pi[i]) <- eta[i]
        eta[i]<-inprod(beta[],X[i,])
        }
    }"
,fill = TRUE)

sink()

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

# Identify parameters
params <- c("beta", "sigma")

BB1 <- jags(data = glogit.data,
                      inits = inits,
                      parameters = params,
                      model.file = "BBI.txt",
                      n.thin = 1,
                      n.chains = 3,
                      n.burnin = 6000,
                      n.iter = 10000)

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

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

 

Output on screen:

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

    3 chains, each with 10000 iterations (first 6000 discarded)

    n.sims = 12000 iterations saved

 

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

beta[1]           -1.264          0.365      -2.006       -0.559      1.001       4700

beta[2]            0.233          0.407      -0.563         1.058      1.004        930

beta[3]            0.453          0.409      -0.347         1.270      1.003       1100

beta[4]          -0.242           0.377      -0.969         0.534      1.001       5900

sigma           22.054         18.066       4.930       71.257      1.001      12000

deviance      60.807           5.394      52.332      73.132      1.001       12000

 

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 = 14.5 and DIC = 75.4

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