## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(vbpm)

## ----simulate-----------------------------------------------------------------
B <- matrix(0, 3, 9)
for (k in 1:3) B[k, (k * 3 - 2):(k * 3)] <- .3   # 3 covariates per factor

sim <- sim_lvm(N = 500, K = 3, J = 18, P = 9, b = B, phx = 0, rseed = 1)
Y <- sim$dat[, 1:18]        # items first ...
X <- sim$dat[, 19:27]       # ... covariates last

## ----design-------------------------------------------------------------------
## Q_A is an AZ (anchor-zero) design: each anchor is specified (1) on its own
## factor and fixed to zero on the other two
Q_A <- matrix(-1L, 18, 3)
for (k in 1:3) {
  a <- which(rep(1:3, each = 6) == k)[1:2]
  Q_A[a, ] <- 0L
  Q_A[a, k] <- 1L
}
Q_B <- matrix(-1L, 3, 9)     # structural selection is the question

## ----fit----------------------------------------------------------------------
fit <- vbmimic(Y, X, Q_A, Q_B)
fit

## ----structural---------------------------------------------------------------
round(fit$B, 2)
round(fit$pi_B, 2)

## ----recovery-----------------------------------------------------------------
selected <- fit$pi_B >= .5
table(truth = B != 0, selected = selected)

## ----measurement--------------------------------------------------------------
round(fit$A, 2)[1:6, ]
round(fit$Phi, 2)            # factor correlations, from the disturbances

## ----missing------------------------------------------------------------------
Ym <- Y
Ym[cbind(1:20, rep(1:4, each = 5))] <- NA
sum(is.na(Ym))

fitm <- vbmimic(Ym, X, Q_A, Q_B)
fitm$preprocess$n_missing

## structural and measurement recovery are essentially unaffected
max(abs(fitm$B - fit$B))
max(abs(fitm$A - fit$A))

## ----missing-x, error = TRUE--------------------------------------------------
try({
Xm <- X; Xm[1, 1] <- NA
vbmimic(Y, Xm, Q_A, Q_B)
})