model
      {
         tau.alpha ~ dgamma(0.001,0.001)
         alpha0 ~ dnorm( 0.0,1.0E-6)
         beta0 ~ dnorm( 0.0,1.0E-6)
         tau.beta ~ dgamma(0.001,0.001)
         for( i in 1 : N ) {
            alpha[i] ~ dnorm(alpha0,tau.alpha)
            beta[i] ~ dnorm(beta0,tau.beta)
            y0[i] ~ dnorm(mu0[i],tau)
            mu0[i] ~ dnorm(theta,psi)
         }
         for( j in 1 : T ) {
            for( i in 1 : N ) {
               Y[i , j] ~ dnorm(mu[i , j],tau)
               mu[i , j] <- alpha[i] + beta[i] * (t[i , j] - 6.5) + 
                  gamma * (mu0[i] - mean(y0[]))
            }
         }
         tau ~ dgamma(0.001,0.001)
         sigma <- 1 / sqrt(tau)
         gamma ~ dnorm( 0.0,1.0E-6)
         theta ~ dnorm( 0.0,1.0E-6)
         psi ~ dgamma(0.001,0.001)
      }