---
title: "Replicability-Gated Hierarchies: A Recommended Workflow"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Replicability-Gated Hierarchies: A Recommended Workflow}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---



The other vignettes in this package document its individual verbs. This one
takes a position: it lays out the workflow we recommend for answering the
questions bass-ackwards analysis exists to answer, and explains why each step
is there. We refer to it as the **replicability-gated** (or Girard) workflow,
because its central move is to let *replication*, not fit, decide how deep the
hierarchy you interpret goes.

## The problem: depth is where bass-ackwards analyses go wrong

A bass-ackwards analysis asks two questions: *how does the structure of this
domain organize hierarchically*, and *how deep is that hierarchy meaningful?*
The first question is what the method computes. The second is where published
applications most often stumble, and the failure mode is consistent:
**overextraction**, followed by substantive interpretation of deep-level
factors that would not re-emerge in a new sample. Forbes (2023) documents this
directly for the bass-ackwards context — non-replicable structure concentrates
in the deeper levels of an overextracted hierarchy.

The standard tooling half-addresses this. Retention criteria like parallel
analysis (`suggest_k()`) estimate a *plausible range* for k from the
eigenstructure, and they are the right first screen. But
they do not measure the thing the caution is actually about: whether the
specific factors in *your* solution would show up again in another sample from
the same population. A level can sit comfortably inside the plausible range
and still contain a factor that is pure sample idiosyncrasy.

