---
title: "Bifactor Structures and PEFA Sweeps"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Bifactor Structures and PEFA Sweeps}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

This vignette covers the bifactor side of `vbpm`: what a bifactor model is,
how a `pefa()` sweep is read when the factor count is uncertain, and how a
later oblique-versus-bifactor comparison is built out of the sweep that
preceded it. The sweep object itself is walked in `vignette("pefa")`; nothing
here re-explains its tables.

## Three related structures

Suppose `J` items fall into `K` clusters.

**Correlated factors.** Each item loads on its cluster factor, and the cluster
factors correlate through $\Phi$.

**Orthogonal bifactor.** Every item loads on a general factor and on one group
factor; general and group factors are mutually orthogonal.

**Higher-order.** Items load on first-order factors, which load in turn on a
second-order factor. Its Schmid-Leiman representation is the orthogonal
bifactor with

$$
b^{gen}_j=\lambda_j\gamma_k,
\qquad
b^{grp}_j=\lambda_j\sqrt{1-\gamma_k^2}
$$

for item $j$ in cluster $k$, so $|b^{grp}_j|/|b^{gen}_j|$ is constant within a
cluster under exact higher-order proportionality. Testlet and method factors
use the same general-plus-special geometry under a different substantive
reading.

In `vbfa()` and `pefa()`, `bifactor = TRUE` makes `K` the number of **group**
factors: the design matrix supplies only the group block and the general column
is added internally, so a fitted matrix has `K + 1` columns. Every
between-candidate comparison a sweep performs removes that general column
first and works on the `K` group columns, which is what makes an ordinary and a
bifactor sweep of the same data structurally comparable.

## A bifactor population and two backbones

Three clusters of six items, each item loading on a general factor and on its
own group factor.

```{r population}
K <- 3
ipf <- 6
J <- K * ipf
cluster <- rep(seq_len(K), each = ipf)

set.seed(7)
mla <- matrix(0, J, K + 1)
mla[, 1] <- round(runif(J, .35, .75), 2)
for (k in seq_len(K)) {
  mla[cluster == k, k + 1] <- round(runif(ipf, .35, .70), 2)
}
while (max(rowSums(mla^2)) > .9) mla <- round(mla * .97, 3)

sim <- sim_fa(N = 800, mla = mla, phi = 0, rseed = 1)
Y <- sim$dat
```

Two anchors per backbone factor, expressed as the two designs `vbpm`
documentation uses throughout.

* **AO (anchor-only)** specifies each anchor's intended cell as `1` and leaves
  that anchor's cells on the *other* backbone columns unspecified (`-1`).
* **AZ (anchor-zero)** additionally fixes those cross-backbone anchor cells to
  `0`.

```{r backbones}
Q0_AO <- matrix(-1L, J, 2)
for (k in 1:2) Q0_AO[which(cluster == k)[1:2], k] <- 1L

Q0_AZ <- Q0_AO
for (k in 1:2) {
  anchors <- which(cluster == k)[1:2]
  Q0_AZ[anchors, ] <- 0L
  Q0_AZ[anchors, k] <- 1L
}

Q0_show <- cbind(Q0_AO[1:8, ], Q0_AZ[1:8, ])
colnames(Q0_show) <- c("AO_1", "AO_2", "AZ_1", "AZ_2")
Q0_show
```

AZ is the stronger claim: it asserts that an anchor of cluster 1 has no loading
on backbone factor 2. AO only asserts where the anchor *does* load. Both
backbones cover two of the three clusters, so the third has to be found by the
sweep.

## Step 1: ordinary oblique sweeps

The calls below use each function's documented defaults; no call in this
vignette passes `v0` or `max_it`. Both `pefa()` and `vbfa()` walk the same
four-step spike ladder `v0 = c(.01, .005, .002, .001)`. What differs is
`max_it`, and it is a cap **per ladder stage**, not per fit or per candidate:
10000 per stage in `pefa()` against 5000 per stage in a direct `vbfa()`, so
the worst-case totals over the four stages are 40000 and 20000. These are two
independent documented defaults, not one shared ceiling.

