Getting started with vbpm: partially confirmatory factor analysis

The design matrix Q

Everything starts with a J x K integer matrix with three codes:

A fully confirmatory model has no -1; a fully exploratory one is all -1 (plus enough anchors to identify the factors). Everything between is the partially confirmatory continuum.

Two anchor conventions recur throughout vbpm’s documentation, and both appear in this vignette. AO (anchor-only) codes each anchor’s intended cell 1 and leaves that anchor’s cells on the other factors -1: it asserts only where the anchor does load. AZ (anchor-zero) additionally fixes those other anchor cells to 0, asserting that an anchor of one factor has no loading on any other. AZ is the stronger claim. vignette("bifactor") compares the two side by side.

Simulate, design, fit

sim_fa() generates factor-analytic data from a loading pattern: K factors, ipf items per factor, primary loadings lam, cross-loadings lac.

sim <- sim_fa(N = 500, K = 3, ipf = 6, lam = .7, lac = .3, rseed = 1)
Y <- sim$dat

## an AZ (anchor-zero) design: two anchors per factor, each specified (1) on
## its own factor and fixed to zero on the other two; every non-anchor row is
## left entirely to the data
Q <- matrix(-1L, ncol(Y), 3)
for (k in 1:3) {
  a <- which(rep(1:3, each = 6) == k)[1:2]
  Q[a, ] <- 0L
  Q[a, k] <- 1L
}

fit <- vbfa(Y, Q)
fit
#> vbfa: VB partially confirmatory factor analysis
#>   18 items, 3 factors, N = 500  (oblique)
#>   converged in 112 iterations over 4 v0 stages (0.12 secs)
#>   ELBO: -11096.957 
#>   active unspecified loadings (PIP >= 0.5): 18 of 36
#> 
#> Components: model, call, nobs, nitem, nfactor, converged, Lam, pi, eta, Phi, PsiInv, ELBO, ELBO_conditional, objective, objective_type, rho, objective_terms, Lam_var, eta_cov, Phi_inv_mean, iter, flag, time, Q, Qe, orthogonal, ld, bifactor, n_general, general, preprocess, sample_cov, Psi, W, q_star, xi_star, tau, path
#> Access them with $ as usual; see ?vbpm_fit.

The six anchor rows are what makes this AZ rather than AO: they contribute six 1 cells and twelve fixed zeros, leaving only the 36 cells of the twelve non-anchor rows for the spike-and-slab prior to choose between — the denominator the summary reports. For the AO version, drop the Q[a, ] <- 0L line and those twelve cross-factor anchor cells stay -1.

Typing the fit’s name gives the compact summary above — the object carries an S3 class (vbpm_fit), but it is still an ordinary list underneath and every component is public API (see ?vbpm_fit):

round(fit$Lam[1:6, ], 2)   # posterior mean loadings
#>      [,1]  [,2] [,3]
#> [1,] 0.70  0.00    0
#> [2,] 0.69  0.00    0
#> [3,] 0.62 -0.01    0
#> [4,] 0.67  0.01    0
#> [5,] 0.66  0.34    0
#> [6,] 0.67 -0.29    0
round(fit$pi[1:6, ], 2)    # PIPs of the unspecified entries
#>      [,1] [,2] [,3]
#> [1,]    1 0.00 0.00
#> [2,]    1 0.00 0.00
#> [3,]    1 0.04 0.04
#> [4,]    1 0.05 0.04
#> [5,]    1 1.00 0.04
#> [6,]    1 1.00 0.04
round(fit$Phi, 2)          # factor correlations (oblique by default)
#>      [,1] [,2] [,3]
#> [1,] 1.00 0.53 0.52
#> [2,] 0.53 1.00 0.44
#> [3,] 0.52 0.44 1.00

A loading is conventionally treated as active when its PIP is at least .5. Compare recovered structure against the truth:

active <- (Q == 1) | (Q == -1 & fit$pi >= .5)
table(truth = sim$MLA != 0, active = active)
#>        active
#> truth   FALSE TRUE
#>   FALSE    30    0
#>   TRUE      0   24

The estimator is deterministic: it uses a fixed initialization and no random numbers, so there is no seed to set. The default v0 is a four-stage, warm-started regularization path. It is the right starting point for ordinary use; ?vbfa documents the advanced controls and the scalar fixed-spike form.

Fit statistics

