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[1],sGammaShRa[2]) 
    sigmaSD ~ dgamma(sGammaShRa[1],sGammaShRa[2])

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[1],agammaShRa[2])
    for ( j in 1:NxWithinLvl ) { aW[j] ~ dnorm( 0.0 , 1/sigmaW^2 ) }
    sigmaW ~ dgamma(agammaShRa[1],agammaShRa[2])
    for ( k in 1:NxSubjectLvl ) { aS[k] ~ dnorm( 0.0 , 1/sigmaS^2 ) }
    sigmaS ~ dgamma(agammaShRa[1],agammaShRa[2])
    for ( i in 1:NxBetweenLvl ) { for ( j in 1:NxWithinLvl ) {
      aBxW[i,j] ~ dnorm( 0.0 , 1/sigmaBxW^2 )
    } }
    sigmaBxW ~ dgamma(agammaShRa[1],agammaShRa[2])

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


Your Answer

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

Browse other questions tagged or ask your own question.