There was, historically, a direct instrument for exactly this. Everett (1983)
proposed **factor comparability coefficients**: split the sample in half, fit
the solution in each half independently, apply *both* halves' scoring weights
to the full sample, and correlate the matched factor scores. A factor that is
real re-emerges in both halves and its two score estimates correlate near 1; a
factor that is noise does not. Split-half replication was a working gate in
the research program that produced the bass-ackwards method itself: Saucier
(1997) screened factor solutions by their stability across sample halves
(computing Everett's coefficients alongside), and Saucier, Georgiades,
Tsaousis, and Goldberg (2005) chose the optimal hierarchical level — the
number of factors — by requiring split-half replication above the
conventional .90 threshold. The modern bass-ackwards literature largely
dropped the practice.

`comparability()` restores that gate as a first-class verb, extended to the
bass-ackwards setting: coefficients for every factor at every level, computed
through the same score-correlation algebra as the hierarchy's own edges.

## Three verbs, three different questions

The workflow below leans on the fact that this package now separates three
questions that are easy to conflate:

| Question | Verb | What it measures |
|---|---|---|
| What depth range is *plausible*? | `suggest_k()` | Eigenstructure retention criteria (consensus range) |
| Which factors *replicate*? | `comparability()` | Split-half score comparability per factor per level |
| Which factors *differentiate*? | `prune()` | Forbes (2023) redundancy: factors that persist without adding resolution |

These are genuinely different. A factor can replicate perfectly and still be
redundant (the same construct restated at every level); a level can sit inside
the plausible range and still fail to replicate. You need all three answers,
and no one of them substitutes for another.

## Setup

We use the `bfi25` Big Five data (25 items, n = 875 after removing incomplete
rows), analyzed with the default PCA engine on the Pearson basis throughout.
That basis is deliberate: `comparability()` and `boot_edges()` screen on
Pearson/Spearman by design — polychoric estimation in every split-half or
bootstrap resample is slow and unstable — so the ordinal-detection warning
`bfi25` triggers is *expected* here and suppressed below for readability. You
screen structure and depth on Pearson, then fit the model you actually *report*
with `cor = "polychoric"`.


``` r
library(ackwards)
bfi <- na.omit(bfi25)
```

## Step 1: Screen the plausible depth range with `suggest_k()`


``` r
sk <- suggest_k(bfi, seed = 1)
#> ℹ Running parallel analysis (20 iterations, PC + FA)...
#> ✔ Running parallel analysis (20 iterations, PC + FA)... [134ms]
#> 
#> ℹ Running MAP and VSS...
#> ✔ Running MAP and VSS... [88ms]
#> 
#> ℹ Running Comparison Data (CD)...
#> ✔ Running Comparison Data (CD)... [5.5s]
#> 
print(sk)
#> 
#> ── Factor / Component Count Suggestion (ackwards) ──────────────────────────────
#> Variables: 25
#> n: 875
#> Basis: pearson
#> Tested k: 1-8
#> 
#> ── Criteria (k = 1-8) ──
#> 
#> k = 1: PA-PC ✔ PA-FA ✔ MAP 0.0254 VSS-1 0.5178 VSS-2 0.0000 CD ✔
#> k = 2: PA-PC ✔ PA-FA ✔ MAP 0.0194 VSS-1 0.5839 VSS-2 0.6719 CD ✔
#> k = 3: PA-PC ✔ PA-FA ✔ MAP 0.0175 VSS-1 0.5913 VSS-2 0.7354 CD ✔
#> k = 4: PA-PC ✔ PA-FA ✔ MAP 0.0164 VSS-1 0.6215* VSS-2 0.7837 CD ✔
#> k = 5: PA-PC ✔ PA-FA ✔ MAP 0.0160* VSS-1 0.5738 VSS-2 0.7950* CD ✔
#> k = 6: PA-PC - PA-FA ✔ MAP 0.0172 VSS-1 0.5594 VSS-2 0.7629 CD ✔*
#> k = 7: PA-PC - PA-FA - MAP 0.0205 VSS-1 0.5613 VSS-2 0.7616 CD -
#> k = 8: PA-PC - PA-FA - MAP 0.0236 VSS-1 0.5600 VSS-2 0.7215 CD -
#> 
#> ── Recommendations ──
#> 
#> • PA-PC: k <= 5
#> • PA-FA: k <= 6
#> • MAP: k = 5
#> • VSS-1: k = 4
#> • VSS-2: k = 5
#> • CD: k = 6
#> Consensus range: k = 4-6
#> ────────────────────────────────────────────────────────────────────────────────
#> Note: k_max in ackwards() is a maximum depth. Setting k_max one or two levels
#> above the consensus to observe factor fragmentation is intentional.
#> Caution: PA-PC tends to overextract; structures may not replicate (Forbes,
#> 2023). PA-FA and CD are more conservative. Use the range.
```



The criteria disagree — they always do — but they bracket a range (here
k = 4–6). Treat the top of that range as a ceiling worth
*probing*, not a depth worth *asserting*: parallel analysis on the PC basis in
particular tends to overextract, and `vignette("ackwards-suggest-k")` explains
each criterion's bias direction in detail. Deliberately looking one or two
levels past the consensus is a normal part of the method — the point of the
next step is that you can now do so safely.

## Step 2: Gate the depth on replicability with `comparability()`

This is the step the workflow is named for. We evaluate every level up to one
past the ceiling, across ten random split-halves:


``` r
cmp <- comparability(bfi, k_max = sk_hi + 1, n_splits = 10, seed = 2026)
#> ℹ Fitting 10 split-half replicates (pca, k = 1-7)...
#> ✔ Fitting 10 split-half replicates (pca, k = 1-7)... [377ms]
#> 
print(cmp)
#> 
#> ── Split-Half Factor Comparability (ackwards) ──────────────────────────────────
#> Engine: pca
#> Basis: pearson
#> n: 875 (437 per half)
#> Splits: 10
#> Levels: 1-7
#> 
#> ── Comparability by level (median across splits) ──
#> 
#> k = 1: median r 1.00, min r 1.00 (m1f1)
#> k = 2: median r .99, min r .98 (m2f2)
#> k = 3: median r .98, min r .92 (m3f3)
#> k = 4: median r .98, min r .95 (m4f4)
#> k = 5: median r .98, min r .97 (m5f4)
#> k = 6: median r .95, min r .88 (m6f5)
#> k = 7: median r .87, min r .58 (m7f7)
#> ────────────────────────────────────────────────────────────────────────────────
#> Per-factor detail (incl. Tucker's φ) in `$summary`; per-split values in
#> `$coefficients`.
#> Conventional benchmarks: ≥ .90 replication floor (Everett, 1983; Saucier et
#> al., 2005), ≥ .95 factors interchangeable (Lorenzo-Seva & ten Berge, 2006) --
#> conventions, not tests. Interpret levels whose factors all replicate.
```



Read the per-level minima, not the medians: a level is only as interpretable
as its least replicable factor. Here every factor through k = 5 clears
the conventional floor (the weakest median comparability at any of those
levels is .92), and the structure degrades
past it — at k = 7, m7f7 reaches a median of only
.58 and dips to -.24 on its worst
split. Those deep factors are not stable dimensions of the data; they are
different factors in every half-sample that happen to occupy the same slot.


``` r
autoplot(cmp)
```

![plot of chunk comparability-plot](assets/ackwards-girard-comparability-plot-1.png)

Two honesty notes. First, the conventional benchmarks marked on the plot
(.90, the replication threshold of Goldberg's lexical research program —
Everett, 1983; Saucier et al., 2005 — and .95, the stricter bound at which
factors are conventionally treated as interchangeable) are conventions,
not tests — `comparability()` reports every coefficient and flags nothing, so
borderline cases stay visible and the judgment stays yours. Second,
comparability is *internal* replication: it asks whether the structure
re-emerges in halves of the *same* sample. That is the right gate for "is this
factor real here," and a necessary condition for — but not a demonstration
of — replication in a new population.

The result is a **hierarchy floor**: the deepest level at which every factor
replicates. For these data that is k = 5 — the Big Five, as it should
be.

## Step 3: Fit the hierarchy to the gated depth


``` r
x <- ackwards(bfi, k_max = 5)
print(x)
#> 
#> ── Bass-Ackwards Analysis (ackwards) ───────────────────────────────────────────
#> Engine: pca
#> Rotation: varimax
#> Basis: pearson
#> n: 875
#> k (max): 5
#> 
#> ── Levels ──
#> 
#> ✔ k = 1: 1 factor, 20.9% variance
#> ✔ k = 2: 2 factors, 32.3% variance
#> ✔ k = 3: 3 factors, 40.7% variance
#> ✔ k = 4: 4 factors, 47.8% variance
#> ✔ k = 5: 5 factors, 53.8% variance
#> 
#> ── Edges ──
#> 
#> 14 of 40 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.
```


``` r
autoplot(x)
```

![plot of chunk hierarchy-plot](assets/ackwards-girard-hierarchy-plot-1.png)

If you prefer to *display* the fragmentation beyond the floor (it can be
substantively interesting to show a factor dissolving), fit one level deeper
and say explicitly in your report which levels passed the replicability gate.
What the gate rules out is interpreting the deeper factors as constructs.

## Step 4: Find what perpetuates without differentiating — `prune()`

Replicability is necessary but not sufficient: a factor can re-emerge in every
half-sample and still add nothing, because it is the same construct restated
level after level. That is Forbes's (2023) redundancy question, and it is
`prune()`'s job, not `comparability()`'s:


``` r
pr <- prune(x, "redundant")
#> ℹ Redundancy pruning (direct criterion, |r| ≥ 0.9) flagged 6 nodes.
#> ℹ Nodes are retained in the object; inspect with `x$prune$nodes` and
#>   `x$prune$chains`.
pr$prune$chains
#>    chain_id   id level r_to_prev phi_to_prev retain endpoint_r
#> 1         1 m2f2     2        NA          NA  FALSE  0.9648941
#> 2         1 m3f2     3 0.9703055   0.9780829  FALSE  0.9648941
#> 3         1 m4f2     4 0.9853622   0.9916584  FALSE  0.9648941
#> 4         1 m5f2     5 0.9990140   0.9995122   TRUE  0.9648941
#> 5         2 m3f1     3        NA          NA   TRUE  0.9978857
#> 6         2 m4f1     4 0.9978857   0.9988054  FALSE  0.9978857
#> 7         3 m4f3     4        NA          NA  FALSE  0.9773027
#> 8         3 m5f3     5 0.9773027   0.9872943   TRUE  0.9773027
#> 9         4 m4f4     4        NA          NA  FALSE  0.9570776
#> 10        4 m5f5     5 0.9570776   0.9671856   TRUE  0.9570776
#>    endpoint_r_agrees
#> 1               TRUE
#> 2               TRUE
#> 3               TRUE
#> 4               TRUE
#> 5               TRUE
#> 6               TRUE
#> 7               TRUE
#> 8               TRUE
#> 9               TRUE
#> 10              TRUE
```

Here the chains tell us that several of the Big Five arrive early and simply
persist — for example one factor emerges at k = 2 and travels essentially
unchanged (r ≥ .97 at every link) down to k = 5. Every node on that chain
*replicates* beautifully; the chain flags that the intermediate appearances
add no resolution. The two verbs answer different questions, and reporting
both gives the honest picture: which levels are real, and which levels are
new. See `vignette("ackwards-forbes")` for chain mechanics, thresholds, and
the retention rule.

### How sure are we about these edges? `boot_edges()`

`prune()`'s redundancy rule and the Forbes practice of reading off the
strongest skip-level correlation both consume a point estimate at a hard
threshold, so the natural follow-up is how precisely each edge is estimated.
`boot_edges()` attaches a nonparametric bootstrap standard error and percentile
confidence interval to every edge. Each replicate resamples respondents,
recomputes the correlations, refits the whole hierarchy, and — importantly —
re-anchors each replicate's factors to the full-sample solution (the same
matching and sign-orientation `comparability()` uses), so factor
label-switching across replicates does not contaminate the intervals. It runs
on the PCA and EFA engines with a Pearson or Spearman basis — a polychoric
matrix is too slow and unstable to re-estimate in every replicate (the same
scope as `comparability()`), so screen edge stability on the Pearson basis even
when the final model is polychoric.


