model{
  for(i in 1:N){
    #measurement equation model
    for(j in 1:P){
       y[i,j]~dnorm(mu[i,j],errorprec[j])
       z[i,j]~dinterval(y[i,j], cutoffs) 
    }
    mu[i,1]<-xi[i]
    mu[i,2]<-lam[1]*xi[i]
    mu[i,3]<-lam[2]*xi[i]
    mu[i,4]<-eta[i]
    mu[i,5]<-lam[3]*eta[i]
    mu[i,6]<-lam[4]*eta[i]

  #structural equation model
    xi[i]  ~ dnorm(0,phi)
    eta[i] ~ dnorm(nu[i],laterrprec)
    nu[i] <- gam * xi[i]

  }# end of i

  for(j in 1:4) {
    lam[j] ~ dnorm(1,1)
  }
  gam ~ dnorm(0.8,1)

  for(j in 1:P){
    errorprec[j]~dgamma(3,1)
    errorvar[j]<-1/errorprec[j]
  }

  laterrprec ~ dgamma(3,1)
  laterrvar <- 1/laterrprec

  phi ~ dgamma(1,1)
} 

