---
title: "MIMIC models with vbmimic: regularizing measurement and structure"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{MIMIC models with vbmimic: regularizing measurement and structure}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

A MIMIC (multiple-indicators multiple-causes) model adds a structural part to
the factor model: observed covariates $X$ predict the factors, which in turn
produce the item responses $Y$. `vbmimic()` implements the extended MIMIC model
of Jin & Chen (2025, *Multivariate Behavioral Research*), whose contribution is
to regularize **both** parts: spike-and-slab priors select unspecified entries
of the measurement design `Q_A` (items on factors, `J x K`) *and* of the
structural design `Q_B` (factors on covariates, `K x P`). Both matrices use the
same `-1/0/1` codes as `vbfa()`.

With the default settings, `vbmimic()` reproduces the published estimator
exactly.

## Simulating MIMIC data

`sim_lvm()` (absorbed from the LAWBL package) generates the data. Supply the
structural coefficients as a matrix — a sparse design where each factor is
predicted by its own covariates is both realistic and estimable:

```{r 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
```

(A scalar `b` would make *every* covariate predict *every* factor, which
quickly leaves no disturbance variance — `sim_lvm()` stops with an informative
error if the design is infeasible.)

## Anchor at least one part

One modelling caution, learned the empirical way: with **both** `Q_A` and
`Q_B` fully exploratory the model converges cleanly but the solution is
rotationally ambiguous — factors can merge or swap. Anchoring either part
resolves it. Here we anchor two items per factor in the measurement part and
leave the whole structural part to the data:

```{r 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
```

AZ (anchor-zero) and AO (anchor-only) are the two anchor conventions `vbpm`'s
documentation uses throughout: AO codes each anchor's intended cell `1` and
leaves its cells on the other factors `-1`, while AZ additionally fixes those
cells to `0`. AZ is the stronger claim, and it is what anchors the measurement
part here; `vignette("bifactor")` compares the two.

## Fit and read

```{r fit}
fit <- vbmimic(Y, X, Q_A, Q_B)
fit
```

The structural question — *which covariates predict which factors?* — is
answered by `B` and its posterior inclusion probabilities `pi_B`:

```{r structural}
round(fit$B, 2)
round(fit$pi_B, 2)
```

Compare PIP-based selection with the generating design:

```{r recovery}
selected <- fit$pi_B >= .5
table(truth = B != 0, selected = selected)
```

Read that table both ways. All nine generating coefficients are selected — no
false negatives — but so are two of the eighteen true zeros, and not
marginally: both sit at PIP 1.00 in the matrix above (`B[2, 3]` and
`B[3, 4]`, each a factor picking up a covariate from a neighbouring block).
Exploratory structural selection at this `N` is therefore not free of false
positives, and a PIP at the ceiling is not on its own evidence that a
coefficient is real. Magnitude is the useful second filter here: the two false
positives are estimated at -0.18 and -0.14, below every one of the nine true
coefficients, whose smallest estimate is 0.25 against a generating value of
0.30.

The measurement side reads exactly as in `vbfa()`:

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

## Missing responses

`vbmimic()` accepts `NA` in `Y`; a handful of missing item responses is
enough to demonstrate it:

```{r 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))
```

The residual covariance here is diagonal (`vbmimic()` has no LD branch), so a
missing response is replaced in-loop by its conditional mean `eta_i A'`, with
the conditional variance `1 / V_j` carried into the residual sum of squares —
the same conditional-moment logic as `vbfa()`'s diagonal case, simplified
because there is no cross-item residual covariance to condition on. As with
`vbfa()`, this in-loop treatment is valid under **missing at random (MAR)**
(Chen, 2021).

`X` must be **complete**. Covariates are conditioned on, not modelled, so
there is no distribution to impute them from:

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

## Options worth knowing

* **`v0` path.** As in `vbfa()`, a decreasing vector gives a warm-started
  regularization path over both parts; the scalar default `0.001` is the
  published fixed spike.
* **`standardize`.** The default `FALSE` is the published estimator's
  behaviour. `vbfa()` by contrast always standardizes internally; set
  `standardize = TRUE` for the analogous behaviour.
* **Determinism.** Like `vbfa()`, the estimator consumes no random numbers —
  no seed argument, bit-identical reruns.
* **Fit statistics.** The published algorithm converges on a residual
  criterion, not an ELBO; `fit$ELBO` is `NA`. `fit_stats()` accepts a
  `vbmimic` fit (it dispatches to a dedicated method) but returns a
  deliberately limited result: `NA` for every SEM-like index (`RMSEA`, `BIC`,
  `ELBO`, `objective`, ...) and a single `n_active_coef` — the count of
  soft-selected measurement plus structural coefficients, not a full
  parameter count, and not to be fed to an information criterion. `pefa()`
  sweeps `vbfa()` from one backbone `Q0`; it does not sweep MIMIC models.

## References

* Jin, Y., & Chen, J. (2025). Regularized variational Bayesian approximations
  for variable selection in extended multiple-indicators multiple-causes
  models. *Multivariate Behavioral Research*.
  https://doi.org/10.1080/00273171.2025.2483253
* Chen, J. (2021). A generalized partially confirmatory factor analysis
  framework with mixed Bayesian Lasso methods. *Multivariate Behavioral
  Research*, 57(6), 879–894. https://doi.org/10.1080/00273171.2021.1925520
* Ročková, V., & George, E. I. (2018). The spike-and-slab LASSO. *Journal of
  the American Statistical Association*, 113(521), 431–444.
  https://doi.org/10.1080/01621459.2016.1260469