``` r
x_boot <- ackwards(bfi, k_max = 5, pairs = "all") |>
  boot_edges(bfi, n_boot = 200, seed = 1)
#> ℹ Fitting 200 bootstrap replicates (pca, k = 1-5)...
#> ✔ Fitting 200 bootstrap replicates (pca, k = 1-5)... [2.9s]
#> 

boot_tbl <- tidy(x_boot, what = "edges", sort = "strength")
skip_boot <- boot_tbl[
  abs(boot_tbl$level_to - boot_tbl$level_from) > 1 & abs(boot_tbl$r) >= 0.5,
  c("from", "to", "r", "lo", "hi")
]
head(skip_boot, 6)
#>    from   to         r        lo        hi
#> 4  m3f2 m5f2 0.9852300 0.9541520 0.9984825
#> 7  m2f2 m4f2 0.9673808 0.8942028 0.9916158
#> 8  m2f2 m5f2 0.9648941 0.8883935 0.9893501
#> 12 m2f1 m4f1 0.8700141 0.7593997 0.9290617
#> 14 m3f3 m5f3 0.8084301 0.5312592 0.9636492
#> 15 m3f1 m5f1 0.8080606 0.6585810 0.9044025
```

The `lo`/`hi` columns are the 95% percentile interval. A skip-level edge whose
interval sits comfortably above `prune()`'s `redundancy_r` threshold (0.9 by
default) is a defensible "same construct" claim; one whose interval straddles
the threshold should not be treated as decisively redundant. The intervals are
per-edge error bars, not a familywise guarantee: they do not correct for having
*searched* the full edge table for the largest value, so a strongest-edge claim
should still be pre-specified rather than selected after the fact.

