PEFA: Factor-Count Evidence from a Fixed Window

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.

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

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, ]
#>      [,1] [,2]
#> [1,]    1   -1
#> [2,]    1   -1
#> [3,]   -1   -1
#> [4,]   -1   -1
#> [5,]   -1   -1
#> [6,]   -1   -1
#> [7,]   -1    1
#> [8,]   -1    1

p <- pefa(Q0, Y, Kmin = 2, Kmax = 7, verbose = FALSE)
p
#> PEFA sweep (ordinary): K = 2:7, K0 = 2
#>   sweep 6 x 11 | transitions 5 x 11 | persistence 3 x 6 x 6
#>   all candidates converged

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

names(p)
#> [1] "sweep"       "transitions" "persistence" "loadings"    "pips"       
#> [6] "Q0"          "settings"
class(p)
#> [1] "pefa"
names(p$settings)
#> [1] "bifactor"    "general"     "v0"          "max_it"      "convChk"    
#> [6] "tolVal"      "tau"         "rank_adjust" "rank_max_J"

$settings holds the resolved fitting and fit-statistic controls only — no threshold or rule. Candidate matrices are stored in two K-named lists.

names(p$loadings)
#> [1] "2" "3" "4" "5" "6" "7"
names(p$pips)
#> [1] "2" "3" "4" "5" "6" "7"
round(p$loadings[["3"]], 2)
#>            F1    F2    F3
#> item_1   0.69  0.02  0.02
#> item_2   0.74 -0.03 -0.02
#> item_3   0.65  0.03  0.01
#> item_4   0.67 -0.01  0.01
#> item_5   0.74  0.24  0.00
#> item_6   0.67 -0.21 -0.01
#> item_7   0.00  0.70 -0.01
#> item_8  -0.02  0.70  0.00
#> item_9   0.00  0.76 -0.01
#> item_10  0.03  0.69  0.02
#> item_11 -0.03  0.68  0.35
#> item_12  0.02  0.72 -0.37
#> item_13  0.02  0.00  0.68
#> item_14 -0.01  0.02  0.68
#> item_15  0.02 -0.01  0.69
#> item_16  0.00  0.00  0.67
#> item_17  0.30 -0.02  0.72
#> item_18 -0.32  0.01  0.66
round(p$pips[["3"]][1:6, ], 2)
#>        F1   F2   F3
#> item_1  1 0.04 0.04
#> item_2  1 0.04 0.04
#> item_3  1 0.05 0.03
#> item_4  1 0.03 0.03
#> item_5  1 1.00 0.03
#> item_6  1 1.00 0.03

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

p$sweep
#>   K      ELBO      AIC      BIC      RMSEA       SRMR       CFI       TLI  t
#> 1 2 -11253.00 22192.69 22373.83 0.10410846 0.07996783 0.8290836 0.7957015 43
#> 2 3 -10936.65 21519.83 21709.40 0.01664782 0.03509306 0.9956978 0.9947759 45
#> 3 4 -10964.55 21541.92 21739.91 0.02462359 0.04185642 0.9907375 0.9885713 47
#> 4 5 -10993.37 21611.93 21822.56 0.04160641 0.04828957 0.9741948 0.9673703 50
#> 5 6 -11021.89 21619.28 21850.97 0.04296856 0.04687183 0.9736147 0.9651988 55
#> 6 7 -11048.46 21624.82 21877.58 0.04403729 0.05871339 0.9734805 0.9634460 60
#>   iter converged
#> 1  478      TRUE
#> 2  300      TRUE
#> 3 1105      TRUE
#> 4 1991      TRUE
#> 5 1046      TRUE
#> 6  906      TRUE

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

p$transitions
#>   K_from K_to ELBO_gain_pct BIC_gain_pct    phi_min       rmsd    rmsd_max
#> 1      2    3    100.000000   100.000000 0.04111435 0.46946389 0.575890437
#> 2      3    4     -8.820479    -4.592867 0.99563527 0.02549756 0.038952081
#> 3      4    5     -9.107501   -12.438138 0.65585390 0.03035042 0.051782096
#> 4      5    6     -9.017905    -4.276202 0.99793493 0.00293848 0.003389433
#> 5      6    7     -8.398175    -4.004473 0.29384297 0.11313264 0.254971534
#>         ari    pip_rmsd unmatched_ssl collision
#> 1 0.5405405 0.674665357   3.049497001     FALSE
#> 2 1.0000000 0.134071252   0.002982102     FALSE
#> 3 1.0000000 0.115030787   0.006564887     FALSE
#> 4 1.0000000 0.005471327   0.004728663      TRUE
#> 5 0.8495575 0.237019711   0.641606891      TRUE

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.

names(p$persistence)
#> [1] "phi"       "rmsd"      "collision"
round(p$persistence$phi, 3)
#>    2     3     4     5     6     7
#> 2 NA 0.041 0.037 0.036 0.034 0.052
#> 3 NA    NA 0.996 0.990 0.989 0.789
#> 4 NA    NA    NA 0.656 0.651 0.421
#> 5 NA    NA    NA    NA 0.998 0.284
#> 6 NA    NA    NA    NA    NA 0.294
#> 7 NA    NA    NA    NA    NA    NA
round(p$persistence$rmsd, 3)
#>    2     3     4     5     6     7
#> 2 NA 0.576 0.576 0.572 0.572 0.563
#> 3 NA    NA 0.039 0.059 0.062 0.253
#> 4 NA    NA    NA 0.052 0.054 0.254
#> 5 NA    NA    NA    NA 0.003 0.255
#> 6 NA    NA    NA    NA    NA 0.255
#> 7 NA    NA    NA    NA    NA    NA
p$persistence$collision
#>    2     3     4     5     6     7
#> 2 NA FALSE FALSE FALSE FALSE FALSE
#> 3 NA    NA FALSE FALSE FALSE FALSE
#> 4 NA    NA    NA FALSE FALSE FALSE
#> 5 NA    NA    NA    NA  TRUE  TRUE
#> 6 NA    NA    NA    NA    NA  TRUE
#> 7 NA    NA    NA    NA    NA    NA

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.

