Bass-ackwards hierarchical structural analysis

Extracts factor or component solutions at levels 1 through k. A factor (or component) is a summary variable that stands in for a group of items that move together. It then characterises the hierarchy by computing the correlations between the factor scores of different levels. A factor score is each person's estimated standing on a factor. Each edge is computed from the scoring weights and the correlation matrix stored as x$r, and it is exact for that matrix. For engine = "esem", lavaan fits a different matrix in two settings. With cor = "spearman" and missing = "pairwise" or "listwise", x$r is a Spearman matrix, but lavaan fits Pearson covariances of the raw data. With cor = "pearson" and missing = "pairwise" on data with missing values, x$r holds pairwise correlations. There, the ML, MLR, and ULSMV estimators fit the complete rows, and WLSMV fits the pairwise covariances. The correlations of those covariances differ slightly from x$r, because each item's variance uses all of its observed rows. The "hierarchy" is descriptive: edges are score correlations, not a fitted higher-order SEM.

Usage

ackwards(
  data,
  k_max,
  engine = "pca",
  cor = "pearson",
  fm = "minres",
  estimator = NULL,
  missing = "pairwise",
  n_obs = NULL,
  align_signs = TRUE,
  keep_scores = FALSE,
  keep_fits = FALSE,
  seed = NULL,
  pairs = "adjacent",
  cut_show = 0.3,
  correct = 0.5,
  rotation = "varimax",
  ...
)

Arguments

data

A data frame or numeric matrix of observed variables (items in columns, observations in rows). Alternatively, a pre-computed correlation matrix may be supplied (a square, symmetric, numeric matrix with unit diagonal). When a correlation matrix is supplied, engine must be "pca" or "efa", because ESEM (exploratory structural equation modeling) needs raw data. The missing and cor arguments are then ignored. See the Correlation-matrix input section below.

k_max

