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