model
   {
      for(p in 1 : N) {
         Y[p] ~ dnorm(mu[p], tau[p])
         mu[p] <- alpha[school[p], 1] + alpha[school[p], 2] * LRT[p] 
            + alpha[school[p], 3] * VR[p, 1] + beta[1] * LRT2[p] 
            + beta[2] * VR[p, 2] + beta[3] * Gender[p] 
            + beta[4] * School.gender[p, 1] + beta[5] * School.gender[p, 2]
            + beta[6] * School.denom[p, 1] + beta[7] * School.denom[p, 2]
            + beta[8] * School.denom[p, 3]
         log(tau[p]) <- theta + phi * LRT[p]
         sigma2[p] <- 1 / tau[p]
         LRT2[p] <- LRT[p] * LRT[p]
      }
      min.var <- exp(-(theta + phi * (-34.6193))) # lowest LRT score = -34.6193
      max.var <- exp(-(theta + phi * (37.3807))) # highest LRT score = 37.3807

   # Priors for fixed effects:
      for (k in 1 : 8) { 
         beta[k] ~ dnorm(0.0, 0.0001) 
      }
      theta ~ dnorm(0.0, 0.0001)
      phi ~ dnorm(0.0, 0.0001)

   # Priors for random coefficients:
      for (j in 1 : M) {
         alpha[j, 1 : 3] ~ dmnorm(gamma[1:3 ], T[1:3 ,1:3 ]); 
         alpha1[j] <- alpha[j,1]
      }

   # Hyper-priors:
      gamma[1 : 3] ~ dmnorm(mn[1:3 ], prec[1:3 ,1:3 ]);
      T[1 : 3, 1 : 3 ] ~ dwish(R[1:3 ,1:3 ], 3)
   }