```{r step1}
AO <- pefa(Q0_AO, Y, Kmin = 2, Kmax = 5, verbose = FALSE)
AZ <- pefa(Q0_AZ, Y, Kmin = 2, Kmax = 5, verbose = FALSE)

AO$sweep[, c("K", "ELBO", "AIC", "BIC", "SRMR", "CFI")]
AZ$sweep[, c("K", "ELBO", "AIC", "BIC", "SRMR", "CFI")]
```

```{r step1-transitions}
AO$transitions[, c("K_from", "K_to", "ELBO_gain_pct", "BIC_gain_pct",
                   "phi_min", "rmsd_max", "unmatched_ssl", "collision")]
AZ$transitions[, c("K_from", "K_to", "ELBO_gain_pct", "BIC_gain_pct",
                   "phi_min", "rmsd_max", "unmatched_ssl", "collision")]
```

Both backbones put their only positive gain on the `2 -> 3` step and both put
`unmatched_ssl` near 2.6 there, which is a full column's worth against the
`ssl()` values below. Both minimize BIC at `K = 3` --- 35212.89 for AO and
35209.52 for AZ --- but only just: the `3 -> 4` BIC step gives back 0.5% (AO)
and 1.0% (AZ) of that first gain, so the criterion is nearly flat between three
and four factors. AIC dissents in both backbones, taking its minimum at `K = 4`
(35015.32 and 35015.94, against 35025.55 and 35022.18 at three).

```{r step1-persistence}
round(AO$persistence$phi, 3)
round(AZ$persistence$phi, 3)
ssl(AO)[c("3", "4")]
ssl(AZ)[c("3", "4")]
```

Persistence separates the two backbones a little more than the criteria do.
AO's three-factor columns survive into `K = 4` and `K = 5` at 0.989; AZ's
survive at 0.985, and AZ's `2 -> 4` and `2 -> 5` cells read 0.905 against AO's
0.925. In both backbones the fourth column is small --- `ssl()` puts it at 0.37
(AO) and 0.45 (AZ) against three columns above 2.4 --- so the four-factor
candidate is the three-factor one plus a thin extra column, not a different
structure.

