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.25 Synthetic data from a binomial model in R
====================================================

set.seed(33559)
nobs = 2000

â€‹

m = 1 + rpois(nobs,5)

x1 = runif(nobs)
x2 = runif(nobs)

â€‹

xb <- -2 - 1.5 * x1 + 3 * x2

exb <- exp(xb)

p <- exb/(1 + exb)                                                         # prob of p=0
y <- rbinom(nobs,prob=p,size=m)

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

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

â€‹

â€‹

Code 5.26 Binomial model in R using JAGS

====================================================
library(R2jags)

â€‹

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

â€‹

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

â€‹

sink("GLOGIT.txt")

â€‹

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

â€‹

# Likelihood
for (i in 1:N){
Y[i] ~ dbin(p[i],m[i])
logit(p[i]) <- eta[i]
eta[i] <- inprod(beta[], X[i,])
}
}
"
,fill = TRUE)

â€‹

sink()

â€‹

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

â€‹

params <- c("beta")

â€‹

BINL <- 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(BINL, intervals=c(0.025, 0.975), digits=3)

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

Anchor 1

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.964           0.177            -2.118                 -1.543                  1.003          1100

beta[2]              -1.558           0.166            -1.733                -1.228                   1.012            190

beta[3]               3.038           0.281              2.350                 3.251                   1.002          2100

deviance       4980.539       274.496       4945.780           5097.757                   1.001          6000

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 = 37686.2 and DIC = 42666.7

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

bottom of page