# Multi-level Bayesian hierarchical regression using rjags

I am trying to to implement a Bayesian hierarchical Model in R. I have a few predictor variables (2 metric and one categorical) and am trying to predict quarterly home sales in the US. Each sales observation is for a state that is a member of a region, which in turn make up the population as a whole. I am using the rjags package.

I want to implement the model so that the coefficients for individual states are drawn from distributions determined by their region, which in turn are determined by the population. Essentially, I want to model this as a three-level hierarchy. However, I only understand how to implement the two-level hierarchy, i.e., state-level coefficients are a function of the population.

I have prepared an example that highlights my question and (hopefully) provides a framework for an answer. Below is my code that generates dummy data and implements my model. I would appreciate any insight into how to accomlish this as well as other suggeestions on how to improve my code, model structure, or anything else.

# load rjags for modeling
library(rjags)

# data
Nregion <- nlevels(state.region)
Nstate <- length(state.abb)
Nqtr <- 20
N <- Nstate * Nqtr # number of obs
y <- runif(N, 4500, 5500) # regresand
x1 <- dplyr::lag(y); x1[is.na(x1)] <- 4500 # lag of y as covaritate
x2 <- rnorm(N, 100, Nqtr) # second covariates
a1 <- rep(seq(1, 4), Nstate * 5) #qtr of the year
s <- as.factor(rep(state.abb, Nqtr)) # state
r <- as.factor(rep(state.region, Nqtr)) # region

# model
modelstring = "
model {
for ( i in 1:N ) {
y[i] ~ dnorm(y.hat[i], tau.i)
y.hat[i] <- b0[s[i]] + b1[s[i]] * x1[i] + b2[s[i]] * x2[i] + b3[s[i], a1[i]]
}

for ( s in 1:Nstate ) {
b0[s] ~ dnorm(b0u, b0t)
b1[s] ~ dnorm(b1u, b1t)
b2[s] ~ dnorm(b2u, b2t)

for ( j in 1:Nqtr ) {
b3[s, j] ~ dnorm(0, tau.sq.alpha)
}

}

b0u ~ dnorm(0, .0001)
b0t ~ dgamma(0.001, 0.001)
b1u ~ dnorm(0, .0001)
b1t ~ dgamma(0.001, 0.001)
b2u ~ dnorm(0, .0001)
b2t ~ dgamma(0.001, 0.001)
tau.i ~ dgamma(0.001, 0.001)
tau.sq.alpha ~ dgamma(0.001, 0.001)
sigma.sq.alpha <- 1 / tau.sq.alpha

}
"
# write to file
writeLines(modelstring, con="model.txt")

# create jags model object
jags <- jags.model('model.txt',
data = list('x1' = x1,
'x2' = x2,
'a1' = a1,
'y' = y,
'N' = N,
's' = s,
'Nstate' = Nstate,
'Nqtr' = Nqtr),
n.chains = 4,


You want a distribution for each quarter (given a state), each state (given a region), and each region. That means you'll need at least some state parameters indexed by s (in your model b0, b1, b2), region parameters indexed by r (which I'll call c0, c1, c2), and global parameters (which I'll call d0, d1, d2). Your model should look something like this:

modelstring = "
model {
for ( i in 1:N ) {
y[i] ~ dnorm(y.hat[i], tau.i)
y.hat[i] <- b0[s[i]] + b1[s[i]] * x1[i] + b2[s[i]] * x2[i] + b3[s[i], a1[i]]
}

for ( s in 1:Nstate ) {
b0[s] ~ dnorm(c0[r[s]], b0t)
b1[s] ~ dnorm(c1[r[s]], b1t)
b2[s] ~ dnorm(c2[r[s]], b2t)

for ( j in 1:Nqtr ) {
b3[s, j] ~ dnorm(0, tau.sq.alpha)
}

}

for ( r in 1:Nregion ) {
c0[r] ~ dnorm(d0, c0t)
c1[r] ~ dnorm(d1, c1t)
c2[r] ~ dnorm(d2, c2t)
}

b0t ~ dgamma(0.001, 0.001)
b1t ~ dgamma(0.001, 0.001)
b2t ~ dgamma(0.001, 0.001)
c0t ~ dgamma(0.001, 0.001)
c1t ~ dgamma(0.001, 0.001)
c2t ~ dgamma(0.001, 0.001)
d0 ~ dnorm(0, .0001)
d1 ~ dnorm(0, .0001)
d2 ~ dnorm(0, .0001)
tau.i ~ dgamma(0.001, 0.001)
tau.sq.alpha ~ dgamma(0.001, 0.001)
sigma.sq.alpha <- 1 / tau.sq.alpha

}
"


where r[s] gives a region index, given a state index.