round(p$loadings[["2"]][c(1, 7, 8, 13, 17, 18), ], 2)
#>            F1   F2
#> item_1   0.66 0.01
#> item_7   0.54 0.02
#> item_8   0.49 0.08
#> item_13  0.01 0.69
#> item_17  0.25 0.70
#> item_18 -0.29 0.67

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.

ssl(p)
#> $`2`
#>       F1       F2 
#> 4.624276 3.090395 
#> 
#> $`3`
#>       F1       F2       F3 
#> 3.082710 3.135258 3.049497 
#> 
#> $`4`
#>          F1          F2          F3          F4 
#> 2.927150934 3.103792975 3.059204940 0.002982102 
#> 
#> $`5`
#>          F1          F2          F3          F4          F5 
#> 3.146855076 3.020348740 3.073076195 0.006659368 0.006564887 
#> 
#> $`6`
#>          F1          F2          F3          F4          F5          F6 
#> 3.151805952 3.014687889 3.080020875 0.004117262 0.004728663 0.004706185 
#> 
#> $`7`
#>          F1          F2          F3          F4          F5          F6 
#> 2.930601005 2.936020782 0.641606891 0.569967670 1.852713426 0.002794206 
#>          F7 
#> 0.002805425

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.

round(p$loadings[["7"]][13:18, 3:5], 2)
#>           F3    F4   F5
#> item_13 0.00 -0.01 0.71
#> item_14 0.01  0.00 0.68
#> item_15 0.00  0.00 0.70
#> item_16 0.01  0.02 0.64
#> item_17 0.70  0.00 0.01
#> item_18 0.00  0.75 0.00

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.

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.

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
#>   profile practical_default Khat_ELBO Khat_BIC C20 support Kp layer g_CP
#> 1  .85/r1             FALSE         3        3 {3}    both  3    L1    0
#> 2  .80/r2              TRUE         3        3 {3}    both  3    L1    0
#> 3  .70/r3             FALSE         3        3 {3}    both  3    L1    0

## 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)
#>   profile declared_source required_endpoint resolution        layer
#> 1  .85/r1               7                 8 unresolved unclassified

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

summary(p)
#> PEFA sweep summary (ordinary): K = 2:7, K0 = 2
#> 
#> Candidates:
#>  K     ELBO     AIC     BIC RMSEA  SRMR   CFI   TLI  t iter converged
#>  2 -11253.0 22192.7 22373.8 0.104 0.080 0.829 0.796 43  478      TRUE
#>  3 -10936.7 21519.8 21709.4 0.017 0.035 0.996 0.995 45  300      TRUE
#>  4 -10964.6 21541.9 21739.9 0.025 0.042 0.991 0.989 47 1105      TRUE
#>  5 -10993.4 21611.9 21822.6 0.042 0.048 0.974 0.967 50 1991      TRUE
#>  6 -11021.9 21619.3 21851.0 0.043 0.047 0.974 0.965 55 1046      TRUE
#>  7 -11048.5 21624.8 21877.6 0.044 0.059 0.973 0.963 60  906      TRUE
#> 
#> SS loadings (colSums(Lam^2)):
#>   K=2: 4.62 3.09
#>   K=3: 3.08 3.14 3.05
#>   K=4: 2.93 3.10 3.06 0.00
#>   K=5: 3.15 3.02 3.07 0.01 0.01
#>   K=6: 3.15 3.01 3.08 0.00 0.00 0.00
#>   K=7: 2.93 2.94 0.64 0.57 1.85 0.00 0.00
#> 
#> Transitions:
#>  K_from K_to ELBO_gain_pct BIC_gain_pct phi_min  rmsd rmsd_max   ari pip_rmsd
#>       2    3       100.000      100.000   0.041 0.469    0.576 0.541    0.675
#>       3    4        -8.820       -4.593   0.996 0.025    0.039 1.000    0.134
#>       4    5        -9.108      -12.438   0.656 0.030    0.052 1.000    0.115
#>       5    6        -9.018       -4.276   0.998 0.003    0.003 1.000    0.005
#>       6    7        -8.398       -4.004   0.294 0.113    0.255 0.850    0.237
#>  unmatched_ssl collision
#>          3.049     FALSE
#>          0.003     FALSE
#>          0.007     FALSE
#>          0.005      TRUE
#>          0.642      TRUE
#> 
#> Persistence (direct K_from x K_to comparisons):
#>   phi (phi_min):
#>   2     3     4     5      6      7
#> 2   0.041 0.037 0.036  0.034  0.052
#> 3         0.996 0.990  0.989  0.789
#> 4               0.656  0.651  0.421
#> 5                     0.998* 0.284*
#> 6                            0.294*
#> 7                                  
#>   rmsd (rmsd_max):
#>   2     3     4     5      6      7
#> 2   0.576 0.576 0.572  0.572  0.563
#> 3         0.039 0.059  0.062  0.253
#> 4               0.052  0.054  0.254
#> 5                     0.003* 0.255*
#> 6                            0.255*
#> 7                                  
#>   * = target reuse (collision)

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.

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