model {
for (i in 1:N) {
      alphaTimeQ[i] ~ dnorm(alphaTimeQ.c[tr[i]],alphaTimeQ.tau)
      alphaTime[i] ~ dnorm(alphaTime.c[tr[i]],alphaTime.tau)
      alpha[i] ~ dnorm(alpha.c[tr[i]],alpha.tau)
}
for (j in 1:T) {
   for (i in 1:N) {
      mu[i , j] <- alpha[i] + (alphaTime[i] * (time[j] - tiBar)) + (alphaTimeQ[i] * pow(time[j] - tiBar,2))
   }
}

for (k in 1:K) {
  for (j in 1:T) {
     for (i in 1:N)  {
            Y[i , j , k] ~ dnorm(mu[i , j],tau.y)
            sresid[i , j , k] <- (Y[i , j , k] - mu[i , j]) * sqrt(tau.y)
            p.inv[i , j , k] <-  (sqrt(2 * PI / tau.y)) * exp(0.5 * pow(sresid[i , j , k],2))
     }
  }
}

tau.y ~ dgamma(0.001,0.001)
alphaTimeQ.c[1] ~ dnorm( 0.0,1.0E-6)
alphaTimeQ.c[2] ~ dnorm( 0.0,1.0E-6)
alphaTimeQ.c[3] ~ dnorm( 0.0,1.0E-6)
alphaTimeQ.c[4] ~ dnorm( 0.0,1.0E-6)
alphaTimeQ.tau ~ dgamma(0.001,0.001)
alphaTime.c[1] ~ dnorm( 0.0,1.0E-6)
alphaTime.c[2] ~ dnorm( 0.0,1.0E-6)
alphaTime.c[3] ~ dnorm( 0.0,1.0E-6)
alphaTime.c[4] ~ dnorm( 0.0,1.0E-6)
alphaTime.tau ~ dgamma(0.001,0.001)
alpha.c[1] ~ dnorm( 0.0,1.0E-6)
alpha.c[2] ~ dnorm( 0.0,1.0E-6)
alpha.c[3] ~ dnorm( 0.0,1.0E-6)
alpha.c[4] ~ dnorm( 0.0,1.0E-6)
alpha.tau ~ dgamma(0.001,0.001)
}