Maximum number of factors or components to extract. It is required, so use suggest_k() if uncertain. It sets the depth of the hierarchy: levels 1 through k_max are all extracted and retained. (The k_max of suggest_k() means something related but different. It is the largest number of factors that function evaluates when recommending a depth, not a depth itself. See that function's docs.)

engine

Extraction engine: "pca" (default), "efa", or "esem". "esem" uses lavaan::efa() with rotation-aware SEs and per-level fit indices. ESEM is exploratory structural equation modeling, and it is recommended for the clinical or HiTOP workflow (Kim & Eaton, 2015, Forbush et al., 2024). It requires lavaan >= 0.6-13.

cor

Correlation basis: "pearson" (default), "spearman", or "polychoric". For PCA and EFA, "polychoric" computes a polychoric correlation matrix via psych::polychoric() (requires psych). Polychoric correlations estimate the correlation between the continuous traits assumed to underlie ordered responses. For ESEM, it triggers WLSMV estimation via lavaan. A cli warning is emitted when ordinal-looking columns are detected and cor is not "polychoric". Ordinal items have a few ordered categories, such as a 1 to 5 rating.

fm

Factor extraction method passed to psych::fa(), and used only when engine = "efa". One of "minres" (the default, an ordinary least squares fit that converges reliably), "ml" (maximum likelihood, gives chi-square fit but converges less reliably at deep levels), or "pa" (principal axis). Ignored for engine = "pca".

estimator

Estimation method for the ESEM engine. NULL (default) auto-selects: "WLSMV" when cor = "polychoric", "ML" otherwise. Pass explicitly to override: "ULSMV" (unweighted WLS), "MLR" (maximum likelihood with standard errors that tolerate non-normal data). cor = "polychoric" with estimator = "ML"/"MLR" errors, because lavaan itself does not support ML or MLR on ordered indicators. "WLSMV"/"ULSMV" with a continuous cor is allowed, and gives a valid, if atypical, continuous WLS or ADF estimator. Ignored for PCA and EFA engines. The effective value (after auto-selection) is recorded in x$meta$estimator (NA for PCA and EFA).

missing

How to handle missing item responses. One of:

"pairwise" (the default) uses all available observations pairwise. For PCA and EFA this feeds stats::cor(use = "pairwise.complete.obs"). For ESEM with WLSMV, and with ULSMV on cor = "polychoric", lavaan uses available.cases, which computes its sample statistics from all rows that contribute to each pair. That is valid under MCAR and uses the full N. For ESEM with ML, MLR, or ULSMV on a continuous cor, lavaan fits the complete rows, while edges are computed from a pairwise correlation matrix. The fields $meta$missing and $meta$n_complete record the setting and the number of complete rows. A warning is emitted when incomplete rows are detected.
"listwise" uses only complete rows. It reduces the data to stats::complete.cases() before fitting, so the correlation matrix, the engine fit, and the edges all use the same rows. Under engine = "esem" with cor = "spearman", lavaan still fits a Pearson matrix (see the description above). n_obs in the result reflects the reduced N.
"fiml" is full information maximum likelihood. For engine = "esem" (with estimator = "ML"/"MLR" and cor = "pearson"), it passes missing = "fiml" to lavaan::efa() and derives edge correlations from lavaan's FIML-estimated saturated model. For engine = "pca"/"efa" (M38), the correlation matrix is estimated via psych::corFiml() and fed to the usual W'RW algebra. That route is announced via a cli message. Both routes require cor = "pearson", because corFiml and lavaan's FIML each estimate a multivariate-normal matrix. FIML errors for WLSMV and ULSMV, and for a non-Pearson basis on any engine. FIML improves estimation under missingness but does not impute item responses. So score materialisation (keep_scores = TRUE) still produces NA rows for incomplete observations. See n_obs for the fit-index sample size on the PCA and EFA path.

n_obs

Number of observations. On the raw-data FIML path it may instead be a string selecting which N feeds the fit indices. FIML is full information maximum likelihood, which uses every observed value without dropping incomplete rows.

Correlation-matrix input: a positive integer. It is required for engine = "efa", because psych::fa() needs N for chi-square, RMSEA, and TLI. It is optional for "pca" (stored as NA_integer_ if omitted). For PCA it is recorded in the result metadata and feeds the N-based sampling-adequacy checks only. PCA level fit is eigenvalue-based and computes no N-dependent fit statistics.
Raw data: N is normally taken from nrow(data) and a numeric n_obs is ignored (with a warning). The exception is missing = "fiml" with engine = "pca"/"efa" (M38). psych::corFiml() estimates the correlation matrix from incomplete rows, so n_obs may be "total" or "complete". The default is "total", meaning every row that contributes to the FIML likelihood, which matches the FIML convention (Enders, 2010). "complete" is the complete-case N, a conservative lower bound. Point estimates do not depend on this choice. Those are the loadings (how strongly each item reflects a factor) and the edges. Only the fit indices of EFA (exploratory factor analysis) do. Those indices are approximate whatever N you pick, because of this two-step route: a FIML matrix fed into normal-theory EFA (Zhang & Savalei, 2020). A string n_obs is accepted only on this path.

align_signs

Logical. Sign-align factors to primary-parent lineage? Default TRUE.

keep_scores

Logical. Store factor scores in the result? Default FALSE (recomputable via augment.ackwards()). When TRUE, per-observation scores are stored in x$scores as a named list of n x k_j matrices, one per level, standardized by real score SDs (see augment.ackwards()).

keep_fits

Logical. Store raw engine fit objects? Default FALSE. When TRUE, the per-level fit objects (psych or lavaan) are stored in x$fits as a named list indexed by level.

seed

Integer seed for steps that start at random: every lavaan rotation (ESEM) and psych's oblimin rotation for EFA. Both rotate from several random starts and keep the best one. Without a seed, two oblimin EFA fits of the same data can reach a different solution at a deep level. It is not used by PCA but is captured for reproducibility metadata. Default NULL.

pairs

Which level pairs to compute edges for. "adjacent" is the default and is the classic Goldberg choice of consecutive levels only. "all" is the Forbes extension, which takes every pair of levels, including skip-level correlations. The prune() function recomputes its own all-pairs edges on demand whatever this setting is, so pruning does not require pairs = "all" here.

cut_show

Edges with |r| >= cut_show are flagged above_cut in tidy() output. Default 0.3.

correct

Continuity correction passed to psych::polychoric() on the PCA/EFA polychoric path (engine = "pca"/"efa" with cor = "polychoric"). Default 0.5 (psych's own default), which adds that value to zero cells before estimating thresholds. Set correct = 0 if psych::polychoric() fails on your data, and its error suggests exactly this. That failure typically happens when an item has a near-empty response category, or when items with unequal category counts produce a sparse cross-cell. It is ignored on other paths: ESEM computes its own polychoric correlations inside lavaan, and the Pearson/Spearman bases do not use it.

rotation

The rotation used at every level with two or more factors. The default, "varimax", keeps the factors within a level uncorrelated. The other values are oblique rotations, which let them correlate. With engine = "pca" or "efa", they are "oblimin" and "promax". With these engines "oblimin" needs the GPArotation package, and so does "promax" with engine = "efa". With engine = "esem", they are "oblimin" and "geomin" (lavaan's oblique geomin). With PCA, "promax" is stats::promax(). With EFA it is psych's Promax() through psych::kaiser(). The two can give slightly different loadings for the same solution.

Under an oblique rotation, each edge r is a total correlation. Primary parents and signs still use r, and tidy.ackwards() reports the partialled beta beside it. Each level's factor_cor holds the engine's factor correlation. EFA and ESEM scores use the oblique form of the ten Berge weights, and PCA scores are the oblique components, so the scores correlate as the factors do. The variance values follow psych's convention for correlated factors. The fit announces what its edges mean, and prune() warns that its thresholds were set for varimax. Two kinds of rotation start at random: psych's oblimin for engine = "efa" and every lavaan rotation. Set seed to reproduce such a fit exactly (see seed).

A rotation that fails is handled by how many starts it uses. PCA rotates once, so a level whose oblique rotation fails to converge, or that psych replaces with promax, ends the hierarchy at the level before it with a warning. So does a PCA level whose oblique fit raises an error. EFA (psych's fa() defaults to several starts) rotates from several random starts and keeps the best one, so an oblique rotation warning there can come from a discarded start. It is shown, and the level is kept. An EFA level whose final rotation step failed still ends the hierarchy, and so does any EFA rotation warning when the installed psych's fa() rotates from one start.

ESEM counts the failed starts under every lavaan rotation, varimax included. At each rotated level, lavaan keeps the best of 30 random starts. This is lavaan's default, which ackwards() does not change. If some starts did not converge, a warning gives the count (for example, "3 of 30"), and the level is kept. If no start converged, the kept rotation did not converge either, so the hierarchy ends at the level before it with a warning.

...

Reserved for future arguments.

Value

An object of class "ackwards". See print.ackwards(), tidy.ackwards(), glance.ackwards(), and augment.ackwards() for output methods.

Defaults and why

engine = "pca" is the original Goldberg (2006) method. PCA (principal component analysis) is the fastest engine and never fails to converge. Waller (2007) showed that the correlations between the component scores of different levels follow from the correlation matrix alone. (An oblique rotation can still fail, see rotation.)
rotation = "varimax" keeps the within-level factors mutually uncorrelated (orthogonal). A rotation re-orients the factors without changing how well they fit, and varimax pushes each item toward one factor. With uncorrelated factors, each between-level edge equals that ancestor's unique contribution to the descendant. That is what lets the edges be read as a lineage diagram. The closed-form W'RW edge algebra is itself exact for the stored x$r under any fixed linear scoring, orthogonal or not. So the choice is interpretive, not a numerical necessity. It matches Goldberg (2006), Kim & Eaton (2015), and Forbush et al. (2024). An oblique rotation is available as a non-default option (see rotation). Its edges are total correlations, which also carry overlap through correlated factors at the same level.
cor = "pearson" means no silent basis switching. If your items look ordinal (<= 7 distinct integer values), a cli warning will suggest cor = "polychoric", which is available for all three engines. Ordinal items have a few ordered categories, such as a 1 to 5 rating. Polychoric correlations estimate the correlation between the continuous traits assumed to underlie those ordered responses.
align_signs = TRUE because unaligned signs make the output unreadable. The anchor is m1f1, which is oriented toward the positive manifold. Each subsequent factor is flipped so its edge to its primary parent is positive.
keep_scores = FALSE / keep_fits = FALSE for memory and privacy. Scores are O(n x Sigmak) and often sensitive, and raw engine fits can be large. Both are recomputable from the stored r matrix.

Performance (ESEM, large item sets)

The ESEM engine fits a separate lavaan model at every level 1..k_max. For ordinal data (cor = "polychoric", WLSMV) lavaan derives costly sample statistics from the raw data. Those are the thresholds, the polychoric correlation matrix, and the asymptotic weight matrix. They depend only on the data, not on the number of factors. So they are computed once at the first level and reused for every deeper level, which gives identical solutions for much less work. This matters most when you have many items (hundreds), where recomputing those statistics at each level dominated the run time.

The per-level model fits are mutually independent and are dispatched through the future framework when future.apply is installed. By default the plan is sequential (no behaviour change). To run the levels in parallel, set a plan once before calling ackwards():


  future::plan(future::multisession, workers = 4)  # or multicore on Unix
  x <- ackwards(items, k_max = 8, engine = "esem", cor = "polychoric")

Parallelism pays off when the per-level fits are heavy (large p, several levels). For small problems the worker startup cost can outweigh it. Results are reproducible across plans when seed is supplied. PCA and EFA already compute their correlation matrix once and are unaffected.

Correlation-matrix input

When data is a pre-computed correlation matrix (square, symmetric, unit diagonal), ackwards() runs entirely from that matrix using the W'RW algebra, so no raw item responses are needed. This is useful when you have a published correlation table or a polychoric matrix computed externally.

Constraints and behaviour when a correlation matrix is supplied:

Engine: only "pca" and "efa" are supported. "esem" requires raw data (for lavaan's own polychoric computation, WLSMV estimation, and per-level fit indices) and will error clearly.
n_obs: required for "efa", because psych needs N for chi-square, RMSEA, and TLI. It is optional for "pca", where it is stored as NA if omitted and is used only for the N-based sampling-adequacy checks and the result metadata. PCA computes no N-dependent fit statistics.
cor argument: ignored, because the basis is already determined by the matrix you supply. A warning is emitted if you set cor explicitly.
missing argument: ignored, because missingness was handled when computing the matrix. A warning is emitted if you set missing explicitly.
Factor scores: keep_scores = TRUE will error. augment() and tidy(what = "scores") will also error because individual-level scores require row-level item responses.
$cor field: stored as NA_character_, and printed as "(user-supplied matrix)".

When to trust the result

ackwards() raises most diagnostics as it fits, while tidy() and summary() raise the warnings about a level's within-level score correlation (the item below says which raises which). They fall into three tiers by what they mean for whether you should trust and report the solution:

Fatal, so fix before trusting. The result is undefined or rests on a broken correlation matrix:

A constant item (no variance) errors. Drop it (see check_items()).
psych::polychoric() fails, usually because of a near-empty response category. Set correct = 0 or collapse rare categories.
A level fails to converge. The hierarchy is truncated to the deepest level that did converge, and you should not interpret beyond it.
A near-singular correlation matrix (smallest eigenvalue < 1e-4, recorded in meta$near_singular / meta$min_eigenvalue and re-surfaced by print()/summary()) means per-level fit indices and factor scores are unreliable. The loadings and edges then rest on a rank-deficient matrix. (For EFA the residual-based fallback inflates TLI and RMSEA, and for ESEM CFI comes back NA.) Trim redundant items, or use missing = "listwise". On the polychoric basis you can instead collapse sparse categories or set correct = 0.

Caution, so interpret carefully. The solution exists but may be unstable:

A Heywood case at a level, meaning a communality > 1 or a negative uniqueness. Check that level's loadings make substantive sense, and consider fewer factors.
A near-constant item (one response category dominates) can drive a meaningless factor. Inspect it with check_items().
Ordinal data on a Pearson basis attenuates correlations. That is fine for a quick look or suggest_k() screening, but report the final model on cor = "polychoric".
A nearly singular within-level score correlation (smallest eigenvalue < 1e-2) means the factors at that level are close to collinear, so the partialled beta for edges from that level is unstable. The tidy() method still reports beta and raises a warning that names the level. Fitting with ackwards() does not raise it. The package chose this cutoff as a numerical guard. It is not a published rule. If the correlation cannot be inverted at all, beta and r2 are NA and a different warning names the level. Both tidy() and summary() raise that one.

Informational, and usually fine. Proceed, just be aware of the pairwise-missing note, a merely sparse (rare-but-present) response category, and the ordinal-detection warning when you did intend Pearson.

References

Goldberg, L. R. (2006). Doing it all Bass-Ackwards: The development of hierarchical factor structures from the top down. Journal of Research in Personality, 40(4), 347–358. doi:10.1016/j.jrp.2006.01.001

Waller, N. G. (2007). A general method for computing hierarchical component structures by Goldberg's bass-ackwards method. Journal of Research in Personality, 41(4), 745–752. doi:10.1016/j.jrp.2006.08.005

Forbes, M. K. (2023). Improving hierarchical models of individual differences: An extension of Goldberg's bass-ackward method. Psychological Methods. doi:10.1037/met0000546

Enders, C. K. (2010). Applied Missing Data Analysis. Guilford Press.

Zhang, X., & Savalei, V. (2020). Examining the effect of missing data on RMSEA and CFI under normal theory full-information maximum likelihood. Structural Equation Modeling, 27(2), 219–239. doi:10.1080/10705511.2019.1642111

Kaiser, H. F. (1958). The varimax criterion for analytic rotation in factor analysis. Psychometrika, 23(3), 187–200. doi:10.1007/BF02289233

Grice, J. W. (2001). Computing and evaluating factor scores. Psychological Methods, 6(4), 430–450. doi:10.1037/1082-989X.6.4.430

Beauducel, A., Hilger, N., & Kuhl, T. (2024). The trade-off between factor score determinacy and the preservation of inter-factor correlations. Educational and Psychological Measurement, 84(2), 289–313. doi:10.1177/00131644231171137

Williams, A. L., Conway, C. C., Olino, T. M., Revelle, W., Zinbarg, R. E., & HiTOP Utility Workgroup. (2025). Testing criterion validity in hierarchical models of psychopathology: Comparison of latent-variable and factor-score approaches. Clinical Psychological Science, 13(1), 128–145. doi:10.1177/21677026231225414

Examples

# sim16 is continuous with a known 1 -> 2 -> 4 hierarchy, so the default
# pearson basis is appropriate. For ordinal items (e.g. bfi25), fit on the
# polychoric basis instead -- see `cor = "polychoric"` and the ordinal
# vignette.
x <- ackwards(sim16, k_max = 4)
print(x)
#> 
#> ── Bass-Ackwards Analysis (ackwards) ───────────────────────────────────────────
#> Engine: pca
#> Rotation: varimax
#> Basis: pearson
#> n: 1,000
#> k (max): 4
#> 
#> ── Levels ──
#> 
#> ✔ k = 1: 1 factor, 28.2% variance
#> ✔ k = 2: 2 factors, 46.5% variance
#> ✔ k = 3: 3 factors, 57.5% variance
#> ✔ k = 4: 4 factors, 67.7% variance
#> 
#> ── Edges ──
#> 
#> 9 of 20 edges have |r| ≥ 0.3
#> ────────────────────────────────────────────────────────────────────────────────
#> Note: This is a series of linked solutions, not a fitted hierarchical model.
#> Cross-level edges are descriptive score correlations. Per-level fit indices
#> (EFA/ESEM) describe how well a k-factor model fits the items at that level --
#> they do not validate the edges or the hierarchy itself.
tidy(x)
#>    from   to level_from level_to            r         beta is_primary above_cut
#> 1  m1f1 m2f1          1        2  0.707296113  0.707296113       TRUE      TRUE
#> 2  m1f1 m2f2          1        2  0.706917399  0.706917399       TRUE      TRUE
#> 3  m2f1 m3f1          2        3 -0.011342395 -0.011342395      FALSE     FALSE
#> 4  m2f1 m3f2          2        3  0.697092981  0.697092981       TRUE      TRUE
#> 5  m2f1 m3f3          2        3  0.716891014  0.716891014       TRUE      TRUE
#> 6  m2f2 m3f1          2        3  0.995593231  0.995593231       TRUE      TRUE
#> 7  m2f2 m3f2          2        3  0.074615767  0.074615767      FALSE     FALSE
#> 8  m2f2 m3f3          2        3 -0.056803218 -0.056803218      FALSE     FALSE
#> 9  m3f1 m4f1          3        4 -0.016458058 -0.016458058      FALSE     FALSE
#> 10 m3f1 m4f2          3        4 -0.004812835 -0.004812835      FALSE     FALSE
#> 11 m3f1 m4f3          3        4  0.663913488  0.663913488       TRUE      TRUE
#> 12 m3f1 m4f4          3        4  0.747612767  0.747612767       TRUE      TRUE
#> 13 m3f2 m4f1          3        4  0.940647736  0.940647736       TRUE      TRUE
#> 14 m3f2 m4f2          3        4  0.031735215  0.031735215      FALSE     FALSE
#> 15 m3f2 m4f3          3        4  0.262765620  0.262765620      FALSE     FALSE
#> 16 m3f2 m4f4          3        4 -0.212435736 -0.212435736      FALSE     FALSE
#> 17 m3f3 m4f1          3        4  0.051257892  0.051257892      FALSE     FALSE
#> 18 m3f3 m4f2          3        4  0.970778427  0.970778427       TRUE      TRUE
#> 19 m3f3 m4f3          3        4 -0.171584098 -0.171584098      FALSE     FALSE
#> 20 m3f3 m4f4          3        4  0.159752219  0.159752219      FALSE     FALSE
glance(x)
#>   engine rotation     cor k_max n_obs deepest_converged n_edges CFI TLI RMSEA
#> 1    pca  varimax pearson     4  1000                 4      20  NA  NA    NA
#>   SRMR BIC
#> 1   NA  NA

# Correlation-matrix input (PCA engine; n_obs optional)
R <- cor(sim16)
x_R <- ackwards(R, k_max = 4)
#> ℹ `n_obs` not supplied; stored as `NA`.
#> ℹ PCA level fit is eigenvalue-based and does not use N; supplying `n_obs = <N>`
#>   records it in the result metadata and enables the N-based sampling-adequacy
#>   checks.
print(x_R)
#> 
#> ── Bass-Ackwards Analysis (ackwards) ───────────────────────────────────────────
#> Engine: pca
#> Rotation: varimax
#> Basis: (user-supplied matrix)
#> n: NA
#> k (max): 4
#> 
#> ── Levels ──
#> 
#> ✔ k = 1: 1 factor, 28.2% variance
#> ✔ k = 2: 2 factors, 46.5% variance
#> ✔ k = 3: 3 factors, 57.5% variance
#> ✔ k = 4: 4 factors, 67.7% variance
#> 
#> ── Edges ──
#> 
#> 9 of 20 edges have |r| ≥ 0.3
#> ────────────────────────────────────────────────────────────────────────────────
#> Note: This is a series of linked solutions, not a fitted hierarchical model.
#> Cross-level edges are descriptive score correlations. Per-level fit indices
#> (EFA/ESEM) describe how well a k-factor model fits the items at that level --
#> they do not validate the edges or the hierarchy itself.