1
$\begingroup$

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.

$\endgroup$

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.