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 8.1 Random intercept Gaussian data generated in R

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

set.seed(1656)

N <- 4500                                                                       # 20 groups, each with 200 observations
NGroups <- 20

x1 <- runif(N)
x2 <- runif(N)

Groups <- rep(1:20, each = 225)

a <- rnorm(NGroups, mean = 0, sd = 0.5)
print(a, 2)

mu <- 1 + 0.2 * x1 - 0.75 * x2 + a[Groups]
y <- rnorm(N, mean=mu, sd=2)

normr <- data.frame(
y = y,
x1 = x1,
x2 = x2,
Groups = Groups,
RE = a[Groups]
)

â€‹

â€‹

â€‹

â€‹

â€‹

Code 8.2 - Random intercept normal model in R using JAGS

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

library(R2jags)

â€‹

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

re <- as.numeric(normr\$Groups)
Nre <- length(unique(normr\$Groups))

â€‹

model.data <- list(
Y = normr\$y,                                               # response
X = X,                                                          # covariates
N = nrow(normr),                                         # rows in model
re = re,                                                         # random effect
b0 = rep(0,K),                                              # parameter priors with initial 0 value
B0 = diag(0.0001, K),                                  # priors for V-C matrix
a0 = rep(0,Nre),                                           # priors for scale parameters
A0 = diag(Nre))                                           # hyperpriors for scale parameters

# Fit
sink("lmm.txt")

â€‹

cat("
model {
# Diffuse normal priors for regression parameters
beta ~ dmnorm(b0[], B0[,])

â€‹

# Priors for random intercept groups
a ~ dmnorm(a0, tau.plot * A0[,])

â€‹

# Priors for the two sigmas and taus
tau.plot <- 1 / (sigma.plot * sigma.plot)
tau.eps <- 1 / (sigma.eps * sigma.eps)

sigma.plot ~ dunif(0.001, 10)
sigma.eps ~ dunif(0.001, 10)

â€‹

# Likelihood
for (i in 1:N) {
Y[i] ~ dnorm(mu[i], tau.eps)
mu[i] <- eta[i]
eta[i] <- inprod(beta[], X[i,]) + a[re[i]]
}
}
"
,fill = TRUE)

â€‹

sink()

â€‹

inits <- function () {
list(beta = rnorm(K, 0, 0.01),
a = rnorm(Nre, 0, 1),
sigma.eps = runif(1, 0.001, 10),
sigma.plot = runif(1, 0.001, 10)
)}

â€‹

params <- c("beta","a", "sigma.plot", "sigma.eps")

â€‹

NORM0 <- jags(data = model.data,
inits = inits,
parameters = params,
model.file = "lmm.txt",
n.thin = 10,
n.chains = 3,
n.burnin = 6000,
n.iter = 10000)

â€‹

# Output
print(NORM0, intervals=c(0.025, 0.975), digits=3)

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

Anchor 1

Output on screen:

â€‹

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

3 chains, each with 10000 iterations (first 6000 discarded), n.thin = 10

n.sims = 1200 iterations saved

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

a[1]                    0.850           0.223               0.393                  1.292             1.000         1200

a[2]                    0.339           0.224              -0.106                  0.781            1.000          1200

a[3]                    0.027           0.221              -0.406                  0.477            1.000          1200

a[4]                    0.393           0.215              -0.028                  0.830            1.000          1200

a[5]                   -0.213           0.216             -0.659                   0.212            1.001          1200

a[6]                   -1.111           0.226              -1.579                 -0.678            1.001          1200

a[7]                   -0.946           0.214             -1.370                  -0.520            1.005          1200

a[8]                   -0.053           0.219             -0.471                  0.368             1.001           1200

a[9]                    0.448           0.220              0.025                   0.883             1.002           1200

a[10]                   0.347          0.221             -0.072                   0.792            1.000            1200

a[11]                  -0.558          0.214             -0.976                  -0.122            1.001            1200

a[12]                   -0.025          0.224            -0.498                   0.412            1.002              740

a[13]                   -1.002          0.219            -1.456                  -0.585            1.001            1200

a[14]                    1.224          0.222              0.784                   1.680            1.000             1200

a[15]                    1.265          0.228              0.829                   1.714            1.000             1200

a[16]                    0.051          0.214             -0.364                   0.450            1.000             1200

a[17]                   -0.408          0.221             -0.846                   0.017            1.001             1200

a[18]                   -0.020          0.222             -0.461                   0.429            1.001             1200

a[19]                   -1.300          0.216             -1.730                  -0.884            1.001            1200

a[20]                    0.600          0.220               0.192                   1.017            1.000            1200

beta[1]                 0.718          0.194               0.321                   1.102            1.000            1200

beta[2]                 0.204          0.104              -0.005                   0.402            1.004              470

beta[3]               -0.689          0.101               -0.889                 -0.490             1.002           1200

sigma.eps            1.994          0.021                1.953                  2.036             1.000            1200

sigma.plot           0.793          0.142                0.574                  1.122             1.003              620

deviance      18978.659         6.706         18967.426          18993.124             1.000            1200

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 = 22.5 and DIC = 19001.2

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

bottom of page