```{r step1-reading}
## Analysis-owned profiles. ELBO and BIC each use the full fitted window,
## a 20% gain cutoff, sustain one, a post-peak strict crossing, and no fallback.
profile_specs <- data.frame(
  profile = c(".85/r1", ".80/r2", ".70/r3"),
  phi_cut = c(.85, .80, .70), r = 1:3,
  practical_default = c(FALSE, TRUE, FALSE)
)
profile_sources <- 2:3

read_sweep <- function(x, phi_cut, r, sources) {
  K <- as.integer(x$sweep$K)
  converged <- x$sweep$converged

  read_count <- function(score) {
    usable <- length(K) >= 2L && all(is.finite(K)) && all(diff(K) == 1L) &&
      length(score) == length(K) && all(is.finite(score)) &&
      is.logical(converged) && all(converged %in% TRUE)
    if (!usable) return(list(usable = FALSE, Khat = NA_integer_))
    gain <- diff(score)
    gain_max <- max(gain)
    if (!(gain_max > 0)) return(list(usable = TRUE, Khat = NA_integer_))
    peak <- max(which(gain == gain_max))
    hit <- which(seq_along(gain) > peak & gain < .20 * gain_max)
    list(usable = TRUE,
         Khat = if (length(hit)) K[hit[1L]] else NA_integer_)
  }

  ELBO <- read_count(x$sweep$ELBO)
  BIC <- read_count(-x$sweep$BIC)
  C20 <- sort(unique(c(ELBO$Khat, BIC$Khat)[
    is.finite(c(ELBO$Khat, BIC$Khat))]))

  phi <- x$persistence$phi
  collision <- x$persistence$collision
  conv <- setNames(converged, K)
  state <- vapply(sources, function(k) {
    targets <- k + seq_len(r)
    edge <- vapply(targets, function(target) {
      from <- as.character(k); to <- as.character(target)
      if (!from %in% rownames(phi) || !to %in% colnames(phi) ||
          !isTRUE(conv[[from]]) || !isTRUE(conv[[to]]) ||
          is.na(collision[from, to])) return(NA_integer_)
      if (isTRUE(collision[from, to])) return(0L)
      value <- phi[from, to]
      if (!is.finite(value)) NA_integer_ else as.integer(value >= phi_cut)
    }, integer(1))
    if (any(edge %in% 0L)) 0L else
      if (length(edge) && all(edge %in% 1L)) 1L else NA_integer_
  }, integer(1))

  persistent <- which(state == 1L)
  if (!length(state)) {
    resolution <- "unresolved"
    Kp <- NA_integer_
  } else if (!length(persistent)) {
    resolution <- if (all(state %in% 0L)) "nonpersistent" else "unresolved"
    Kp <- NA_integer_
  } else {
    top <- max(persistent)
    higher <- if (top < length(state)) state[(top + 1L):length(state)] else 0L
    resolution <- if (all(higher %in% 0L)) "persistent" else "unresolved"
    Kp <- if (resolution == "persistent") sources[top] else NA_integer_
  }

  hit_ELBO <- is.finite(Kp) && is.finite(ELBO$Khat) && ELBO$Khat == Kp
  hit_BIC <- is.finite(Kp) && is.finite(BIC$Khat) && BIC$Khat == Kp
  support <- if (hit_ELBO && hit_BIC) "both" else if (hit_ELBO) "ELBO only" else
    if (hit_BIC) "BIC only" else "none"
  any_count_usable <- ELBO$usable || BIC$usable
  layer <- if (resolution == "persistent" && Kp %in% C20) "L1" else
    if (resolution == "persistent" && any_count_usable) "L2" else
    if (resolution == "nonpersistent") "L3" else
      "unclassified"
  gap <- if (resolution == "persistent" && length(C20))
    min(abs(C20 - Kp)) else NA_real_

  list(Khat_ELBO = ELBO$Khat, Khat_BIC = BIC$Khat, C20 = C20,
       count_usable = c(ELBO = ELBO$usable, BIC = BIC$usable),
       support = support, Kp = Kp, source_state = setNames(state, sources),
       resolution = resolution, layer = layer, g_CP = gap)
}

profile_row <- function(x, backbone, spec) {
  z <- read_sweep(x, spec$phi_cut, spec$r, profile_sources)
  data.frame(backbone = backbone, profile = spec$profile,
             practical_default = spec$practical_default,
             Khat_ELBO = z$Khat_ELBO, Khat_BIC = z$Khat_BIC,
             C20 = if (length(z$C20)) paste0("{", paste(z$C20,
               collapse = ","), "}") else "{}",
             support = z$support, Kp = z$Kp, layer = z$layer,
             g_CP = z$g_CP, row.names = NULL)
}

sweeps <- list(AO = AO, AZ = AZ)
profiles <- do.call(rbind, lapply(names(sweeps), function(backbone)
  do.call(rbind, lapply(seq_len(nrow(profile_specs)), function(i)
    profile_row(sweeps[[backbone]], backbone, profile_specs[i, ])))))
profiles
```

ELBO and BIC are co-primary count readings here: each uses the 20% rule with
sustain one over all of `K = 2:5`, and finite readings form $C_{20}$. The
`support` column says descriptively whether both, only ELBO, only BIC, or
neither corroborates $K_p$; it does not rank the paths. The three structural
profiles are analysis-owned and scenario-bound. `.80/r2` was retained as the
practical default in one simulation study, not made a **vbpm** default.
The profiles share one declared source set, so the comparison is like-for-like;
here it is `2:3` because this demonstration window is short. Because this short window ends at five, `.70/r3` lacks the
`3 -> 6` edge and is therefore unclassified rather than silently searching a
smaller source universe. Fitting through `K = 6` would resolve that row.

Within either backbone, agreement among the simultaneously displayed resolved
profiles is **within-family robustness** to these declared cutoff/depth choices,
not evidence of population truth, validation, or a package-selected count.
Disagreement identifies sensitivity within the family; it does not say which
profile is true. The `.70/r3` row is also a concrete short-window unresolved
example, not structural nonpersistence.

