HSI
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)
===============================================
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).