--- title: "PEFA: Factor-Count Evidence from a Fixed Window" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{PEFA: Factor-Count Evidence from a Fixed Window} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>") library(vbpm) ``` `pefa()` fits one partially exploratory factor model for every integer `K` in a fixed consecutive window and reports what it measured: a candidate table, the adjacent criterion gains, and the direct structural comparison of every ordered pair of candidates. It **chooses nothing**: there is no count, cut, threshold, or stopping rule inside the package, and no argument that would add one. This vignette walks one sweep end to end, deliberately an *overextraction* example: the population has three factors and the window runs to seven, so the point is what the three tables look like when the extra columns are surplus. Bifactor sweeps and the step that follows a count are in `vignette("bifactor")`. ## Decision inventory Six things must be decided before a factor count comes out of a sweep. The package decides none of them; each line below names the field that carries the evidence. * **The window.** `Kmin:Kmax` bounds every number under it; a reading at the edge of `$sweep$K` is a reading about the window. * **The criterion.** `$sweep` carries `ELBO`, `AIC`, `BIC`, and four descriptive indices side by side and ranks none of them. * **The count rule and its variant.** `$transitions$ELBO_gain_pct` and `BIC_gain_pct` are the raw material; the cut, the look-ahead depth, and the boundary fallback are yours. * **The collision policy.** `$transitions$collision` and `$persistence$collision` record target reuse; whether reuse disqualifies a candidate is your call. * **The persistence horizon.** A row of `$persistence$phi` says how far a candidate's columns survive; how far they must survive is not in the object. * **The tie-break.** When two criterion paths disagree, nothing arbitrates: `$sweep` reports both. ## One sweep Three factors correlated at 0.5, six items each with main loadings of 0.7, and two cross-loadings of 0.3 per factor with alternating signs --- the `sim_fa()` defaults. The backbone is **anchor-only (AO)**: two anchors per backbone factor get their intended cell fixed to `1`, and every other cell stays `-1` (regularized, otherwise unspecified). ```{r sweep} sim <- sim_fa(N = 500, K = 3, ipf = 6, lam = .7, lac = .3, rseed = 11) Y <- sim$dat groups <- rep(1:3, each = 6) Q0 <- matrix(-1L, ncol(Y), 2) for (k in 1:2) Q0[which(groups == k)[1:2], k] <- 1L Q0[1:8, ] p <- pefa(Q0, Y, Kmin = 2, Kmax = 7, verbose = FALSE) p ``` `persistence 3 x 6 x 6` in that header means three 6 x 6 matrices, one per measure, not a three-dimensional array. The backbone has `K0 = 2` columns, so the window may start no lower than 2. For each candidate `K`, `pefa()` appends `K - K0` fully regularized columns to `Q0` and fits that design with `vbfa()`; the first two group columns keep the backbone's identity by position, the rest are exploratory. The call omits `v0` and `max_it`, so every candidate is fitted at the documented `pefa()` defaults: the four-step spike ladder `c(.01, .005, .002, .001)` and `max_it = 10000`, a cap **per ladder stage** rather than per fit. `vbfa()` defaults to the same ladder and differs only in that cap, at 5000. The object has exactly seven components and class `"pefa"`. ```{r object} names(p) class(p) names(p$settings) ``` `$settings` holds the resolved fitting and fit-statistic controls only --- no threshold or rule. Candidate matrices are stored in two K-named lists. ```{r matrices} names(p$loadings) names(p$pips) round(p$loadings[["3"]], 2) round(p$pips[["3"]][1:6, ], 2) ``` A PIP is the posterior probability that a regularized cell is selected. Fixed backbone cells are not regularized, so their PIPs are structural rather than informative. ## The candidate table ```{r sweep-table} p$sweep ``` One row per `K`, always the same eleven columns. `AIC` and `BIC` are the hard-selection values at the resolved `tau` and `t` is their effective parameter count. `RMSEA`, `SRMR`, `CFI`, and `TLI` are descriptive and may be `NA` when their own definition does not apply --- `RMSEA` and `TLI` need positive degrees of freedom. `converged` reflects the last ladder stage only: a candidate whose last stage ran to `max_it` is retained with `converged = FALSE` and one aggregate warning, and `iter` totals the iterations over all four stages. Every candidate converged here. ELBO peaks at `K = 3` (-10936.65) and falls afterwards. BIC is minimized at the same candidate, 21709.40, and rises at every later one. AIC is minimized there as well, at 21519.83, and also rises monotonically --- at the default spike ladder the three criteria agree on this data set. The descriptive indices point the same way rather than adding discrimination: RMSEA (0.017) and SRMR (0.035) are lowest at `K = 3`, and CFI and TLI peak there at 0.996 and 0.995. ## Adjacent transitions ```{r transitions} p$transitions ``` One row per adjacent pair. `ELBO_gain_pct` and `BIC_gain_pct` are oriented so that positive favors `K_to`: the raw gains are `ELBO_to - ELBO_from` and `BIC_from - BIC_to`, and each is divided by the largest positive gain on its own path and multiplied by 100. Here the largest positive gain on both paths is the `2 -> 3` step, so that row is 100 and every later row is negative. If a path has no positive gain at all, its whole percentage column is `NA`. The remaining seven columns describe structural change. Backbone columns are paired by position. Each exploratory source column is paired with the closest target after allowing a whole-column sign flip. See `?pefa` for the exact distance and tie rules. The assignment is independent, not one-to-one, so two sources may pick the same target. That reuse is a **collision**, recorded rather than repaired; the `5 -> 6` and `6 -> 7` rows carry it here. ## Persistence, read directly There is no accessor function. The three upper-triangular matrices over the fitted `K` values are read with `$`: `phi` (holding `phi_min`), `rmsd` (holding **`rmsd_max`**, the weakest-column value, not the pooled `rmsd` of `$transitions`), and the logical collision mask. ```{r persistence} names(p$persistence) round(p$persistence$phi, 3) round(p$persistence$rmsd, 3) p$persistence$collision ``` Every cell above the diagonal is a **direct** comparison of its two endpoint loading matrices, not a product of the adjacent steps between them: `3 -> 7` compares the three- and seven-factor solutions themselves. The diagonal and lower triangle are typed `NA`. Read the `K = 3` row across: 0.996 against `K = 4`, 0.990 against `K = 5`, 0.989 against `K = 6`, and then 0.789 against `K = 7`. Because `phi_min` is a weakest-link statistic, 0.989 says the *worst*-matched of the three columns still has a near-copy three candidates on, with the weakest-column RMSD no worse than 0.062. That is what a high persistence row buys: one number retires "the structure moved when I added a column" for the whole range out to `K = 6`. Only the seven-factor cell breaks the row, at RMSD 0.253. The `K = 2` row never exceeds 0.052 anywhere, and the two-factor loading matrix says why. ```{r k2} round(p$loadings[["2"]][c(1, 7, 8, 13, 17, 18), ], 2) ``` Backbone columns are compared by *position*, and at `K = 2` the second backbone column is not holding what it holds from `K = 3` on. Items 13-18 have taken it over --- the three shown load 0.69, 0.70, and 0.67 on it --- while its own anchors, items 7 and 8, sit at 0.02 and 0.08, absorbed into the first column alongside cluster 1. Its congruence with the second column of the three-factor solution is 0.041, and that one weak pair is the whole of `phi_min` for every `2 -> K` cell. Specifying an anchor cell leaves it free to be estimated; it does not force the column to be about that cluster when there are too few columns to go round. The collision mask is `TRUE` at `5 -> 6`, `5 -> 7`, and `6 -> 7`, where two or more source columns selected the same target. Those cells sit at `phi_min` 0.998, 0.284, and 0.294 --- a first sign that a collision is not by itself bad news about the structure. ## Column strength `ssl()` is the descriptive companion to those numbers: `colSums(Lam^2)` for every stored candidate matrix. It screens nothing and enters no correspondence rule. ```{r ssl} ssl(p) ``` This is what a *low* `phi_min` needs to be read against. `phi_min` is a **weakest-link** statistic over *every* source column, including one the fit has driven to near zero. Such a column has no shape: it matches whichever target has the smallest norm, contributes an arbitrary congruence, and, because two of them pick the same target, raises `collision`. The `4 -> 5` row is exactly that case --- `phi_min` 0.656 from a source column of squared length 0.003 matching a target of 0.007 --- while its pooled `rmsd` is 0.030 and its `ari` is 1: not one item changed its dominant column. The `5 -> 6` and `6 -> 7` collisions come from the same near-empty columns. `unmatched_ssl` is on the same scale: the largest sum of squared loadings among the exploratory target columns that no exploratory source selected, the biggest column the smaller candidate did *not* account for. Across the adjacent rows it reads 3.049, 0.003, 0.007, 0.005, 0.642. At `2 -> 3` the surplus column carries 3.05, indistinguishable from the three genuine columns of the three-factor solution (3.08, 3.14, 3.05): the two-factor candidate really was missing a factor. At `3 -> 4`, `4 -> 5`, and `5 -> 6` it is under 0.01: each of those candidates adds a column with nothing in it and leaves the three real columns alone. `6 -> 7` breaks that pattern, and `ssl(p)` shows why. The seven-factor candidate does append two empty columns --- `F6` and `F7`, both at 0.003 --- but neither of them is the *surplus* one. `unmatched_ssl` here is 0.642, which is exactly `ssl(p)[["7"]]["F3"]`, and the third cluster has been broken into pieces of 1.853, 0.642, and 0.570. ```{r k7} round(p$loadings[["7"]][13:18, 3:5], 2) ``` Items 17 and 18 have each been given a column of their own --- a genuinely different structure, and what drives the `3 -> 7` congruence down to 0.789 after three cells at 0.989 or better. The `6 -> 7` value of 0.294 is a separate matter: like `4 -> 5`, it is two empty columns being compared. ## What the three readings say here The criterion path, the persistence row, and `unmatched_ssl` against `ssl()` point the same way in this example. * ELBO, AIC, and BIC all turn at `K = 3`, and every later adjacent gain is negative on both percentage paths. * The three-factor columns persist into `K = 4`, `5`, and `6` at congruence 0.989 or better: the extra columns are added beside that structure. * The column added at `K = 3` carries 3.05, a full factor's worth; those added at `K = 4`, `5`, and `6` carry under 0.01 each. The population had three factors, so the reading is right. **It is still the reader's inference, not a package verdict.** The sentence "the structure is three-dimensional and the later columns are surplus" was composed here, in the vignette, out of six candidate rows, five transition rows, and three triangles. Two of those numbers point the other way if read alone: `phi_min` is 0.656 at `4 -> 5` and 0.294 at `6 -> 7`, structural instability until `ssl()` shows that both belong to columns with nothing in them. ## Analysis-side dual-count and persistence profiles The following compact reader illustrates one declared analysis. ELBO and BIC are co-primary descriptive count paths: each uses a 20% gain cutoff, sustain one, a strict scan after the last maximum gain, the **full fitted window**, and no boundary fallback. Their finite readings form $C_{20}=\{\widehat K_{ELBO},\widehat K_{BIC}\}$ after absent values are removed. Neither path breaks a tie or outranks the other. The same reader applies three scenario-bound minimum-congruence profiles: `.85/r1`, `.80/r2`, and `.70/r3`, with `.80/r2` labelled the practical default for the motivating study. Each source is compared directly with every endpoint through `K + r`; no adjacent chain is substituted. A known collision vetoes that edge. The common declared source set is `2:4`, chosen so the deepest profile has all required endpoints in the fitted window. $K_p$ is the highest persistent source only after every higher declared source has resolved nonpersistent. ```{r reading} 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:4 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, spec) { z <- read_sweep(x, spec$phi_cut, spec$r, profile_sources) data.frame(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) } profile_results <- do.call(rbind, lapply(seq_len(nrow(profile_specs)), function(i) profile_row(p, profile_specs[i, ]))) profile_results ## A declared source at the upper edge cannot persist without its endpoint. short_window <- read_sweep(p, phi_cut = .85, r = 1L, sources = 7L) data.frame(profile = ".85/r1", declared_source = 7L, required_endpoint = 8L, resolution = short_window$resolution, layer = short_window$layer) ``` Here both gain paths read three, so $C_{20}=\{3\}$. The three structural rows also select three and receive L1 with descriptive support from both criteria. The support label records whether one or both count paths corroborate $K_p$; it does not rank the paths or select the structure. Agreement across all three declared profiles is **within-family robustness** to these particular cutoff/depth combinations. It is not evidence that three is the population truth, a validation result, or a package decision. Conversely, profile disagreement would describe sensitivity within this family rather than identify which row is true. The final printed row is the short-window counterexample. Source 7 under `.85/r1` needs the direct `7 -> 8` endpoint, but this fit stops at 7; its persistence resolution and layer are therefore `unresolved` and `unclassified`, not L3 and not a reason to shrink the declared source set. L1 means that a selected $K_p$ belongs to finite $C_{20}$. L2 means that $K_p$ is selected but does not belong to $C_{20}$ while at least one count path is usable. L3 requires every declared source to resolve nonpersistent regardless of count-path usability. All remaining cases are unclassified. The displayed $g_{CP}=\min_{\widehat K\in C_{20}}|\widehat K-K_p|$ is corroboration only, never a persistence or delivery criterion, and is undefined when there is no finite count reading. Algebraically, $g_{CP}=0$ merely restates the L1 membership condition $K_p\in C_{20}$; only positive gap magnitudes add a descriptive measure of how far an L2 count disagreement lies. If an analysis continues to Step 2, only an L1/L2 delivery is eligible and then $K_{\mathrm{step2}}=K_p$. Thus the Step-2 count is $K_p$ conditional on delivery, not another count selector; `vignette("bifactor")` shows the design construction. These definitions, including the three profiles and the practical-default label, belong to this analysis rather than **vbpm**. A different declared source set, horizon, cutoff, collision policy, or window can change the result, and an unlabelled reading copied from a vignette can quietly become a de facto default. Gain fractions are relative to the largest positive gain in the full declared window. Widening that window can change the denominator and the reading; **vbpm** 0.9.1 has no reassembly API, so a wider window requires another complete `pefa()` call. ## Display methods ```{r summary} summary(p) ``` `summary()` returns exactly `window`, `sweep`, `transitions`, `persistence`, `ssl`, `settings`, and `nonconverged_K`. Its print method marks a persistence value with `*` wherever the parallel collision cell is `TRUE`; that is display formatting only, and the stored matrices are untouched. The plot methods are descriptive trajectories, with no threshold line and no marked count. ```{r plots, fig.width=7, fig.height=4.5} plot(p, type = "objective", criterion = "ELBO") plot(p, type = "gain", criterion = "BIC", pct = TRUE) plot(p, type = "fit") ``` The ELBO path turns visibly at `K = 3`, and both fit panels turn with it. The gain view is indexed by `K_to`, so its point at `x = 3` is the `2 -> 3` step, not the three-factor candidate; `pct = TRUE` reads the stored percentage column, and the default `pct = FALSE` recomputes the raw oriented gains from `$sweep`. ## Scope What `pefa()` returns is measurement. The named profiles above illustrate analysis-owned readings; they do not become package selectors or defaults, and there is no automatic window extension. A failed or malformed candidate aborts the call with an error naming its `K`; a nonconverged but well-formed one is retained and reported. The reading above is a property of this simulated example --- 500 observations, main loadings of 0.7, three clean clusters, an anchored backbone --- and not evidence about how these tables behave at weaker signal. ## References * Chen, J. (2023). Fully and partially exploratory factor analysis with bi-level Bayesian regularization. *Behavior Research Methods*, 55(4), 2125-2142. \doi{10.3758/s13428-022-01884-7} * 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.