## Step 5: Interpret with the built-in guardrails


``` r
top_items(x, level = 5, cut = 0.5)
#> 
#> ── Salient items by factor (ackwards) ──────────────────────────────────────────
#> Engine: pca
#> Cut: |loading| >= 0.5
#> Top-n: all
#> 
#> ── Level 5 (5 factors) ──
#> 
#> m5f1
#> E2 [-0.720]
#> E4 [0.704]
#> E1 [-0.685]
#> E3 [0.659]
#> E5 [0.584]
#> 
#> m5f2
#> N3 [-0.799]
#> N1 [-0.787]
#> N2 [-0.782]
#> N5 [-0.656]
#> N4 [-0.626]
#> 
#> m5f3
#> C2 [0.702]
#> C4 [-0.692]
#> C1 [0.671]
#> C3 [0.653]
#> C5 [-0.627]
#> 
#> m5f4
#> A3 [0.698]
#> A2 [0.672]
#> A1 [-0.672]
#> A5 [0.569]
#> A4 [0.513]
#> 
#> m5f5
#> O5 [-0.676]
#> O3 [0.617]
#> O2 [-0.586]
#> O1 [0.554]
#> ────────────────────────────────────────────────────────────────────────────────
#> Loadings reflect primary-parent sign alignment. Use tidy(x, what = "loadings")
#> for the full matrix.
```

Three interpretation habits keep the analysis honest (all covered in depth in
`vignette("ackwards-interpret")`):

* **Negative loadings do not mean "low"** — sign is anchored to the primary
  parent, and a negative secondary edge is information, not error.
* **Name factors within levels, borrowing names down the lineage** rather than
  treating each level as a fresh exploratory result.
* **Say what the hierarchy is.** A bass-ackwards result is a series of linked
  solutions whose edges are score correlations — it is descriptive, not a
  fitted hierarchical model (no Schmid–Leiman, no higher-order SEM), and
  should be reported as such.

## Step 6: Validate downstream, out of sample

