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

Code 5.33 Explicitly given beta–binomial data in R

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

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)

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

Code 5.34 Beta–binomial model in R using JAGS  for explicitly given data and the zero trick

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

library(R2jags)

 

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

K <- ncol(X)

 

model.data <- list(Y = bindata$y,
                              N = nrow(bindata),
                              X =X,
                              K =K,
                              m = m,
                              Zeros = rep(0, nrow(bindata))
)

sink("BBL.txt")

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

 

    # Prior for sigma
    sigma ~ dunif(0, 100)

 

    C <- 10000
    for (i in 1:N){
        Zeros[i] ~ dpois(Zeros.mean[i])
        Zeros.mean[i] <- -LL[i] + C

 

        #mu[i] <- 1/(1+exp(-eta[i])) # can use for logit(mu[i]) below


        logit(mu[i]) <- max(-20, min(20, eta[i]))
        L1[i] <- loggam(m[i]+1) - loggam(Y[i]+1) - loggam(m[i]-Y[i]+1)
        L2[i] <- loggam(1/sigma) + loggam(Y[i]+mu[i]/sigma)
        L3[i] <- loggam(m[i] - Y[i]+(1-mu[i])/sigma) - loggam(m[i]+1/sigma)
        L4[i] <- loggam(mu[i]/sigma) + loggam((1-mu[i])/sigma)
        LL[i] <- L1[i] + L2[i] + L3[i] - L4[i]
        eta[i] <- inprod(beta[], X[i,])
    }
}
"
,fill = TRUE)

 

sink()

 

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

 

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

 

BBIN0 <- jags(data = model.data,
                          inits = inits,
                         parameters = params,
                         model.file = "BBL.txt",
                         n.thin = 3,
                         n.chains = 3,
                         n.burnin = 10000,
                         n.iter = 15000)

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

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

 

Output on screen:

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

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

    n.thin = 3 n.sims = 5001 iterations saved

 

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

beta[1]                  -1.222        0.416              -2.055               -0.375       1.001        4800

beta[2]                   0.226        0.478              -0.725                 1.161       1.001        5000

beta[3]                   0.433        0.470              -0.506                 1.349       1.001        5000

beta[4]                  -0.249        0.436              -1.143                 0.622       1.001        2700

sigma                     0.106        0.076               0.021                  0.316       1.001       5000

deviance      240079.042        4.471      240073.054       240089.825       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 = 10.0 and DIC = 240089.0

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