Bifactor Structures and PEFA Sweeps

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.

A bifactor population and two backbones

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

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.

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
#>      AO_1 AO_2 AZ_1 AZ_2
#> [1,]    1   -1    1    0
#> [2,]    1   -1    1    0
#> [3,]   -1   -1   -1   -1
#> [4,]   -1   -1   -1   -1
#> [5,]   -1   -1   -1   -1
#> [6,]   -1   -1   -1   -1
#> [7,]   -1    1    0    1
#> [8,]   -1    1    0    1

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.

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")]
#>   K      ELBO      AIC      BIC       SRMR       CFI
#> 1 2 -18054.14 35786.49 35992.56 0.07899512 0.8578269
#> 2 3 -17668.64 35025.55 35212.89 0.04212484 0.9853987
#> 3 4 -17687.93 35015.32 35216.70 0.04121520 0.9876296
#> 4 5 -17719.45 35026.34 35246.45 0.04213192 0.9864462
AZ$sweep[, c("K", "ELBO", "AIC", "BIC", "SRMR", "CFI")]
#>   K      ELBO      AIC      BIC       SRMR       CFI
#> 1 2 -18062.99 35798.46 36009.21 0.07998786 0.8559784
#> 2 3 -17669.53 35022.18 35209.52 0.04179432 0.9859667
#> 3 4 -17684.05 35015.94 35217.32 0.04176480 0.9875248
#> 4 5 -17715.08 35028.18 35248.30 0.04283681 0.9861359
AO$transitions[, c("K_from", "K_to", "ELBO_gain_pct", "BIC_gain_pct",
                   "phi_min", "rmsd_max", "unmatched_ssl", "collision")]
#>   K_from K_to ELBO_gain_pct BIC_gain_pct   phi_min    rmsd_max unmatched_ssl
#> 1      2    3    100.000000  100.0000000 0.9268021 0.173931747   2.618378591
#> 2      3    4     -5.005149   -0.4890918 0.9889516 0.056553150   0.366505391
#> 3      4    5     -8.174384   -3.8163641 0.9998816 0.006025838   0.001363168
#>   collision
#> 1     FALSE
#> 2     FALSE
#> 3     FALSE
AZ$transitions[, c("K_from", "K_to", "ELBO_gain_pct", "BIC_gain_pct",
                   "phi_min", "rmsd_max", "unmatched_ssl", "collision")]
#>   K_from K_to ELBO_gain_pct BIC_gain_pct   phi_min   rmsd_max unmatched_ssl
#> 1      2    3    100.000000  100.0000000 0.9267992 0.16111415   2.653897041
#> 2      3    4     -3.691467   -0.9759852 0.9853839 0.06545695   0.454986448
#> 3      4    5     -7.884652   -3.8733361 0.9999756 0.00507446   0.002154502
#>   collision
#> 1     FALSE
#> 2     FALSE
#> 3     FALSE

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

round(AO$persistence$phi, 3)
#>    2     3     4     5
#> 2 NA 0.927 0.925 0.925
#> 3 NA    NA 0.989 0.989
#> 4 NA    NA    NA 1.000
#> 5 NA    NA    NA    NA
round(AZ$persistence$phi, 3)
#>    2     3     4     5
#> 2 NA 0.927 0.905 0.905
#> 3 NA    NA 0.985 0.985
#> 4 NA    NA    NA 1.000
#> 5 NA    NA    NA    NA
ssl(AO)[c("3", "4")]
#> $`3`
#>       F1       F2       F3 
#> 3.306630 2.647015 2.618379 
#> 
#> $`4`
#>        F1        F2        F3        F4 
#> 3.3056530 2.5181604 2.5432713 0.3665054
ssl(AZ)[c("3", "4")]
#> $`3`
#>       F1       F2       F3 
#> 3.353168 2.663625 2.653897 
#> 
#> $`4`
#>        F1        F2        F3        F4 
#> 3.3337463 2.4235747 2.5432108 0.4549864

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.

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

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.

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)
#> K_step2      K0   depth 
#>       3       2       2
Q2[rowSums(Q2 == 1L) > 0L, , drop = FALSE]
#>         F1 F2 F3
#> item_1   1 -1 -1
#> item_2   1 -1 -1
#> item_7  -1  1 -1
#> item_8  -1  1 -1
#> item_13 -1 -1  1
#> item_18 -1 -1  1

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.

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)
#>          RMSEA  SRMR   CFI   TLI      AIC      BIC
#> oblique  0.029 0.042 0.985 0.983 35025.40 35212.74
#> bifactor 0.004 0.020 1.000 1.000 34954.23 35207.14

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.