The final step closes the loop the same way it opened: with data the model
has not seen. Scores for new observations use the fit-time metric, so
training and test scores are directly comparable:


``` r
train <- bfi[1:600, ]
test <- bfi[601:875, ]

x_train <- ackwards(train, k_max = 5)
test_scores <- predict(x_train, newdata = test)
round(head(test_scores[, grep("^\\.m5", names(test_scores))], 3), 2)
#>   .m5f1 .m5f2 .m5f3 .m5f4 .m5f5
#> 1  0.21 -2.77 -1.21 -1.45 -0.27
#> 2  0.19 -1.20  0.28 -1.21 -1.41
#> 3 -0.91 -1.46  1.34  0.06 -0.98
```

If a downstream claim (a correlation with an outcome, a group difference)
holds for training-sample scores but not test-sample scores, the hierarchy was
fine and the claim was overfit — a different failure than anything in Steps
1–5, and one only held-out data can catch. The train/test mechanics are
covered in `vignette("ackwards-intro")`.

## Common mistakes this workflow prevents

**Interpreting non-replicable deep factors.** The classic failure. Without a
replicability instrument the deep levels are just *there*, printed with the
same authority as the robust ones. Step 2 makes their status measurable and
visible.

**Trusting a single split.** One split-half can flatter or slander a factor by
luck of the draw. `comparability()` defaults to ten random splits and reports
the spread; the minima across splits are as informative as the medians.

**Cherry-picking the strongest edge.** With `pairs = "all"`, the number of
between-level correlations grows fast, and reporting the most dramatic one
capitalizes on chance. Report full edge tables (or the complete diagram), and
treat any "strongest link" claim as descriptive rather than inferential.
`boot_edges()` (Step 4 above) attaches a bootstrap confidence interval to every
edge, which makes each estimate's precision visible — useful for judging whether
an edge clears `prune()`'s redundancy threshold — but per-edge intervals do not
undo the selection bias of scanning many edges for the largest.

**Mistaking persistence for structure.** A chain of near-1.0 correlations down
the hierarchy looks impressive but means the levels are restating one
construct. That is a `prune()` finding, and it should shorten your reported k
range, not lengthen your list of constructs.

**Mistaking replicability for validity.** A comparability of .99 means the
factor is a stable feature of these items in this population — not that it is
a good construct, and not that it generalizes beyond the population sampled.
Cross-population replication and external validation remain separate,
subsequent steps.

## Summary: the workflow in six calls

1. **`suggest_k(data)`** — screen the plausible depth range; treat its top as
   a ceiling to probe.
2. **`comparability(data, k_max = ceiling + 1)`** — find the hierarchy floor:
   the deepest level at which every factor replicates across split-halves.
3. **`ackwards(data, k_max = floor)`** — fit the hierarchy to the gated depth.
4. **`prune(x, "redundant")`** — identify factors that perpetuate without
   differentiating; report which levels add resolution.
5. **`top_items(x)` + `autoplot(x)`** — interpret with the sign, naming, and
   descriptive-hierarchy guardrails.
6. **`predict(x_train, newdata = test)`** — validate downstream claims on
   held-out data in the fit-time metric.

Steps 1, 3, 4, and 5 are the established method. Step 2 restores the gate the
method's own inventor applied to it, and Step 6 extends the same logic — *show
me it holds in data you did not fit* — to whatever you do with the scores.

## References

Everett, J. E. (1983). Factor comparability as a means of determining the
number of factors and their rotation. *Multivariate Behavioral Research*,
*18*(2), 197–218. https://doi.org/10.1207/s15327906mbr1802_5

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

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. https://doi.org/10.1016/j.jrp.2006.01.001

Lorenzo-Seva, U., & ten Berge, J. M. F. (2006). Tucker's congruence coefficient
as a meaningful index of factor similarity. *Methodology*, *2*(2), 57–64.
https://doi.org/10.1027/1614-2241.2.2.57

Saucier, G. (1997). Effects of variable selection on the factor structure of
person descriptors. *Journal of Personality and Social Psychology*, *73*(6),
1296–1312. https://doi.org/10.1037/0022-3514.73.6.1296

Saucier, G., Georgiades, S., Tsaousis, I., & Goldberg, L. R. (2005). The
factor structure of Greek personality adjectives. *Journal of Personality and
Social Psychology*, *88*(5), 856–875. https://doi.org/10.1037/0022-3514.88.5.856