fit_stats() computes SEM-like statistics for a fit, with hard selection (PIP >= tau) and the nominal parameter count as the defaults:

round(fit_stats(fit), 3)
#>      t_nom          t        t_S      RMSEA       SRMR        CFI        TLI 
#>     45.000     45.000     45.797      0.013      0.025      0.997      0.997 
#>        AIC        BIC      AIC_S      BIC_S       ELBO  objective 
#>  21849.901  22039.469  21851.169  22044.094 -11096.957 -11096.957 
#> attr(,"class")
#> [1] "vbpm_fit_stats" "numeric"       
#> attr(,"model")
#> [1] "vbfa"
#> attr(,"objective_type")
#> [1] "elbo"

t_nom and t_S are the nominal and soft counts; by default t = t_nom. This deterministic default does not depend on whether an optional package is installed. The function reads orthogonal and ld from the fit object itself, so a bifactor fit needs no extra argument — and supplying a contradictory one is an error, not a silent miscount.

For a specifically motivated sensitivity analysis, request the numerical Jacobian-rank count explicitly:

fit_stats(fit, rank_adjust = TRUE, rank_max_J = 100)

That branch requires the suggested package numDeriv. It stops rather than silently reverting to the nominal count when numDeriv is unavailable or when the number of items exceeds rank_max_J; raise the guard deliberately only when the computational cost is acceptable. Because the count can affect degrees of freedom, AIC, BIC, and derived fit indices, record the setting and use one policy throughout a comparison. pefa() exposes the same two arguments and computes the requested counts before discarding full candidate fits.

Factor-count evidence over a window (returned by pefa(), which applies no count rule itself) is covered in vignette("pefa"); bifactor and higher-order models are covered in vignette("bifactor").

Local dependence

Correlated residuals (e.g. testlets, shared stems) are handled by a graphical spike-and-slab prior on the residual precision, solved by QUIC (Jin, Chen, Yan, & Zhang, 2026). Turn it on with ld = TRUE; by default the search is fully exploratory, or restrict it with a J x J design Qe using the same -1/0/1 codes.

The LD calls below use a shorter iteration limit and a looser tolerance only to keep vignette rendering quick. For an analysis, omit those two arguments and start from the documented defaults.

simLD <- sim_fa(N = 500, K = 3, ipf = 6, lam = .7, lac = .3, ecr = .3,
                rseed = 2)
fLD <- vbfa(simLD$dat, Q, ld = TRUE, max_it = 300, tolVal = 1e-3)

## the largest recovered residual edges, vs the planted pairs
Poff <- abs(fLD$Psi); Poff[lower.tri(Poff, diag = TRUE)] <- 0
which(Poff >= sort(Poff, decreasing = TRUE)[3], arr.ind = TRUE)
#>      row col
#> [1,]   3   4
#> [2,]   8  13
#> [3,]   1  14
simLD$ofd_ind
#>      row col
#> [1,]  14   1
#> [2,]   7   2
#> [3,]   4   3
#> [4,]  13   8
#> [5,]  10   9
#> [6,]  16  15

All three of the largest estimated edges — items 3-4, 8-13, and 1-14 — are genuine planted pairs, so nothing spurious outranks a real residual correlation here. The other three planted pairs (2-7, 9-10, 15-16) have smaller estimated edges and fall outside the top three at this sample size: the graphical prior orders the residual dependencies it finds by strength, and reading only the top of that ordering will not exhaust them.

For LD fits, fit_stats() uses the estimated residual covariance W. The returned objective is the terminal VECM criterion; ELBO remains NA because the residual precision is point-updated rather than assigned a full variational distribution.

J <- ncol(simLD$dat)
Qe <- matrix(-1L, J, J)
Qe[1:3, 1:3] <- 1L    # a known testlet: freely estimated among these 3 items
Qe[4:6, 7:9] <- 0L    # a block known to be residually independent
Qe[7:9, 4:6] <- 0L
isSymmetric(unname(Qe))
#> [1] TRUE

fRestricted <- vbfa(simLD$dat, Q, ld = TRUE, Qe = Qe, max_it = 300,
                    tolVal = 1e-3)

