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 10.14

require(R2jags)
require(jagstools)

# Data
path_to_data =

# Prepare data to JAGS
N <- nrow(GC_dat)
x <- GC_dat\$MV_T
y <- GC_dat\$N_GC
X <- model.matrix(~ x, data=GC_dat)
K = ncol(X)

JAGS_data <- list(
Y = y,
X = X,
N = N,
K = K)

Code 10.15 Negative binomial model, in R using JAGS, for modeling the relationship between globular cluster  population and host galaxy visual magnitude

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

# Fit
model.NB <- "model{
# Diffuse normal priors betas
for (i in 1:K) { beta[i] ~ dnorm(0, 1e-5)}

# Prior for theta
theta ~ dgamma(1e-3,1e-3)

for (i in 1:N){
eta[i] <- inprod(beta[], X[i,])
mu[i] <- exp(eta[i])
p[i] <- theta/(theta+mu[i])
Y[i] ~ dnegbin(p[i],theta)

# Discrepancy
expY[i] <- mu[i] # mean
varY[i] <- mu[i] + pow(mu[i],2)/theta # variance
PRes[i] <- ((Y[i] - expY[i])/sqrt(varY[i]))^2
}

Dispersion <- sum(PRes)/(N-3)
}"

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

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

# Start JAGS
NB_fit <- jags(data = JAGS_data ,
inits = inits,
parameters = params,
model = textConnection(model.NB),
n.thin = 1,
n.chains = 3,
n.burnin = 3500,
n.iter = 7000)

# Output
# Plot posteriors

MyBUGSHist(out,c("Dispersion",uNames("beta",K),"theta"))

# Dump results on screen
print(NB_fit, intervals=c(0.025, 0.975), digits=3)

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

Anchor 1

Output on screen:

Inference for Bugs model at "3", fit using jags,

3 chains, each with 7000 iterations (first 3500 discarded)

n.sims = 10500 iterations saved

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

Dispersion           1.919         0.207         1.547             2.351     1.007     350

beta              -11.704         0.359      -12.401         -11.026     1.053       50

beta                -0.878        0.018         -0.913           -0.844     1.052       51

theta                    1.097         0.073          0.963             1.244     1.003     920

deviance       5191.717         2.584     5188.701       5198.450     1.015     190

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 = 3.3 and DIC = 5195.0

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

bottom of page