--- title: "Ordinal Data: Polychoric Correlations and WLSMV" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Ordinal Data: Polychoric Correlations and WLSMV} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- Most psychological scales use ordinal response formats: Likert-type items with a handful of ordered categories, binary (yes/no) items, and everything in between. An ordinal item is recorded on an ordered scale with few categories. Treating these items as continuous and computing Pearson correlations between them is common in practice. But it introduces systematic bias that can distort both the loadings and the between-level edges that bass-ackwards analysis depends on. A loading is the correlation between an item and a factor, and a factor is an unobserved dimension that explains why a set of items correlate. This vignette explains the problem and how `ackwards` addresses it. ## Why Pearson underestimates for ordinal items When a continuous latent trait is sliced into a small number of ordered categories, the observed correlation between two items drops. It is always lower than the correlation between their underlying latent variables. The coarser the rating scale, the larger the attenuation. For a 5-point scale, Pearson correlations can be 10 to 20% lower than the true latent correlations. For a 2-point (binary) scale, the underestimate can exceed 30%. This matters for bass-ackwards analysis because: 1. **Loadings are attenuated.** Factors appear weaker than they are. 2. **Between-level edges are attenuated.** The hierarchy looks flatter. Factors that should be nearly perfectly correlated across levels may appear only moderately correlated. 3. **The number of factors is underestimated.** Attenuation compresses the eigenvalue spectrum, so selection criteria may suggest fewer factors than truly exist. Parallel analysis, which keeps a factor only when it explains more variance than random data of the same size would, is one such criterion, and MAP is another. ## The polychoric solution A **polychoric correlation** between two ordinal items estimates the Pearson correlation that *would* exist between their underlying continuous latent variables, if those variables had been measured without discretization. It does this by fitting a bivariate normal model with thresholds at each category boundary. For binary (two-category) items this same estimate is traditionally called a *tetrachoric* correlation, but it is not a separate method. Tetrachoric is just the two-category special case of polychoric. `cor = "polychoric"` handles binary items automatically: `psych::polychoric()` detects the two-level case and returns the tetrachoric estimate, so you never request tetrachoric separately. Setting `cor = "polychoric"` in `ackwards()` replaces step 1 of the pipeline, computing the item correlation matrix `R`, with a polychoric matrix. All downstream computation then operates on that polychoric `R`: factor extraction, rotation, tenBerge scoring, and the between-level algebra. A rotation re-orients the factors within a level without changing how well they fit. ## Automatic detection `ackwards()` checks for ordinal-looking data using a simple heuristic. If any column is integer-valued with 7 or fewer distinct levels, the data is flagged, and a warning is issued when `cor = "pearson"` (the default). ``` r library(ackwards) bfi <- na.omit(bfi25) ``` ``` r # The warning fires with the default Pearson basis x_pearson <- ackwards(bfi, k_max = 5) ``` The warning is ackwards reminding you to make an active choice. Specifying `cor = "polychoric"` acknowledges the data type and suppresses it: ``` r x_poly <- ackwards(bfi, k_max = 5, cor = "polychoric") ``` ## Seeing the difference ### Correlation matrices The most direct way to see the impact of the basis is to compare the item correlations under each basis. Rather than read two 5×5 matrices side by side, the chart below plots the ten unique item-pair correlations among the five neuroticism items (`N1` to `N5`), Pearson against polychoric.
Dodged bar chart comparing Pearson and polychoric correlations for the ten unique N1-N5 item pairs; every polychoric bar is taller than its Pearson counterpart.

plot of chunk r-compare