L1 means that selected $K_p$ belongs to finite $C_{20}$. L2 means that selected
$K_p$ does not belong to $C_{20}$ while at least one count path is usable. L3
requires all declared sources to resolve nonpersistent regardless of count-path
usability; otherwise the row is unclassified. The displayed
$g_{CP}=\min_{\widehat K\in C_{20}}|\widehat K-K_p|$ is corroboration only,
never a decision criterion, and is undefined without a finite count. Here
$g_{CP}=0$ algebraically restates L1 membership; only a positive L2 gap has
descriptive magnitude. A horizon is always direct: $r=2$ asks separately for
`K -> K + 1` and `K -> K + 2`,
never for a chain through the intervening fit. A known collision vetoes its
edge, and the reader selects only the highest persistent source whose higher
declared sources have all resolved nonpersistent.

## Step 2: a comparison derived from Step 1

The two criterion paths supply descriptive finite set $C_{20}$. A persistence
profile separately supplies $K_p$ and the layer. Only L1 or L2 delivers a
candidate to Step 2, and then $K_{\mathrm{step2}}=K_p$. Thus
$K_{\mathrm{step2}}$ is literally $K_p$ conditional on L1/L2 delivery; it is
neither another selector nor the true population count. The group Q below is
derived from the AO candidate at that delivered count. All of this is ordinary
analysis-side code: `pefa()` neither performs it nor stores it.

Let `K0 = ncol(Q0)`. The **specification depth** is the number of markers in
each backbone column. Here both columns have depth two. If the depths differed,
the analyst would need to state a rule of their own or supply a substantive
`J x K_step2` group-loading Q.

The Step-2 design keeps `Q0` unchanged and appends one all-`-1` column for each
new factor. In every new column, the `d` largest absolute loadings become
markers (`1`). Exact ties go to the lower item index. There is no PIP or loading
cutoff, and the same item may mark more than one new factor.

```{r step2-design}
default_AO <- profiles[profiles$backbone == "AO" &
                       profiles$profile == ".80/r2", , drop = FALSE]
if (nrow(default_AO) != 1L || !default_AO$layer %in% c("L1", "L2"))
  stop("The declared profile did not deliver a Step-2 proposal.")
K_step2 <- as.integer(default_AO$Kp)
K0      <- ncol(AO$Q0)
Lam     <- AO$loadings[[as.character(K_step2)]]
d     <- unique(colSums(AO$Q0 == 1L))
if (length(d) != 1L || d[1] < 1L) {
  stop("Choose one marker depth before building Q2.")
}

Q2 <- cbind(AO$Q0, matrix(-1L, J, K_step2 - K0))
colnames(Q2) <- paste0("F", seq_len(K_step2))
if (K_step2 > K0) {
  for (k in (K0 + 1L):K_step2) {
    top <- head(order(-abs(Lam[, k]), seq_len(J)), d)
    Q2[top, k] <- 1L
  }
}

c(K_step2 = K_step2, K0 = K0, depth = d)
Q2[rowSums(Q2 == 1L) > 0L, , drop = FALSE]
```

The added column inherits the backbone depth of two and uses items 13 and 18,
both from the cluster the two-column backbone did not cover.

The same group Q now goes to both models. The bifactor general column is added
inside `vbfa()` and never participates in marker selection.

```{r step2-fit}
f_obl <- vbfa(Y, Q2)
f_bif <- vbfa(Y, Q2, bifactor = TRUE)

idx <- c("RMSEA", "SRMR", "CFI", "TLI", "AIC", "BIC")
round(rbind(
  oblique  = fit_stats(f_obl)[idx],
  bifactor = fit_stats(f_bif)[idx]
), 3)
```

The bifactor arm is ahead on every index, which is the right answer for this
population, but by very unequal margins: 71.2 AIC points and only 5.6 BIC
points. It is an operational comparison of two fitted models on one design, not
a test against an unrestricted alternative, and it is conditional on the count
and the markers that Step 1 delivered. A thin margin, or a margin that flips
when the Step-1 count changes, should be reported rather than resolved.

An analyst may instead supply a substantive `J x K_step2` group-loading Q of their
own, never one containing the bifactor general column. That is a different
analysis design: it should be labelled as such, recorded, validated for codes
and dimnames, and used identically in both arms. It is not a silent default and
it is not a `pefa()` decision.

