N <- 1000
P <- 6

xi <- rnorm(N, 0, 1)

x1 <- 1*xi + rnorm(N, 0, 0.1)
x2 <- 1.1*xi + rnorm(N, 0, 0.1)
x3 <- 0.9*xi + rnorm(N, 0, 0.1)

eta <- 0.8*xi + rnorm(N, 0, 0.1)

y1 <- 1*eta + rnorm(N, 0, 0.1)
y2 <- 1.1*eta + rnorm(N, 0, 0.1)
y3 <- 0.8*eta + rnorm(N, 0, 0.1)

y <- cbind(x1, x2, x3, y1, y2, y3)

# Missing data mechanism
p.miss <- x1/20 + x2/20 + x3/20 
p.miss <- replace(p.miss, p.miss<0, 0)
p.miss <- replace(p.miss, p.miss>1, 1)

for (i in 1:N) {
  for (j in 1:P) {
    if (rbinom(1, 1, p.miss[i]) == 1) {
      y[i,j] <- NA
    }
  }
}

R <- replace(is.na(y), FALSE, 0)

# Output to file
dump(c('N', 'P', 'y', 'R'), 'example10data.jags.txt')
dput(list(N=N, P=P, y=aperm(y), R=aperm(R)), 'example10data.openbugs.txt')
dump(c('y'), 'example10data.ml.txt')

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

library(lavaan)

model <- 'xi =~ x1 + x2 + x3
eta =~ y1 + y2 + y3
eta ~ xi'

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