special_effects(f_bif)
#> Special effects of a bifactor fit (3 group factors; tau = 0.5)
#>  group n_items mean_ratio sd_ratio    cv gamma     D
#>     F1       6      1.007    0.208 0.206 0.705 0.273
#>     F2       6      1.367    0.552 0.404 0.591 0.289
#>     F3       6      0.959    0.270 0.282 0.722 0.195
#> Item-weighted mean CV of the group/general loading ratio: 0.297
#>   (near-zero CV = Schmid-Leiman proportional = the fit re-expresses a
#>   higher-order/testlet model, and gamma/lambda are faithful; values
#>   under ~.1 are conventionally read that way -- suggestive, not a test.
#>   With clearly non-constant ratios the higher-order parameters are a
#>   lossy summary of a genuinely richer bifactor.)

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.

BI <- pefa(Q0_AO, Y, Kmin = 2, Kmax = 5, bifactor = TRUE, verbose = FALSE)
BI$sweep[, c("K", "ELBO", "AIC", "BIC", "SRMR", "CFI")]
#>   K      ELBO      AIC      BIC       SRMR       CFI
#> 1 2 -17673.10 35014.02 35243.51 0.04362376 0.9888589
#> 2 3 -17686.17 35040.73 35270.21 0.04744536 0.9843579
#> 3 4 -17699.60 35013.95 35243.43 0.04358478 0.9888711
#> 4 5 -17711.41 35038.14 35262.94 0.04748237 0.9846260
BI$transitions[, c("K_from", "K_to", "ELBO_gain_pct", "BIC_gain_pct",
                   "phi_min", "rmsd_max", "unmatched_ssl", "collision")]
#>   K_from K_to ELBO_gain_pct BIC_gain_pct   phi_min   rmsd_max unmatched_ssl
#> 1      2    3            NA    -99.72806 0.9984761 0.02023850  2.676671e-01
#> 2      3    4            NA    100.00000 0.4436801 0.12193700  3.530298e-09
#> 3      4    5            NA    -72.83562 0.9904581 0.01713767  5.085772e-07
#>   collision
#> 1     FALSE
#> 2     FALSE
#> 3      TRUE
round(BI$persistence$phi, 3)
#>    2     3     4     5
#> 2 NA 0.998 1.000 0.999
#> 3 NA    NA 0.444 0.519
#> 4 NA    NA    NA 0.990
#> 5 NA    NA    NA    NA
round(BI$persistence$rmsd, 3)
#>    2    3     4     5
#> 2 NA 0.02 0.000 0.017
#> 3 NA   NA 0.122 0.122
#> 4 NA   NA    NA 0.017
#> 5 NA   NA    NA    NA
ssl(BI)
#> $`2`
#>        G       F1       F2 
#> 4.294704 2.247513 2.086900 
#> 
#> $`3`
#>         G        F1        F2        F3 
#> 4.3649111 2.1724048 2.0881266 0.2676671 
#> 
#> $`4`
#>            G           F1           F2           F3           F4 
#> 4.294599e+00 2.247053e+00 2.085096e+00 3.530298e-09 4.783997e-09 
#> 
#> $`5`
#>            G           F1           F2           F3           F4           F5 
#> 4.358735e+00 2.175613e+00 2.088907e+00 4.319029e-07 1.538884e-07 5.085772e-07

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

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