• Simple mediation model data generation and estimation through re-parameterization

• A working example

## Simulate multilevel mediation model data

model{
  for (i in 1:n){
      re[i, 1:5]~dmnorm(mure[i, 1:5], pre_cov[1:5,1:5])
      mure[i,1]<-b[1,1]+b[1,2]*w[i]
      mure[i,2]<-b[2,1]+b[2,2]*w[i]
      mure[i,3]<-b[3,1]+b[3,2]*w[i]
      mure[i,4]<-b[4,1]+b[4,2]*w[i]
      mure[i,5]<-b[5,1]+b[5,2]*w[i]

     ## Generate w[i]
     w[i]~dnorm(0, 1)

      for (j in 1:mi[i]){
          y[i,j] ~ dnorm(muy[i,j], pre_phi[1])
          muy[i,j]<-re[j,1]+re[j,2]/re[j,3]*m[i,j]+re[j,4]*x[i,j]
          m[i,j] ~ dnorm(mum[i,j], pre_phi[2])
          mum[i,j]<-re[j,5]+re[j,3]*x[i,j]

      ## generate x[i,j]
      x[i,j]~dnorm(0, 1)
      }
  }
  for (i in 1:2){
     pre_phi[i]<-phi[i]
  }

  pre_cov[1:5,1:5]<-inverse(cov[1:5,1:5])
} 


## Population parameter values
list(n=100, phi=c(1,1), b=structure(.Data =c(
0.0,  1.0,  
0.0,  1.0,  
0.0,  1.0,  
0.0,  1.0,  
0.0,  1.0), .Dim = c(5,2)),
cov=structure(.Data =c(1, 0.5, 0.5, 0.5, 0.5, 0.5, 1, 0.5, 0.5, 0.5, 0.5, 
0.5, 1, 0.5, 0.5, 0.5, 0.5, 0.5, 1, 0.5, 0.5, 0.5, 0.5, 0.5, 1), .Dim = c(5,5)),
mi=c(100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 
100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 
100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 
100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 
100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 
100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 
100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 
100, 100, 100, 100, 100, 100, 100, 100, 100, 100))