round(fRestricted$Psi[1:3, 1:3], 3)   # freely estimated: off-diagonal is not forced
#>        [,1]   [,2]   [,3]
#> [1,]  3.063 -0.105 -0.095
#> [2,] -0.105  2.441  0.000
#> [3,] -0.095  0.000  2.843
round(fRestricted$Psi[4:6, 7:9], 3)   # fixed absent: driven to (numerical) zero
#>      [,1] [,2] [,3]
#> [1,]    0    0    0
#> [2,]    0    0    0
#> [3,]    0    0    0

## fully exploratory loadings isolate the comparison to the residual side
Qexp <- matrix(-1L, J, 3)

fDiag <- vbfa(simLD$dat, Qexp)
fLDc  <- vbfa(simLD$dat, Qexp, ld = TRUE, max_it = 300, tolVal = 1e-3)

round(rbind(diagonal = fit_stats(fDiag)[c("BIC", "RMSEA")],
            ld       = fit_stats(fLDc)[c("BIC", "RMSEA")]), 3)
#>               BIC RMSEA
#> diagonal 21360.78 0.129
#> ld       20643.80 0.065

fCtrl <- vbfa(simLD$dat, Q, ld = TRUE, max_it = 300, tolVal = 1e-3,
              ld_control = list(xi0 = c(0.1, 0.5, 1)))
fCtrl$converged   # TRUE if the final v0 stage met the tolerance
#> [1] TRUE

Missing data

vbfa() accepts NA values in continuous response matrices. At each iteration it updates all missing entries within a person jointly from their Gaussian conditional mean (given that person’s observed entries and the current factor-model fit) and adds the conditional covariance to the expected residual cross-product — the deterministic VB counterpart of the LAWBL/PCFA MCMC data augmentation. This in-loop treatment is valid under missing at random (MAR): missingness may depend on observed data (other items, covariates, or the current factor estimates) but not on the missing value itself. Missing-data handling under MAR in the (G)PCFA framework was established by Chen (2021). Entirely missing rows or items are rejected because their location or scale is not identified.

simM <- sim_fa(N = 400, K = 3, ipf = 6, lam = .7, lac = .3, rseed = 1)
Ym0  <- simM$dat

## items 5 and 11 go missing depending on the OBSERVED value of anchor items
## 1 and 7 (higher values make missingness more likely); the anchors
## themselves stay fully observed. This is MAR, not MCAR: the probability of
## missingness varies systematically with observed, not missing, data.
set.seed(42)
Ymar <- Ym0
p5  <- ifelse(Ym0[, 1] > stats::median(Ym0[, 1]), .40, .05)
p11 <- ifelse(Ym0[, 7] > stats::median(Ym0[, 7]), .40, .05)
Ymar[stats::rbinom(nrow(Ym0), 1, p5)  == 1, 5]  <- NA
Ymar[stats::rbinom(nrow(Ym0), 1, p11) == 1, 11] <- NA
sum(is.na(Ymar))
#> [1] 177

fMar <- vbfa(Ymar, Q)
fMar$preprocess$n_missing
#> [1] 177
round(fMar$Lam[c(1, 5, 7, 11), ], 2)   # loadings recovered despite the missingness
#>      [,1] [,2]  [,3]
#> [1,] 0.71 0.00  0.00
#> [2,] 0.63 0.35 -0.01
#> [3,] 0.00 0.73  0.00
#> [4,] 0.03 0.65  0.28

fClean <- vbfa(Ym0, Q)
max(abs(fMar$Lam - fClean$Lam))    # largest discrepancy
#> [1] 0.06340608
mean(abs(fMar$Lam - fClean$Lam))   # typical discrepancy is much smaller
#> [1] 0.006088052

data(nlsy27)
Yn <- as.matrix(nlsy27$dat)
dim(Yn)
#> [1] 3458   27
sum(is.na(Yn))            # incomplete cells, handled in-loop
#> [1] 1050

