ackwards: Bass-ackwards hierarchical structural analysis

In ackwards: Bass-Ackwards Hierarchical Structural Analysis

Bass-ackwards hierarchical structural analysis

Description

Extracts factor/component solutions at levels 1 through k, then characterises the hierarchy by computing between-level factor-score correlations. 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,
  ...
)

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" (ESEM requires raw data), and the missing and cor arguments are ignored. See the Correlation-matrix input section below.
k_max	Maximum number of factors/components to extract. Required; use suggest_k() if uncertain. Sets the depth of the hierarchy: levels 1 through k_max are all extracted and retained. (Note: suggest_k()'s k_max means something related but distinct – the max number of factors/components it 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; recommended for the clinical/HiTOP workflow (Kim & Eaton, 2015; Forbush et al., 2024). Requires lavaan >= 0.6-13.
cor	Correlation basis: "pearson" (default), "spearman", or "polychoric". For PCA/EFA, "polychoric" computes a polychoric correlation matrix via psych::polychoric() (requires psych). For ESEM, it triggers WLSMV estimation via lavaan. A cli warning is emitted when ordinal-looking columns are detected and cor is not "polychoric".
fm	Factor extraction method passed to psych::fa(); only used when engine = "efa". One of "minres" (default, robust OLS), "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" (robust ML). cor = "polychoric" with estimator = "ML"/"MLR" errors (lavaan itself does not support ML/MLR on ordered indicators); "WLSMV"/ "ULSMV" with a continuous cor is allowed (a valid, if atypical, continuous WLS/ADF estimator). Ignored for PCA and EFA engines. The effective value (after auto-selection) is recorded in x$meta$estimator (NA for PCA/EFA).
missing	How to handle missing item responses. One of: "pairwise" (default) – use all available observations pairwise. For PCA/EFA this feeds stats::cor(use = "pairwise.complete.obs"). For ESEM with WLSMV/ULSMV (ordinal), lavaan uses available.cases, which computes polychoric thresholds and correlations from all rows that contribute to each pair – MCAR-valid and uses the full N. For ESEM with ML/MLR (continuous), lavaan uses listwise deletion internally while edges are computed from a pairwise correlation matrix; this minor inconsistency is documented in ⁠$meta⁠. A warning is emitted when incomplete rows are detected. "listwise" – only complete rows are used. Reduces data to stats::complete.cases() before fitting, so the correlation matrix, the engine fit, and the edges are all consistent. n_obs in the result reflects the reduced N. "fiml" – Full Information Maximum Likelihood. For engine = "esem" (with estimator = "ML"/"MLR"), 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; this requires cor = "pearson" (corFiml estimates a multivariate-normal matrix) and the route is announced via a cli message. Errors for WLSMV/ULSMV and for a non-Pearson PCA/EFA basis. Note: FIML improves estimation under missingness but does not impute item responses; score materialisation (keep_scores = TRUE) still produces NA rows for incomplete observations. See n_obs for the fit-index sample size on the PCA/EFA path.
n_obs	Number of observations, or (on the raw-data FIML path) a string selecting which N feeds the fit indices. Correlation-matrix input: a positive integer. Required for engine = "efa" (psych::fa() needs N for chi-square / RMSEA / TLI); 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): because psych::corFiml() estimates the correlation matrix from incomplete rows, n_obs may be "total" (default – every row contributing to the FIML likelihood, matching the FIML convention; Enders, 2010) or "complete" (complete-case N, a conservative lower bound). Point estimates (loadings, edges) do not depend on this choice; only the EFA fit indices do, and those are approximate under this two-step (FIML matrix into normal-theory EFA) route regardless of N (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 stochastic engines (not used by PCA but captured for reproducibility metadata). Default NULL.
pairs	Which level pairs to compute edges for: "adjacent" (default, classic Goldberg – only consecutive levels) or "all" (Forbes extension – every pair of levels, including skip-level correlations). prune() recomputes its own all-pairs edges on demand regardless of this setting, 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 (its error suggests exactly this) – this typically happens when an item has a near-empty response category or when items with unequal category counts produce a sparse cross-cell. Ignored on other paths: ESEM computes its own polychoric correlations inside lavaan, and the Pearson/Spearman bases do not use it.
...	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" – the original Goldberg (2006) method; fastest; never fails to converge; the Waller (2007) algebra is exact for components.

• rotation = "varimax" – the ⁠T'=T^-1⁠ property of orthogonal rotation enables the closed-form ⁠W'RW⁠ edge algebra and keeps within-level factors uncorrelated so cross-level edges reflect only the hierarchical signal. Matches Goldberg (2006), Kim & Eaton (2015), and Forbush et al. (2024). Varimax is the only supported rotation; oblique rotation would confound the between-level signal that is the method's core output.

• cor = "pearson" – 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.

• align_signs = TRUE – unaligned signs make the output unreadable. Anchor: m1f1 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 – memory and privacy. Scores are O(n x Sigmak) and often sensitive; 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) the costly sample statistics lavaan derives from the raw data – thresholds, the polychoric correlation matrix, and the asymptotic weight matrix – 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 (identical solutions, 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 – 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" (psych needs N for chi-square / RMSEA / TLI); optional for "pca" (stored as NA if omitted; used for the N-based sampling-adequacy checks and result metadata only – PCA computes no N-dependent fit statistics).

• cor argument: ignored – the basis is already determined by the matrix you supply. A warning is emitted if you set cor explicitly.

• missing argument: ignored – 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_; printed as "(user-supplied matrix)".

When to trust the result

ackwards() raises diagnostics as it fits. They fall into three tiers by what they mean for whether you should trust and report the solution:

Fatal – 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 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; do 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 and the loadings/edges rest on a rank-deficient matrix. (For EFA the residual-based fallback inflates TLI/RMSEA; for ESEM CFI comes back NA.) Trim redundant items, use missing = "listwise", or – on the polychoric basis – collapse sparse categories or set correct = 0.

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

• A Heywood case (communality ⁠> 1⁠ / negative uniqueness) at a level – check that level's loadings make substantive sense; 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 – fine for a quick look or suggest_k() screening, but report the final model on cor = "polychoric".

Informational – usually fine. Proceed, just be aware: 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. \Sexpr[results=rd]{tools:::Rd_expr_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. \Sexpr[results=rd]{tools:::Rd_expr_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. \Sexpr[results=rd]{tools:::Rd_expr_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. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1080/10705511.2019.1642111")}

See Also

print.ackwards(), tidy.ackwards(), glance.ackwards(), prune() (Forbes-extension redundancy/artifact flagging, piped off the result of this function)

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)
tidy(x)
glance(x)

# Correlation-matrix input (PCA engine; n_obs optional)
R <- cor(sim16)
x_R <- ackwards(R, k_max = 4)
print(x_R)

ackwards documentation built on July 25, 2026, 1:08 a.m.