model{
  for(i in 1:N1){
    #measurement equation model
    for(j in 1:P){
      y1[i,j]~dnorm(mu1[i,j],errorprec1[j])
      z1[i,j]~dinterval(y1[i,j], cutoffs) 
    }
    mu1[i,1]<-xi1[i]
    mu1[i,2]<-lam1[1]*xi1[i]
    mu1[i,3]<-lam1[2]*xi1[i]
    mu1[i,4]<-eta1[i]
    mu1[i,5]<-lam1[3]*eta1[i]
    mu1[i,6]<-lam1[4]*eta1[i]

  #structural equation model
    xi1[i]  ~ dnorm(0,phi)
    eta1[i] ~ dnorm(nu1[i],laterrprec1)
    nu1[i] <- gam1 * xi1[i]

  }# end of i
  for(i in 1:N2){
    #measurement equation model
    for(j in 1:P){
      y2[i,j]~dnorm(mu2[i,j],errorprec2[j])
      z2[i,j]~dinterval(y2[i,j], cutoffs)
    }
    # Equality constraints on lambdas
    mu2[i,1]<-xi2[i]
    mu2[i,2]<-lam1[1]*xi2[i]
    mu2[i,3]<-lam1[2]*xi2[i]
    mu2[i,4]<-eta2[i]
    mu2[i,5]<-lam1[3]*eta2[i]
    mu2[i,6]<-lam1[4]*eta2[i]

  #structural equation model
    xi2[i]  ~ dnorm(0,phi)
    eta2[i] ~ dnorm(nu2[i],laterrprec2)
    nu2[i] <- gam2 * xi2[i]

  }# end of i

  for(j in 1:4) {
    lam1[j] ~ dnorm(1,1)
  }
  gam1 ~ dnorm(0.8,1)
  gam2 ~ dnorm(0.8,1)

  for(j in 1:P){
    errorprec1[j]~dgamma(3,1)
    errorvar1[j]<-1/errorprec1[j]
    errorprec2[j]~dgamma(3,1)
    errorvar2[j]<-1/errorprec2[j]
  }

  laterrprec1 ~ dgamma(3,1)
  laterrvar1 <- 1/laterrprec1
  laterrprec2 ~ dgamma(3,1)
  laterrvar2 <- 1/laterrprec2

  phi ~ dgamma(1,1)
} 

