MIMIC models with vbmimic: regularizing measurement and structure

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:

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:

## 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

fit <- vbmimic(Y, X, Q_A, Q_B)
fit
#> vbmimic: regularized VB MIMIC model
#>   18 items, 3 factors, 9 covariates, N = 500
#>   converged in 19 iterations (0.19 secs)
#>   ELBO: NA (not available for this model) 
#>   active unspecified loadings    (PIP >= 0.5): 12 of 36
#>   active structural coefficients (PIP >= 0.5): 11 of 27
#> 
#> Components: model, call, nobs, nitem, nfactor, converged, A, B, pi_A, pi_B, Q_A, Q_B, eta, Phi, Sig, U, V, rho, theta, A_var, B_var, iter, flag, time, ELBO, objective, objective_type, standardize, preprocess, path
#> Access them with $ as usual; see ?vbpm_fit.

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

round(fit$B, 2)
#>      [,1]  [,2]  [,3]  [,4]  [,5]  [,6] [,7]  [,8]  [,9]
#> [1,] 0.25  0.30  0.32  0.00  0.01  0.00 0.02  0.00 -0.01
#> [2,] 0.00  0.01 -0.18  0.39  0.34  0.30 0.02 -0.01  0.00
#> [3,] 0.03 -0.02  0.01 -0.14 -0.01 -0.01 0.38  0.30  0.30
round(fit$pi_B, 2)
#>      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
#> [1,] 1.00 1.00 1.00 0.03 0.03 0.03 0.04 0.03 0.03
#> [2,] 0.03 0.03 1.00 1.00 1.00 1.00 0.04 0.03 0.03
#> [3,] 0.05 0.04 0.03 1.00 0.03 0.04 1.00 1.00 1.00

Compare PIP-based selection with the generating design:

selected <- fit$pi_B >= .5
table(truth = B != 0, selected = selected)
#>        selected
#> truth   FALSE TRUE
#>   FALSE    16    2
#>   TRUE      0    9

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():

round(fit$A, 2)[1:6, ]
#>      [,1]  [,2]  [,3]
#> [1,] 0.64  0.00  0.00
#> [2,] 0.63  0.00  0.00
#> [3,] 0.62 -0.02 -0.01
#> [4,] 0.62  0.00  0.01
#> [5,] 0.59  0.01 -0.01
#> [6,] 0.71  0.01 -0.03
round(fit$Phi, 2)            # factor correlations, from the disturbances
#>       [,1]  [,2]  [,3]
#> [1,]  1.00 -0.05 -0.03
#> [2,] -0.05  1.00 -0.01
#> [3,] -0.03 -0.01  1.00

Missing responses

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

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

fitm <- vbmimic(Ym, X, Q_A, Q_B)
fitm$preprocess$n_missing
#> [1] 20

## structural and measurement recovery are essentially unaffected
max(abs(fitm$B - fit$B))
#> [1] 0.002613604
max(abs(fitm$A - fit$A))
#> [1] 0.005325698

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:

Xm <- X; Xm[1, 1] <- NA
vbmimic(Y, Xm, Q_A, Q_B)
#> Error: X contains NA. vbmimic() does not support missing covariates: X is conditioned on, not modelled. Impute or drop them beforehand.

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