How to specify nested model

I am using runjags to model some hierarchical data. I can model one level of the hierarchy but I do not know how to extend it to more levels. I am trying to do this using method 3 from page 24 of "Bayesian Hierarchical Modelling using WinBUGS", by Nicky Best et al which uses a nested loop (as opposed to nested indexing).

For one level I can model

filestring <-
  "model{ 
    for(j in 1:Ninner){
      for(i in 1:N){
        y[j,i] ~ dnorm(beta + alpha[j], py)
      }
      alpha[j] ~ dnorm(0, taua)
    }

  beta ~ dnorm(0, 0.001)
  taua ~ dgamma(0.01, 0.01)
  py ~ dgamma(0.01, 0.1)
}"

INITS <- list(list(.RNG.seed=1, .RNG.name="base::Wichmann-Hill"), 
              list(.RNG.seed=2, .RNG.name="base::Wichmann-Hill"))
results <- run.jags(filestring, monitor=c("py", "beta", "alpha"), data=jags_data, sample=1e3, 
                    n.chains=2, inits=INITS, summarise=FALSE)

I then tried to add another level using

for(k in 1:Nouter){
 for(j in 1:Ninner){
  for(i in 1:N){
    y[j,i] ~ dnorm(beta + alpha_in[j] + alpha_out[k], py)
} } }

but receive the error

Compilation error on line 5.
Attempt to redefine node y[1,1]

How do I extend this to model another level of which the first one is nested? Thank you.

Below is some reproducible data which shows the structure of the data. I wish to estimate random estimates for both outer_grp and the inner_grp.

library(data.table)
library(runjags)

set.seed(12345)
dat <- data.table(outer_grp=rep(1:5, each=10), inner_grp=rep(1:10, each=5), y=rnorm(50), x=rnorm(50), time=1:5)

wdat = dcast(dat, inner_grp + outer_grp ~ time, value.var=c("y", "x"))
jags_data = c(setNames(
  lapply(split.default(wdat, substr(names(wdat), 1, 1)),as.matrix), 
  c("inner_grp", "outer_grp","x", "y")),
  N=5, Nouter=5, Ninner=10)

EDIT

Perhaps it is enough to model like??

filestring <-
  "model{ 

     for(j in 1:Ninner){
      for(i in 1:N){
        y[j,i] ~ dnorm(beta + alpha_in[j] + alpha_out[outer_grp[j]], py)
      }
     }

  for(i in 1:Ninner){ alpha_in[i] ~ dnorm(0, taua) }
  for(i in 1:Nouter){ alpha_out[i] ~ dnorm(0, taub) }
  beta ~ dnorm(0, 0.001)
  taua ~ dgamma(0.01, 0.01)
  taub ~ dgamma(0.01, 0.01)
  py ~ dgamma(0.01, 0.1)
}"

It is possible to add the outer group intercept by using nested indexing while still using the loop format. I'll use the Pastes dataset from lme4 for comparison.

filestring <-
  "model{ 
     for(j in 1:Ninner){
      for(i in 1:N){
        y[j,i] ~ dnorm(beta + alpha_in[j] + alpha_out[batch[j]], py)
      }
     }

  for(i in 1:Ninner){ alpha_in[i] ~ dnorm(0, taua) }
  for(i in 1:Nouter){ alpha_out[i] ~ dnorm(0, taub) }
  beta ~ dnorm(0, 0.001)
  taua <- 1/(sa*sa)
  sa ~ dunif(0,100)
  taub <- 1/(sb*sb)
  sb ~dunif(0,100)
  py ~ dgamma(0.001, 0.001)
}"
INITS <- list(list(.RNG.seed=1, .RNG.name="base::Wichmann-Hill"), 
              list(.RNG.seed=2, .RNG.name="base::Wichmann-Hill"))
results <- run.jags(filestring, monitor=c("py", "beta", "alpha_in", "alpha_out", "sa", "sb"), 
                    data=jags_data, burnin=1e4, sample=1e4, n.chains=2, 
                    inits=INITS, summarise=0)
summary(results, vars=c("py", "beta", "sa", "sb"))

Compare to lme4

fm1 <- lmer(strength ~ (1|batch) + (1|sample), Pastes)
print(summary(fm1), corr=FALSE)

Data used

library(lme4); library(data.table); library(runjags)
data(Pastes); setDT(Pastes)
Pastes[,time := sequence(.N), by=sample]

# Change format to match question
wdat = dcast(Pastes, batch + sample ~ time, value.var="strength")
jags_data = list(y=as.matrix(wdat[,3:4]), batch=wdat$batch, N=2, Ninner=length(unique(wdat$sample)), Nouter=length(unique(wdat$batch)))