N2 <- 500
N1 <- 1000
cutoffs <- c(-200, -2.5, -1, 1, 2.5, 200)

# Generate first sample

xi1 <- rnorm(N1, 0, 1)

x11 <- 1*xi1 + rnorm(N1, 0, 0.1)
x21 <- 1.1*xi1 + rnorm(N1, 0, 0.1)
x31 <- 0.9*xi1 + rnorm(N1, 0, 0.1)

eta1 <- 0.8*xi1 + rnorm(N1, 0, 0.1)

y11 <- 1*eta1 + rnorm(N1, 0, 0.1)
y21 <- 1.1*eta1 + rnorm(N1, 0, 0.1)
y31 <- 0.8*eta1 + rnorm(N1, 0, 0.1)

x11c <- findInterval(x11, cutoffs)
x21c <- findInterval(x21, cutoffs)
x31c <- findInterval(x31, cutoffs)
y11c <- findInterval(y11, cutoffs)
y21c <- findInterval(y21, cutoffs)
y31c <- findInterval(y31, cutoffs)

y1 <- cbind(x11, x21, x31, y11, y21, y31)
z1 <- cbind(x11c, x21c, x31c, y11c, y21c, y31c)

# Generate second sample

xi2 <- rnorm(N2, 0, 1)

x12 <- 1*xi2 + rnorm(N2, 0, 0.1)
x22 <- 1.1*xi2 + rnorm(N2, 0, 0.1)
x32 <- 0.9*xi2 + rnorm(N2, 0, 0.1)

eta2 <- 1.2*xi2 + rnorm(N2, 0, 0.1)

y12 <- 1*eta2 + rnorm(N2, 0, 0.1)
y22 <- 1.1*eta2 + rnorm(N2, 0, 0.1)
y32 <- 0.8*eta2 + rnorm(N2, 0, 0.1)

x12c <- findInterval(x12, cutoffs)
x22c <- findInterval(x22, cutoffs)
x32c <- findInterval(x32, cutoffs)
y12c <- findInterval(y12, cutoffs)
y22c <- findInterval(y22, cutoffs)
y32c <- findInterval(y32, cutoffs)

y2 <- cbind(x12, x22, x32, y12, y22, y32)
z2 <- cbind(x12c, x22c, x32c, y12c, y22c, y32c)

P <- 6
dump(c('N1', 'N2', 'P', 'cutoffs', 'z1', 'z2'), 'example9data.jags.txt')
dump(c('y1', 'y2'), 'example9inits1.jags.txt')

dput(list(N1=N1, N2=N2, P=P, cutoffs=cutoffs, z1=aperm(z1), z2=aperm(z2)), 'example9data.openbugs.txt')
dput(list(y1=aperm(y1), y2=aperm(y2)), 'example9inits1.openbugs.txt')

z1 <- cbind(z1, rep(1,N1))
z2 <- cbind(z2, rep(2,N2))

data <- data.frame(rbind(z1, z2))
colnames(data) <- c('x1', 'x2', 'x3', 'y1', 'y2', 'y3', 'group')

dump(c('data'), 'example9data.ml.txt')

data <- dget('example9data.ml.txt')

library(lavaan)

model <- 'xi =~ x1 + c(lam1, lam1)*x2 + c(lam2, lam2)*x3
eta =~ y1 + c(lam3, lam3)*y2 + c(lam4, lam4)*y3
eta ~ xi'

summary(sem(model, data, group='group', meanstructure=FALSE), standardized=TRUE)