Every polychoric bar sits above its Pearson counterpart: the correlations are consistently higher. The `N1` and `N2` pair, for instance, goes from about 0.73 (Pearson) to 0.79 (polychoric). Across the block the increase runs to roughly 0.05, enough to shift loadings and change eigenvalues. ### Loadings The table below compares primary loadings, the loading of each Neuroticism item on its dominant factor, under both correlation bases at k = 5. The Δ column is the attenuation: how much larger each loading becomes (in absolute value) once the polychoric basis removes the coarse-category suppression. Using |Polychoric| − |Pearson| keeps the attenuation positive regardless of loading sign.
Neuroticism-item loadings at k = 5: Pearson vs polychoric
Positive Δ = attenuation removed by the polychoric basis
Item Factor
Primary loading
Pearson Polychoric Δ (|Poly| − |Pearson|)1
N1 m5f2 −0.79 −0.81 0.02
N2 m5f2 −0.78 −0.81 0.02
N3 m5f2 −0.80 −0.82 0.03
N4 m5f2 −0.63 −0.65 0.02
N5 m5f2 −0.66 −0.69 0.03
1 Factor assignment and sign verified to match between bases. Delta uses |Poly| - |Pearson| so attenuation is positive regardless of loading sign.
Polychoric loadings for the Neuroticism items are noticeably higher. The latent structure is sharper when the attenuating effect of coarse categories is removed. ### Between-level edges The same comparison for the primary-parent edges. The Δ column is the change in connection strength (|r|). It is positive where the polychoric basis recovers a stronger connection by undoing attenuation, and negative in the few cases where the recovered structure is slightly looser. Using absolute values makes the direction read correctly even for the negatively-signed edge.
Primary-parent edges: Pearson vs polychoric
Δ = change in |r|; positive = stronger connection under polychoric basis
From To
Edge strength (r)
Pearson Polychoric Δ (|Poly| − |Pearson|)1
m1f1 m2f1 0.88 0.89 0.01
m1f1 m2f2 0.48 0.46 −0.02
m2f1 m3f1 0.90 0.87 −0.02
m2f2 m3f2 0.97 0.99 0.02
m2f1 m3f3 0.43 0.48 0.05
m3f1 m4f1 1.00 0.99 0.00
m3f2 m4f2 0.99 0.98 −0.01
m3f3 m4f3 0.80 0.73 −0.07
m3f3 m4f4 0.59 0.68 0.09
m4f1 m5f1 0.82 0.84 0.02
m4f2 m5f2 1.00 1.00 0.00
m4f3 m5f3 0.98 0.98 0.01
m4f1 m5f4 0.58 0.55 −0.03
m4f4 m5f5 0.96 0.99 0.03
1 NA in either column means the bases disagree on the primary parent for that factor.
The edges are broadly similar in sign and rank order, so the hierarchy is the same, but polychoric edges are stronger. Factors that represent genuinely stable dimensions show near-perfect correlations with their counterparts at adjacent levels when the true latent structure is properly recovered. ## WLSMV for ESEM with ordinal items When `engine = "esem"` and `cor = "polychoric"` are combined, `ackwards()` automatically switches the lavaan estimator to **WLSMV** (diagonally weighted least squares, mean- and variance-adjusted). ESEM (exploratory structural equation modeling) fits the factor model in lavaan and adds standard errors. WLSMV is the standard estimator for structural equation models with categorical or ordinal indicators, and it is the same method Mplus uses by default. WLSMV genuinely operates on the polychoric basis. The lavaan package estimates the item thresholds and the polychoric correlations among the ordinal indicators. It then fits the factor model to *that* matrix using a diagonal weight matrix. So combining `cor = "polychoric"` with the ESEM engine is not stacking two separate corrections. The polychoric structure is precisely what WLSMV is built to model, and the between-level edges are computed from lavaan's own latent correlation matrix (not a separately estimated `psych` polychoric matrix). ``` r x_esem <- ackwards(bfi, k_max = 3, engine = "esem", cor = "polychoric") x_esem #> #> ── Bass-Ackwards Analysis (ackwards) ─────────────────────────────────────────── #> Engine: esem #> Rotation: varimax #> Basis: polychoric #> n: 875 #> k (max): 3 #> #> ── Levels ── #> #> ✔ k = 1: 1 factor, 23.5% variance #> ✔ k = 2: 2 factors, 32.9% variance #> ✔ k = 3: 3 factors, 39.5% variance #> #> ── Edges ── #> #> 5 of 8 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. ``` The WLSMV fit indices (CFI, RMSEA, SRMR) are now computed on a model that respects the ordinal measurement level. They will typically look better than Pearson-based fit indices on the same data. That reflects both the better model specification and the different reference distribution WLSMV uses. ``` r tidy(x_esem, what = "fit", format = "wide") #> level chi dof p_value CFI TLI RMSEA SRMR BIC #> 1 2 3616.834 251 0 0.7117569 0.6554864 0.1238665 0.09547997 NA #> 2 3 2448.356 228 0 0.8098533 0.7498069 0.1055573 0.07172502 NA ``` See `vignette("ackwards-engines")`, section "Per-level fit: what it tells you (and what it doesn't)", for how to interpret these indices in the bass-ackwards context and how to produce the `autoplot(what = "fit")` trajectory chart. ## Practical guidance: when to use polychoric | Scale characteristics | Recommendation | |-----------------------|---------------| | Continuous or many categories (> 7) | `cor = "pearson"` (default) | | 5 to 7 ordered categories (typical Likert) | `cor = "polychoric"` (recommended) | | 3 to 4 categories | `cor = "polychoric"` (strongly recommended) | | Binary (2 categories) | `cor = "polychoric"` (essential, this is the tetrachoric case) | Polychoric estimation requires fitting a bivariate normal model for every item pair. With p items that is p(p−1)/2 pairs, about 300 for a 25-item scale. Computation is usually fast (seconds), but it scales as O(p²) and can be slow for very large item pools. **Screen skewed or sparse scales first.** Clinical symptom scales and other right-skewed items often have a response category endorsed by only a handful of people. That can make `psych::polychoric()` fail under its default continuity correction (set `correct = 0`). With many such items it can also drive the correlation matrix *near-singular*, in which case per-level fit indices (especially CFI, which comes back `NA`) and factor scores become unreliable. A factor score is each person's estimated position on a factor. Run [`check_items()`][ackwards::check_items] before fitting to see which items are degenerate or sparse. Then read the **"When to trust the result"** section of `?ackwards` for how each warning maps to whether you should report the solution. When a matrix is near-singular, `ackwards()` says so once and records it on the object, so `print()` and `summary()` keep reminding you. [`factorability()`][ackwards::factorability] is the companion screen for the correlation matrix as a whole. It reports the KMO measure of sampling adequacy (overall and per item), Bartlett's test of sphericity, the `N:p` ratio, and the Ledermann bound on how many factors `p` items can identify. Call it on the same `cor` basis you plan to fit. Read its output as *conventions, not verdicts*, because the KMO bands and `N:p` rules of thumb are contested. But a very low KMO or a handful of low-MSA items is a useful early signal. The `ackwards()` function runs a light version of this screen internally. It warns only at the consequential extreme (KMO below .5, `N:p` below 5, or `k_max` past the Ledermann bound for EFA/ESEM). EFA (exploratory factor analysis) models the items as a small number of shared factors plus item-specific noise. > **A note on score computation.** The short answer: yes, you can use the > factor scores from an ordinal analysis in downstream work. The caveat below > is about *scaling*, not validity. When `cor = "polychoric"`, the weight > matrices stored in the object are derived from the polychoric `R`. But > `augment()` materializes scores by applying those weights to > Pearson-standardized raw data (`.standardize(data)`). The resulting scores > are calibrated to *model-implied* unit variance rather than exactly > empirical unit variance. So their empirical standard deviations will be > close to but not precisely 1. That scaling detail is the only reason > `augment()` warns when `cor != "pearson"`. The warning does not mean the > scores are biased or unusable. The between-level edges from > `tidy(what = "edges")` are unaffected either way, because they come from > the `W′RW` algebra on the polychoric `R`, not from materialized scores. The > edges are exact for that `R`. ## References Flora, D. B., & Curran, P. J. (2004). An empirical evaluation of alternative methods of estimation for confirmatory factor analysis with ordinal data. *Psychological Methods*, *9*(4), 466--491. Muthén, B. O. (1984). A general structural equation model with dichotomous, ordered categorical, and continuous latent variable indicators. *Psychometrika*, *49*(1), 115--132. Olsson, U. (1979). Maximum likelihood estimation of the polychoric correlation coefficient. *Psychometrika*, *44*(4), 443--460. Pearson, K. (1900). Mathematical contributions to the theory of evolution. VII. On the correlation of characters not quantitatively measurable. *Philosophical Transactions of the Royal Society of London, Series A*, *195*, 1--47.