`special_effects()` then describes the retained bifactor fit: within-cluster
Schmid-Leiman ratios and their CVs, the higher-order parameters those ratios
imply, and group effect sizes. Print the object rather than pulling out
columns --- the print method carries the reading rule the numbers need.

```{r step2-effects}
special_effects(f_bif)
```

The item-weighted mean CV is 0.297, roughly three times the ~.1 convention,
and that is the expected answer here: the population drew each item's general
and group loading independently, so the within-cluster ratios cannot be
proportional by construction. These data are genuinely bifactor, not a
higher-order model in Schmid-Leiman clothing. The `gamma` column is still
printed, but for this fit it is the deliberately lossy summary the print
method warns about and must not be reported as a set of second-order
loadings.

## The same window as a bifactor sweep

`pefa(..., bifactor = TRUE)` sweeps the number of group factors. The backbone
supplies group columns only, the general column is added internally, and the
comparator strips it before matching, so the correspondence and the seven pair
facts have exactly the definitions used above.

```{r bifactor-sweep}
BI <- pefa(Q0_AO, Y, Kmin = 2, Kmax = 5, bifactor = TRUE, verbose = FALSE)
BI$sweep[, c("K", "ELBO", "AIC", "BIC", "SRMR", "CFI")]
BI$transitions[, c("K_from", "K_to", "ELBO_gain_pct", "BIC_gain_pct",
                   "phi_min", "rmsd_max", "unmatched_ssl", "collision")]
round(BI$persistence$phi, 3)
round(BI$persistence$rmsd, 3)
ssl(BI)
```

In this example the bifactor sweep is much less sharp than the ordinary one.

* `ELBO_gain_pct` is `NA` for every row. ELBO is highest at the window's lower
  edge and falls at every step, so that path contains no positive gain and its
  whole percentage column is `NA` by construction. The declared reader above
  treats the complete finite raw ELBO path as usable but returns no finite
  `Khat_ELBO`; that convention belongs to the analysis, not the package.
  `BIC_gain_pct` does have a positive step, at `3 -> 4`, and normalizes to it.
* `phi_min` is 0.998 or better in every cell that involves `K = 2` and 0.990 at
  `4 -> 5`, then drops to 0.444 and 0.519 in the two cells that start at
  `K = 3`. No `rmsd_max` cell exceeds 0.122, so no column moved far in absolute
  terms; those two low congruences are the weakest-link statistic scoring a
  group column that has been emptied.
* `ssl()` says why. The general column holds 4.29-4.37 and the two backbone
  group columns hold about 2.09-2.25 at every candidate, while the third group
  column reaches only 0.27 at `K = 3` and collapses to the order of 1e-9 to
  1e-7 at `K = 4` and `K = 5`. With the general factor already carrying the
  shared variance, the third cluster's specific variance never accumulates into
  a group column here. The four- and five-factor candidates are the two-group
  solution again --- `$persistence$phi` reads 1.000 and 0.999 against `K = 2`
  --- which is also why the four-factor candidate's two empty columns pick the
  same target and raise `collision` at `4 -> 5`.

Read that as a property of this example --- this population, this two-column
backbone, this window, this sample size --- and not as a general ranking. A
bifactor sweep measures the same quantities as an ordinary one, and whether it
separates candidates as well depends on how much group-specific variance is
left once a general factor is in the model. A study that needs both readings
should run both on the same data, backbone, and window rather than assume they
must agree.

## References

* Schmid, J., & Leiman, J. M. (1957). The development of hierarchical factor
  solutions. *Psychometrika*, 22(1), 53-61.
* Yung, Y.-F., Thissen, D., & McLeod, L. D. (1999). On the relationship
  between the higher-order factor model and the hierarchical factor model.
  *Psychometrika*, 64(2), 113-128.
* Zhang, Y., & Chen, J. (2024). Accommodating and extending various models for
  special effects within the generalized partially confirmatory factor
  analysis framework. *Applied Psychological Measurement*, 48(4-5), 208-229.
  \doi{10.1177/01466216241261704}
* 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.
