Bayesian Split-Plot Design with heterogenous Variances / Kruschke / JAGS / ANOVA

I try to analyze a split-plot design and try to get my head around the code. My problem is that we were only able to gather a convenience sample of patients grouped into 11 different disorders at 2 measurement times, so a 2x11 design with nominal predictors, from which one is a repeated measurement. The y-variable is metric.

The script Jags-Ymet-XnomSplitPlot-MnormalHom.R from Kruschke’s DBDA book for homogenous variances and a normally distributed y seems to work fine. But, since sample sizes for each group are very, very different and it may be expected that different groups reacted differently from pre to post, I think a robust estimation for y is reasonable and more importantly, to allow for the estimation of heterogenous variances for each group within each measurement time. I only found the code for homogenous variances for the split-plot design and therefore tried to understand the differences in the codes for simple and robust 2-factorial designs without repeated measurements. I can see that y is estimated with a t-distribution of course (no problem to adapt it for the split plot design I think), but I am not sure about the estimation for the heterogenous variances in the Jags-Ymet-Xnom2fac-MrobustHet.R script :

for ( j1 in 1:Nx1Lvl ) { for ( j2 in 1:Nx2Lvl ) {
sigma[j1,j2] ~ dgamma( sigmaSh , sigmaRa )
# Prevent from dropping too close to zero:
ySigma[j1,j2] <- max( sigma[j1,j2] , medianCellSD/1000 )
} }
sigmaSh <- 1 + sigmaMode * sigmaRa
sigmaRa <- ( sigmaMode + sqrt( sigmaMode^2 + 4*sigmaSD^2 ) ) /(2*sigmaSD^2)
sigmaMode ~ dgamma(sGammaShRa,sGammaShRa)
sigmaSD ~ dgamma(sGammaShRa,sGammaShRa)

This is the only difference between the two scripts for homogenous and heterogenous variances, respectively.

What would be the analogous code for the split-plot design in the Jags-Ymet-XnomSplitPlot-MnormalHom.R script to allow for heterogenous variances? I am really not sure, since the partitioning of the variances in the split-plot design is quite different and I don’t want to mess it up. The model specification in the original split-plot script is the following:

for ( r in 1:Ntotal ) {
y[r] ~ dnorm( mu[r] , 1/sigma^2 )
mu[r] <- ( a0 + aB[xBetween[r]] + aW[xWithin[r]] + aS[xSubject[r]]
+ aBxW[xBetween[r],xWithin[r]] )
}
sigma ~ dunif( ySD/100 , ySD*10 )
a0 ~ dnorm( yMean , 1/(ySD*5)^2 )
for ( i in 1:NxBetweenLvl ) { aB[i] ~ dnorm( 0.0 , 1/sigmaB^2 ) }
sigmaB ~ dgamma(agammaShRa,agammaShRa)
for ( j in 1:NxWithinLvl ) { aW[j] ~ dnorm( 0.0 , 1/sigmaW^2 ) }
sigmaW ~ dgamma(agammaShRa,agammaShRa)
for ( k in 1:NxSubjectLvl ) { aS[k] ~ dnorm( 0.0 , 1/sigmaS^2 ) }
sigmaS ~ dgamma(agammaShRa,agammaShRa)
for ( i in 1:NxBetweenLvl ) { for ( j in 1:NxWithinLvl ) {
aBxW[i,j] ~ dnorm( 0.0 , 1/sigmaBxW^2 )
} }
sigmaBxW ~ dgamma(agammaShRa,agammaShRa)

I hope you can help me out and many, many thanks in advance for your effort.