Search code examples
jagsrjagshierarchical-bayesian

Hierarchical model using JAGS for R shows wrong number of means

enter image description here

I have 2 columns of data, y and grp and I am trying to create a JAGS model which is shown above. grp is group and I have 5 groups. The following code is from here. I am using this code because the description there under the heading Model and Data looks almost like this hiererchical model.

But I get only one mu when I view the summary. There should be 5 mu's, one for each group. Can someone correct the code ? You may also point to a similar example available elsewhere and I could try to modify it. I am missing something in the code and I believe the code could be like the question but when I modify it like that I don't seem to get proper means even though there are 5 means.

Not sure if this question belongs in math stackexchange.

mod_string = " model {
   for (i in 1:length(y) {
    theta[i] ~ dnorm(mu[grp[i]], invTau2)
    y[i] ~ dnorm(theta[i], 1/sig)
  }
  mu ~ dnorm(0, 1e6)
  invTau2 ~ dgamma(1.0/2.0, 1.3/2.0)
  tau2 <- 1/invTau2

  invgamma2 ~ dgamma(1.0/2.0, 2.1/2.0)
  sig = 1/invgamma2
    } "


summary(mod_sim)

Iterations = 2001:52000
Thinning interval = 1 
Number of chains = 3 
Sample size per chain = 50000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

         Mean    SD  Naive SE Time-series SE
mu  5.639e-07 0.001 2.582e-06      2.582e-06
sig 1.570e+00 1.888 4.874e-03      7.068e-03
Solution

  • It looks like you are using nested indexing, where grp is a vector of length y and describes what group each element y belongs to. If each group needs it's own mean then you need to generate 5 theta priors. This is most easily done in it's own loop. Thus, your model string would look something like:

    mod_string <- " model {
   # loop for y
   for (i in 1:length(y)) {
     y[i] ~ dnorm(theta[grp[i]], sig)
   }
   # loop for theta
   for (g in 1:ntheta){
     theta[g] ~ dnorm(mu, tau2)
   }
   # other priors. JAGS uses precision so you need to use the reciprocal 
   # of 1e6 for the variance of mu.
   mu ~ dnorm(0, 0.000001)
   invTau2 ~ dgamma(1.0/2.0, 1.3/2.0)
   tau2 <- 1/invTau2

   invgamma2 ~ dgamma(1.0/2.0, 2.1/2.0)
   sig <- 1/invgamma2
} "

    Note that you also need to supply the scalar ntheta to the model, which would be the number of unique groups (in this case ntheta = 5). Now if you track theta you will end up with 5 estimates whereas mu has 1 estimate (it is the grand estimated mean of all groups).