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

Code 5.31 Simulated beta–binomial data in R

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

# Required packages
require(R2jags)
require(boot)
require(VGAM)

# Simulation
set.seed(33559)
nobs = 2500

m = 1 + rpois(nobs,5)
x1 = runif(nobs)


beta1 <- -2
beta2 <- -1.5


eta <- beta1 + beta2 * x1
sigma <- 20

p <- inv.logit(eta)

shape1 = sigma*p
shape2 = sigma*(1 - p)

y <- rbetabinom.ab(n=nobs, size=m, shape1=shape1, shape2=shape2)

bindata = data.frame(y=y,m=m,x1)

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

Code 5.32  Beta–binomial synthetic model in R using JAGS

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

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

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

sink("GLOGIT.txt")

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

    # Prior for theta
    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()

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

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

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

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

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

Code 5.32
GET SOURCE

Output on screen:

Inference for Bugs model at "GLOGIT.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

beta[1]               -1.948       0.063           -2.075           -1.827     1.008     1500

beta[2]               -1.595       0.129           -1.833           -1.322      1.004      950

sigma                20.539       3.528           14.671          28.577      1.013      380

deviance       3337.284      74.220      3188.404      3477.857      1.011       380

 

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 = 2741.5 and DIC = 6078.8

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

bottom of page