---
title: "Getting started with vbpm: partially confirmatory factor analysis"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Getting started with vbpm: partially confirmatory factor analysis}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

`vbpm` estimates psychometric measurement models by regularized mean-field
variational Bayes. Its models are *partially confirmatory*: instead of choosing
between a confirmatory analysis (every loading prespecified) and an exploratory
one (nothing prespecified), you specify what you know and let continuous
spike-and-slab priors select the rest from the data.

This vignette walks the factor-analysis side of the package: designing a `Q`
matrix, fitting with `vbfa()`, reading the fit, computing fit statistics with
`fit_stats()`, and local-dependence (residual-correlation) estimation.
Factor-count evidence over a window of `K` — which `pefa()` reports without
applying any count rule of its own — is covered in `vignette("pefa")`; bifactor
and higher-order structures are covered in `vignette("bifactor")`. It ends with
an empirical example, a look at missing data, and a recipe for turning *graded*
prior knowledge into a `Q` matrix.

## The design matrix Q

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

| code | meaning |
|------|---------|
| `1`  | **specified** (anchored): the loading is estimated freely, with no spike — you assert this item measures this factor |
| `0`  | **fixed zero**: the loading is constrained to zero — you assert it does not |
| `-1` | **unspecified**: the data decide, via a spike-and-slab prior; the fit reports a posterior inclusion probability (PIP) |

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`.

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

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

```{r access}
round(fit$Lam[1:6, ], 2)   # posterior mean loadings
round(fit$pi[1:6, ], 2)    # PIPs of the unspecified entries
round(fit$Phi, 2)          # factor correlations (oblique by default)
```

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

```{r recovery}
active <- (Q == 1) | (Q == -1 & fit$pi >= .5)
table(truth = sim$MLA != 0, active = active)
```

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:

```{r vbfit}
round(fit_stats(fit), 3)
```

`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:

```{r rank-adjustment, eval=FALSE}
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.

```{r ld}
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)
simLD$ofd_ind
```

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.

### Restricting the search: Qe

`Qe` mirrors `Q`'s three codes, but on the residual side: `1` frees an edge
(estimated, penalized only by the slab), `0` fixes it at zero, `-1` leaves it
to the spike-and-slab search. Suppose items 1-3 are a known testlet and a
specific block of items is known a priori to be residually independent —
everything else stays exploratory:

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

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
round(fRestricted$Psi[4:6, 7:9], 3)   # fixed absent: driven to (numerical) zero
```

The fixed-zero block comes back exactly zero; the fixed-one block is free to
take on whatever value the data support. This is the residual-side analogue of
`Q`: `1`/`0` remove an entry from the search entirely (anchored in or out), and
`-1` is the only code the spike-and-slab prior actually chooses between.

### Diagonal vs. local dependence: does modeling residual correlation matter?

Fit both a diagonal and an LD model to the same residually-correlated data and
compare with `fit_stats()`. It builds each fit's model-implied covariance from
its own estimated residual covariance (`W` under LD, the diagonal `1/PsiInv`
otherwise), so the `BIC`/`RMSEA` comparison below is fair even though the two
fits' `objective` values are not: `objective` is comparable only across fits
sharing one `objective_type` (`"elbo"` here, `"vecm"` under LD) — which is
exactly why the comparison uses the covariance-based statistics instead.

```{r ld-compare}
## 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)
```

On data simulated with planted residual correlations, the LD fit wins on both
counts, as it should.

### The ld_control knob

`ld_control` exposes the local-dependence path and penalty settings without
changing the call's shape. For instance, a shorter, coarser `xi0` path (the
spike-penalty schedule, in units of `N`):

```{r ld-control}
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
```

`$converged` is the readable form of `$flag`, the raw `1`/`0` convergence
indicator every fit also carries; both refer to the last `v0` stage, since
`max_it` is a cap per stage rather than per fit. The coarser `xi0` path still
converges here, which is exactly what to check whenever you shorten a path.

`diag_penalty` (whether the residual precision diagonal is penalized by `xi1`;
default `1`) lives in the same list, along with `xi1`, `quic_eps`,
`quic_max_it`, and the Beta prior parameters `a1`/`b1` on the LD proportion —
see `?vbfa` for the full set.

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

### A simulated illustration

To make the MAR mechanism concrete: two items are made missing with a
probability that depends on the *observed* value of an anchor item from a
different factor block — never on the missing item's own value.

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

fMar <- vbfa(Ymar, Q)
fMar$preprocess$n_missing
round(fMar$Lam[c(1, 5, 7, 11), ], 2)   # loadings recovered despite the missingness
```

Against a fit to the same data with no missingness at all, the loadings are
close, with the largest discrepancy well within the range expected from
losing part of two items' data:

```{r mar-sim-clean}
fClean <- vbfa(Ym0, Q)
max(abs(fMar$Lam - fClean$Lam))    # largest discrepancy
mean(abs(fMar$Lam - fClean$Lam))   # typical discrepancy is much smaller
```

### Empirical: NLSY 1997

`nlsy27` ships with the package: 3,458 respondents, 27 mixed-type items, and an
initial three-factor design with two to three anchors per factor. Its `Q` is an
AO design — the anchor cells are `1`, and every other cell, including the
anchors' cells on the other two factors, is `-1`.

The data are 1.12% incomplete, which exercises the missing-response path from
above: the `NA`s are passed straight to `vbfa()` and imputed in-loop, with no
listwise deletion.

```{r nlsy}
data(nlsy27)
Yn <- as.matrix(nlsy27$dat)
dim(Yn)
sum(is.na(Yn))            # incomplete cells, handled in-loop

fn <- vbfa(Yn, nlsy27$Q)
fn
fn$preprocess$n_missing   # recorded on the fit
round(fn$Lam, 2)
```

What to look for in that matrix: every anchored item returns on the factor it
was anchored to (items 1-2 on factor 1, 9-10 on factor 2, 22-24 on factor 3),
and the unspecified rows fill in around them: items 4, 6, and 7 join factor 1,
items 12 and 15-17 join factor 2, items 25-26 join factor 3.
The solution is genuinely sparse rather than uniformly loaded: items 3 and
18-21 stay under .25 in absolute value on all three factors, having found no
strong home, and a scattering of moderate negative entries (items 5, 8, 11, 13,
14, 27) marks items running opposite in direction to their factor's anchors.

One caveat remains, and it is about response type rather than missingness: 17
of the 27 items are polytomous and are treated here as continuous. A threshold
model for categorical and mixed responses is out of scope for this release
(see "Known limitations" in the README).

For comparison, the complete-case analysis discards every respondent with any
missing answer:

```{r nlsy-cc}
nrow(Yn) - sum(stats::complete.cases(Yn))   # respondents listwise deletion drops
```

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

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

fsoft <- vbfa(Y, Qsoft)
fsoft$converged
```

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