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 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$datTwo anchors per backbone factor, expressed as the two designs
vbpm documentation uses throughout.
1 and leaves that anchor’s cells on the
other backbone columns unspecified (-1).0.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 1AZ 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.
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.9861359AO$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 FALSEBoth 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.4549864Persistence 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 NAELBO 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.
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 1The 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.14The 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.
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-07In 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.