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.16 Random-intercept–random-slopes Poisson data in R
=================================================
set.seed = 1056 # set seed to replicate example
N <- 5000 # 10 groups, each with 500 observations
NGroups <- 10
​
x1 <- runif(N)
x2 <- ifelse( x1<=0.500, 0, NA)
x2 <- ifelse(x1> 0.500, 1, x2)
​
Groups <- rep(1:10, each = 500)
​
a <- rnorm(NGroups, mean = 0, sd = 0.1)
b <- rnorm(NGroups, mean = 0, sd = 0.35)
​
eta <- 2 + 4 * x1 - 7 * x2 + a[Groups] + b[Groups]*x1
mu <- exp(eta)
y <- rpois(N, lambda = mu)
​
pric <- data.frame(
y = y,
x1 = x1,
x2 = x2,
Groups = Groups)
=================================================
​
​
Code 8.17 Random-intercept–random-slopes Poisson model in R using JAGS
===========================================================
library(R2jags)
​
X <- model.matrix(~ x1 + x2, data = pric)
K <- ncol(X)
re <- as.numeric(pric$Groups)
NGroups <- length(unique(pric$Groups))
​
model.data <- list(Y = pric$y,
X = X,
N = nrow(pric),
b0 = rep(0, K),
B0 = diag(0.0001, K),
re = re,
a0 = rep(0, NGroups),
A0 = diag(1, NGroups))
​
sink("RICGLMM.txt")
​
cat("
model{
#Priors
beta ~ dmnorm(b0[], B0[,])
a ~ dmnorm(a0[], tau.ri * A0[,])
b ~ dmnorm(a0[], tau.rs * A0[,])
tau.ri ~ dgamma( 0.01, 0.01 )
tau.rs ~ dgamma( 0.01, 0.01 )
sigma.ri <- pow(tau.ri,-0.5)
sigma.rs <- pow(tau.rs,-0.5)
​
# Likelihood
for (i in 1:N){
Y[i] ~ dpois(mu[i])
log(mu[i])<- eta[i]
eta[i] <- inprod(beta[], X[i,]) + a[re[i]] + b[re[i]] * X[i,2]
}
}
",fill = TRUE)
​
sink()
​
# Initial values
inits <- function () {
list(
beta = rnorm(K, 0, 0.01),
a = rnorm(NGroups, 0, 0.1),
b = rnorm(NGroups, 0, 0.1))}
​
# Identify parameters
params <- c("beta", "sigma.ri", "sigma.rs","a","b")
​
# Run MCMC
PRIRS <- jags(data = model.data,
inits = inits,
parameters = params,
model.file = "RICGLMM.txt",
n.thin = 10,
n.chains = 3,
n.burnin = 3000,
n.iter = 4000)
​
print(PRIRS, intervals=c(0.025, 0.975), digits=3)
===========================================================
Output on screen:
​
Inference for Bugs model at "RICGLMM.txt", fit using jags,
3 chains, each with 4000 iterations (first 3000 discarded), n.thin = 10
n.sims = 300 iterations saved
mu.vect sd.vect 2.5% 97.5% Rhat n.eff
a[1] 0.164 0.053 0.066 0.266 1.000 300
a[2] 0.016 0.054 -0.101 0.116 1.006 300
a[3] 0.075 0.052 -0.032 0.174 1.005 300
a[4] -0.135 0.054 -0.244 -0.028 1.002 300
a[5] 0.019 0.050 -0.074 0.121 1.007 300
a[6] -0.166 0.053 -0.276 -0.074 0.998 300
a[7] -0.028 0.052 -0.129 0.071 1.007 300
a[8] 0.138 0.050 0.036 0.231 1.022 300
a[9] -0.135 0.055 -0.238 -0.028 1.002 300
a[10] 0.067 0.056 -0.039 0.177 0.999 300
b[1] -0.270 0.172 -0.625 0.043 1.004 300
b[2] -0.497 0.181 -0.850 -0.139 1.008 190
b[3] 0.468 0.164 0.157 0.814 1.009 300
b[4] 0.586 0.174 0.248 0.935 1.009 210
b[5] 0.030 0.172 -0.315 0.338 1.000 300
b[6] 0.611 0.166 0.313 0.969 1.023 270
b[7] -0.205 0.171 -0.520 0.150 1.007 200
b[8] -0.426 0.174 -0.773 -0.103 1.004 300
b[9] -0.426 0.172 -0.777 -0.115 0.999 300
b[10] 0.180 0.170 -0.187 0.524 1.003 300
beta[1] 2.062 0.045 1.973 2.155 1.002 300
beta[2] 3.899 0.154 3.612 4.203 1.005 290
beta[3] -6.993 0.048 -7.084 -6.904 1.057 45
sigma.ri 0.139 0.038 0.089 0.225 1.001 300
sigma.rs 0.483 0.138 0.289 0.811 1.010 190
deviance 17137.256 6.483 17126.668 17150.339 1.008 210
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 = 21.0 and DIC = 17158.2
DIC is an estimate of expected predictive error (lower deviance is better).