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