3
$\begingroup$

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,
                   n.adapt = 100)
$\endgroup$
4
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.