model {
  for(i in 1:nbeer) {
    for(j in 1:nobs) {
      sodium[i,j] ~ dnorm(theta[i], sigma2inv)
    }
    theta[i] ~ dnorm(mu, tau2inv)
  }
  mu ~ dnorm(0, 1.0E-6)
  tau2 ~ dunif(0, 1000)
  tau2inv <- pow(tau, -2)
  tau <- sqrt(tau2)
  sigma2inv ~ dgamma(0.001, 0.001)
  sigma <- 1 / sqrt(sigma2inv)
}