fn <- vbfa(Yn, nlsy27$Q)
fn
#> vbfa: VB partially confirmatory factor analysis
#>   27 items, 3 factors, N = 3458  (oblique)
#>   converged in 416 iterations over 4 v0 stages (11.72 secs)
#>   ELBO: -124045.2 
#>   active unspecified loadings (PIP >= 0.5): 24 of 74
#> 
#> Components: model, call, nobs, nitem, nfactor, converged, Lam, pi, eta, Phi, PsiInv, ELBO, ELBO_conditional, objective, objective_type, rho, objective_terms, Lam_var, eta_cov, Phi_inv_mean, iter, flag, time, Q, Qe, orthogonal, ld, bifactor, n_general, general, preprocess, sample_cov, Psi, W, q_star, xi_star, tau, path
#> Access them with $ as usual; see ?vbpm_fit.
fn$preprocess$n_missing   # recorded on the fit
#> [1] 1050
round(fn$Lam, 2)
#>        [,1]  [,2]  [,3]
#>  [1,]  0.65 -0.04 -0.01
#>  [2,]  0.64 -0.03 -0.01
#>  [3,] -0.06 -0.09 -0.01
#>  [4,]  0.46  0.01 -0.01
#>  [5,] -0.06 -0.27 -0.01
#>  [6,]  0.45  0.05 -0.01
#>  [7,]  0.52  0.01  0.03
#>  [8,] -0.29 -0.07  0.03
#>  [9,] -0.03  0.69 -0.02
#> [10,] -0.07  0.77 -0.02
#> [11,] -0.28  0.06 -0.01
#> [12,]  0.17  0.35  0.04
#> [13,] -0.39  0.00 -0.01
#> [14,] -0.28  0.01  0.04
#> [15,] -0.01  0.80 -0.01
#> [16,]  0.05  0.63  0.01
#> [17,]  0.06  0.65  0.02
#> [18,] -0.04  0.03  0.13
#> [19,]  0.03  0.19  0.05
#> [20,] -0.04  0.07  0.00
#> [21,]  0.03  0.24  0.14
#> [22,]  0.14  0.05  0.45
#> [23,]  0.17  0.07  0.35
#> [24,]  0.00 -0.03  0.90
#> [25,]  0.20  0.07  0.25
#> [26,] -0.01 -0.03  0.90
#> [27,] -0.32  0.01 -0.05

nrow(Yn) - sum(stats::complete.cases(Yn))   # respondents listwise deletion drops
#> [1] 379

From graded knowledge to a Q matrix

Prior knowledge is often graded rather than crisp — say, a membership score S[j, k] in [0, 1] for item j on factor k, from expert ratings, text similarity, or a clustering consensus. This kind of graded membership matrix is what the measurement literature calls a “soft” Q matrix. The illustrative analyst-chosen thresholds below turn such scores into a partially confirmatory (hard) Q; they are not package defaults:

• anchor (1) where the evidence is strong: S >= .80;
• leave to the data (-1) where it is equivocal: .20 <= S < .80;
• exclude (0) where it is absent: S < .20.

The three intervals partition [0, 1], so the code below is a direct transcription of the rule and every score lands in exactly one code.

## a toy graded membership matrix for the simulated items
set.seed(3)
S <- matrix(runif(18 * 3, 0, .2), 18, 3)          # baseline noise
S[cbind(1:18, rep(1:3, each = 6))] <- runif(18, .55, .95)  # true memberships

Qsoft <- matrix(0L, 18, 3)          # S < .20 stays 0
Qsoft[S >= .20 & S < .80] <- -1L
Qsoft[S >= .80] <- 1L
table(Qsoft)
#> Qsoft
#> -1  0  1 
#> 10 36  8

fsoft <- vbfa(Y, Qsoft)
fsoft$converged
#> [1] TRUE

This is simply one transparent way to build the hard -1/0/1 design that vbfa() accepts. The thresholds should come from the application, not from this example.

References

• Chen, J., Guo, Z., Zhang, L., & Pan, J. (2021). A partially confirmatory approach to scale development with the Bayesian Lasso. Psychological Methods, 26(2), 210–235. https://doi.org/10.1037/met0000293
• 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
• Chen, J. (2023). Fully and partially exploratory factor analysis with bi-level Bayesian regularization. Behavior Research Methods, 55(4), 2125–2142. https://doi.org/10.3758/s13428-022-01884-7
• Jin, Y., & Chen, J. (2025). Regularized variational approximation for partially confirmatory factor analysis. Structural Equation Modeling, 32(3), 437–449. https://doi.org/10.1080/10705511.2024.2432612
• Jin, Y., Chen, J., Yan, Z., & Zhang, Y. (2026). Sparse residual estimation in partially confirmatory factor analysis. PsyArXiv preprint. https://doi.org/10.31234/osf.io/dehtv_v2
• Chen, J., & Jin, Y. (2026). Recovering latent structures after variational Bayesian variable selection: Fit assessment and factor-number selection in partially exploratory factor analysis. arXiv preprint arXiv:2607.07159.
• 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