---
title: "Axes Reliability"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Axes Reliability}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---




``` r
library(circumplex)
```



## 1. What is axis reliability?

A circumplex instrument places its scales around a circle and summarizes a
person (or a profile of correlations) by their position on two orthogonal
**axes** — here communion (the X axis, at 0°) and agency (the Y axis, at 90°).
Because those axis scores drive everything downstream — a person's projected
location, the displacement and amplitude of a Structural Summary Method
profile — it is worth asking how *reliably* the instrument measures each axis.

`axes_reliability()` answers that question with the estimator of Strack, Jacobs,
and Grosse Holtforth (2013). It fits an item-level measurement model that
decomposes each item's variance into orthogonal pieces — a general factor common
to all items, the two circumplex **axes**, a scale-specificity component, item
error, and, for an instrument administered in blocks, a block-specificity
component — and reads axis reliability off the **axes** component alone. Reliability
is then the Spearman–Brown "list-length" reliability of a composite of that
length built from items that share only their axes variance.

This is a different question from the other reliability-adjacent tools in the
package. `ssm_sem()` disattenuates a *scale-level* SSM profile for measurement
error; `fit_structure()` evaluates whether a correlation matrix has circumplex
*structure* at all. `axes_reliability()` instead reports a single, interpretable
number per axis: how well the instrument measures communion and agency.

## 2. A worked example

The package ships `simulated_items`, a synthetic dataset of 1–7 Likert
responses from 500 respondents on 32 items — four items on each of the eight
octant scales, in the order that `octants()` returns. The items were drawn from
a five-component population (a general factor, two equal axes with axes variance
.18, a shared scale-specificity component of .10, no block specificity — the
instrument is not blockwise, so that component is zero here — and free item
error) that implies an axis reliability of about .78.

`axes_reliability()` needs three things: the data, a map from items to scales,
and the scales' angles. The map is a list with one character vector of item
names per scale, in the same angle order as the angles you pass:


``` r
data("simulated_items")

# Four items per octant scale, in octants() order (PA, BC, DE, ..., NO).
items <- split(names(simulated_items), rep(1:8, each = 4))

res <- axes_reliability(simulated_items, items = items, angles = octants())
#> axes_reliability(): 500 complete case(s) used.
res
#> 
#> Circumplex Axes Reliability (Strack, Jacobs & Grosse Holtforth, 2013)
#> Input:        item data
#> Items:        32 (8 scales)
#> Complete N:   500
#> SEm scale:    std
#> 
#> # Per-axis reliability
#> 
#>  Axis item_n Reliability SEm   NB_Reliability
#>  X    16     0.773       0.476 0.822         
#>  Y    16     0.773       0.476 0.823         
#> 
#>   Note: the two axes share one axes-variance estimate and, with equal
#>   items per axis, carry the same reliability -- expected, not an error.
#> 
#>   Note: the model is fit to the item correlation matrix as if it were a
#>   covariance matrix (Cudeck, 1989), and both sides of that mismatch are
#>   corrected -- so these numbers differ from LISREL's, and from lavaan's
#>   own, by design.
#>   The component standard errors are adjusted to the correlation metric
#>   and are calibrated; they are typically smaller than the values printed
#>   by Strack et al. (2013), whose LISREL output carries no correction.
```

If your items belong to one of the package's built-in instruments, you can pass
the `instrument` object instead of `items` and `angles` — it supplies both the
scale angles and the item membership, exactly as `score()` does. (The example
above uses the explicit map because `simulated_items` is not a registered
instrument.)

The header confirms how many complete cases were used, and the per-axis table
reports, for each axis, the effective test length (`item_n`), the Strack axis
`Reliability`, its standard error of measurement (`SEm`), and the
Nunnally–Bernstein reliability (`NB_Reliability`) for comparison. For a balanced
instrument the two axes share one axes-variance estimate and carry equal
`item_n`, so they report the same reliability — expected, not an error.

The recovered reliability (about .77) lands close to the .78 built into the
simulated population, and the axes-variance estimate (below) recovers the
population value of .18.

## 3. Reading the components

`summary()` adds the estimated variance components and the model's global fit:


``` r
summary(res)
#> 
#> Circumplex Axes Reliability (Strack, Jacobs & Grosse Holtforth, 2013)
#> Input:        item data
#> Items:        32 (8 scales)
#> Complete N:   500
#> SEm scale:    std
#> 
#> # Per-axis reliability
#> 
#>  Axis item_n Reliability SEm   NB_Reliability
#>  X    16     0.773       0.476 0.822         
#>  Y    16     0.773       0.476 0.823         
#> 
#>   Note: the two axes share one axes-variance estimate and, with equal
#>   items per axis, carry the same reliability -- expected, not an error.
#> 
#>   Note: the model is fit to the item correlation matrix as if it were a
#>   covariance matrix (Cudeck, 1989), and both sides of that mismatch are
#>   corrected -- so these numbers differ from LISREL's, and from lavaan's
#>   own, by design.
#>   The component standard errors are adjusted to the correlation metric
#>   and are calibrated; they are typically smaller than the values printed
#>   by Strack et al. (2013), whose LISREL output carries no correction.
#> 
#> # Variance components
#> 
#>  Component         Estimate SE   
#>  general           0.051    0.005
#>  axes              0.175    0.009
#>  scale_specificity 0.093    0.008
#>  item              0.680    --   
#> 
#> # Global fit
#> 
#>   chi-square(493) = 488.27,  RMSEA = 0.000,  CFI = 1.000
#> 
#>   The global fit statistics chisq, pvalue, rmsea and cfi are scaled to
#>   that metric (Satorra & Bentler, 1994), which removes a distortion that
#>   flatters fit; df and srmr are unchanged. The scaled test can modestly
#>   over-reject at typical sample sizes -- it over-flags misfit rather than
#>   flattering it; see ?axes_reliability for the measured rates.
#>   They follow lavaan's *.scaled definitions, not its *.robust ones, and
#>   differ from what fitMeasures() reports for an equivalent ML fit.
```

The **variance components** show the decomposition the reliability rests on:
the `axes` component is the only one that feeds reliability, while `general`,
`scale_specificity`, `block_specificity` (when blocks were supplied), and item
error are isolated from it. This is precisely why
the Nunnally–Bernstein figure printed alongside runs **higher** than the Strack
reliability: N–B charges scale-specificity variance to the axis rather than
isolating it, so it **overestimates** axis reliability whenever scale
specificity is non-trivial (Strack et al., 2013, Figure 3). The gap between the
two numbers is a direct read-out of how much scale-specific variance the simpler
formula would have miscredited to the axes.

## 4. Starting from a published correlation matrix

You do not always have the raw data. A paper may print an item correlation
matrix and nothing else, and that matrix is enough: pass it as `cormat` together
with the sample size it was computed from, in place of `data`.


``` r
R <- cor(simulated_items)
axes_reliability(
  cormat = R, items = items, angles = octants(), n = nrow(simulated_items)
)
#> 
#> Circumplex Axes Reliability (Strack, Jacobs & Grosse Holtforth, 2013)
#> Input:        correlation matrix
#> Items:        32 (8 scales)
#> Sample N:     500
#> SEm scale:    std
#> 
#> # Per-axis reliability
#> 
#>  Axis item_n Reliability SEm   NB_Reliability
#>  X    16     0.773       0.476 --            
#>  Y    16     0.773       0.476 --            
#> 
#>   Note: the two axes share one axes-variance estimate and, with equal
#>   items per axis, carry the same reliability -- expected, not an error.
#> 
#>   Note: the Nunnally-Bernstein comparison needs the raw item scores
#>   (scale alphas and the axis-composite variance), so it is NA on the
#>   correlation-matrix path.
#> 
#>   Note: the model is fit to the item correlation matrix as if it were a
#>   covariance matrix (Cudeck, 1989), and both sides of that mismatch are
#>   corrected -- so these numbers differ from LISREL's, and from lavaan's
#>   own, by design.
#>   The component standard errors are adjusted to the correlation metric
#>   and are calibrated; they are typically smaller than the values printed
#>   by Strack et al. (2013), whose LISREL output carries no correction.
```

The estimates are identical to the raw-data run above — the raw-data path builds
exactly this matrix internally and fits it the same way. `items` selects and
orders the matrix's rows by name, so its own column ordering does not matter,
and it must be symmetric, positive definite, and have a unit diagonal (the model
assumes unit-variance items).

Two things are unavailable here, and both for the same reason: they are
properties of the respondents, not of their correlations. The
Nunnally–Bernstein comparison is reported as `NA` — it needs each scale's alpha
and the axis composite's variance, neither of which a correlation matrix
carries. And `sd = "raw"` is refused, because there are no scale scores to take
an observed SD from; supply the axis SDs numerically if you want SEm on a raw
scale. Both are reported rather than silently omitted, so a matrix-based result
cannot be mistaken for a raw-data one.

## 5. Caveats to keep in mind

Four properties of the method shape how its output should be read.

**Both the standard errors and the global fit statistics are corrected for the
correlation metric.**
Following the paper's own practice, the model is fit to the item
**correlation** matrix as though it were a covariance matrix. That metric
mismatch leaves the component point estimates and the reliabilities correct,
but it misprices anything computed from the input's sampling variability
(Cudeck, 1989).

The component standard errors are **corrected** for it. Normal-theory maximum
likelihood prices its standard errors for a sample **covariance** input, while
this estimator consumes a sample **correlation** matrix, whose diagonal cannot
vary at all. Left uncorrected that mismatch **overstates** sampling variability
by about 40% for an instrument whose axes carry a lot of variance, and
**understates** it slightly for weak-axes, strong-general instruments — it
changes sign across the range of instruments the function accepts, which is why
it could not be left to a caveat. The reported SEs are adjusted to the
correlation metric and are calibrated uncertainty. They are typically
**smaller** than the standard errors printed in Strack et al. (2013), whose
LISREL values carry the uncorrected approximation. What the estimator reported
before the correction is kept in `details$se_uncorrected`.

The **global fit statistics are corrected too**, for the same mismatch running
the other way. Sample correlations vary *less* than the covariances the
reference chi-square distribution is derived for, so an uncorrected test
statistic comes out too small and fit is flattered. `chisq`, `pvalue`, `rmsea`
and `cfi` are therefore reported as Satorra–Bentler-type **scaled** values —
the chi-square divided by a factor computed at the fitted matrix, with `cfi`
also scaling its own baseline model (Satorra & Bentler, 1994). The factor is
recomputed for every fit rather than being a constant, which is the whole
point: how much fit was flattered depends on the instrument. `df` and `srmr`
are unchanged, being a count of restrictions and a residual summary rather
than test statistics with a reference distribution. What the estimator reported
before the scaling is kept in `details$fit_uncorrected`, and the factors
themselves in `details$scaling_factor`.

If you go looking for these numbers in lavaan, match the variant. They are
built with the definitions lavaan calls `chisq.scaled`, `pvalue.scaled`,
`rmsea.scaled` and `cfi.scaled` — the mean-adjusted Satorra–Bentler forms — and
**not** the `*.robust` ones (`cfi.robust`, `rmsea.robust`), which apply a
different (Brosseau-Liard/Savalei) adjustment and give different numbers.
Because the model is estimated with plain ML and the scaling is applied by this
package rather than by lavaan, `fitMeasures()` on an equivalent fit hands back
the *uncorrected* values under the bare names `chisq`, `pvalue`, `rmsea` and
`cfi` — and reports no `*.scaled` or `*.robust` measure at all, since lavaan
supplies those only for a genuinely scaled estimator such as `"MLM"` or
`"MLR"`, and quietly returns a shorter vector when asked for one it does not
have. So a cross-check against lavaan's bare `cfi` will disagree with
`$fit$cfi`, a request for `cfi.robust` will come back empty rather than
disagreeing, and neither outcome is a defect.

Read the scaled statistic for what it is. The correction makes the test
statistic match its reference distribution *in expectation*; it is a
calibration, not an exactness guarantee, and it will not rescue a model that is
badly wrong. The usual cutoffs (Hu & Bentler, 1999) are no more or less
appropriate than before — what has changed is that the number being compared to
them is no longer systematically kind.

**How well calibrated is the test, and at what sample size?** The scaling fixes
the metric error, and the χ² test built on it is asymptotically exact: its
rejection rate approaches the nominal α as the number of distinct moments,
p\* = p(p+1)/2, falls relative to N. Simulated at one population (eight octant
scales, three items each, axes variance .35), the rejection rate at α = .05 runs

| p\*/N | 0.50 | 0.25 | 0.12 | 0.06 |
|---|---|---|---|---|
| rejection rate | .092 | .079 | .062 | .054 |

reaching the nominal band by a p\*/N of about 0.06. That is a sweep at a single
population, not a general threshold.

You do not have to work out where your own fit sits on that row. The object
reports both halves of the ratio: `details$n_moments` is p\*, and `details$n`
is the N the fit was priced at — on the raw-data paths the count of rows the
estimator was actually handed, smaller than `details$n_total` wherever cases
were dropped; on the correlation-matrix path the `n` you supplied. Divide the
first by the second and read across.

At **N = 600** the test **over-rejects** — measured .06 to .11 at three
populations chosen to bracket the range of instruments the function accepts.
The uncorrected statistic under-rejects over the same range, at .02 to .03, and
moves *further* from nominal as N grows, because its error is asymptotic while
the scaled statistic's is a finite-sample one that shrinks away. The
over-rejection at a fixed N grows with instrument size (larger `df`) and shrinks
with N.

Two practical consequences. A p-value near whatever threshold you are using
deserves caution at moderate N with many items — but note the direction: the
scaled test **over-flags** misfit rather than flattering it, which is the safer
error, and the opposite of what the uncorrected statistic did. And all of this
evidence is **complete-data**: under `missing = "fiml"` the scaled statistic is
calibrated in mean, but its rejection rate has not been measured, so none of
these rates should be read as applying to that path.

If you check the implied matrix, you will also find that it does
not reproduce the unit diagonal exactly. That too is expected: with the loadings
fixed, the condition a free item error satisfies is the *weighted* diagonal
rather than the raw one, so off-diagonal sampling misfit leaks into the implied
diagonal at roughly the sampling standard error of a correlation.

**Missing data: listwise by default, FIML on request.** The default,
`missing = "listwise"`, uses complete cases only and reports how many there
were. `missing = "fiml"` instead estimates from every respondent who answered
at least one item, and reports the total N with the complete-case count beside
it. Pairwise-deletion correlations are never used on either setting.

The choice is a trade of assumptions, not a free upgrade. FIML requires the
data to be **missing at random** — missingness may depend on values you
observed, but not on the unobserved values themselves — **and multivariate
normal**. Under MCAR, where missingness is unrelated to anything, listwise
deletion is already consistent and merely wasteful, so FIML buys precision
rather than correctness; under MAR, listwise deletion is genuinely biased and
FIML is not. Under FIML the standard errors are observed-information standard
errors on the standardized metric, conditional on the standardization
constants. They carry the same correlation-metric correction as every other
path, applied multiplicatively so the observed information's own pricing of the
missing data survives it. What the correction does not reach is the uncertainty
in those standardization constants, which stays small at mild missingness but
grows with the rate — at 15% cellwise missingness the reported SEs average
about 7% *below* the estimator's actual sampling variability, so treat heavy
missingness as the regime where they are least trustworthy. Two results become unavailable, both needing items observed by every
respondent: the Nunnally-Bernstein comparison is `NA` with a stated reason, and
`sd = "raw"` is refused in favour of numeric axis SDs.

Note the provenance, because it differs from the rest of this vignette: Strack
et al. (2013) report no missing-data analyses at all, so nothing about the FIML
path rests on their results. It is certified against this package's own
synthetic oracle, where the true variance components are known by construction.

**A boundary fit returns `NA`, not a clipped value.** If the model estimates an
axes variance outside the interval (0, 1) — at or below zero the axes carry no
variance to be reliable, at or above one they carry all of it — or any negative
variance component, the reliability and SEm are reported as `NA` with a warning
and a boundary flag, rather than a clipped, negative, or missing number. An `NA`
here is a signal that the model did not identify a usable axes-variance
component in your data — not a defect to be worked around.

**A blockwise instrument needs its blocks declared.** Some circumplex
instruments are administered in blocks — items grouped by something other than
their scale — which carries a block-specificity component of its own, reported
as high as 6.7% by Strack et al. (2013, Table 3). Pass `blocks` (a list of item
columns, one element per block, exactly as `items` is a list per scale) and that
component is estimated too: the component table gains a `zeta2` row and
`details$zeta2_fitted` is `TRUE`. Nothing in the bundled instrument objects
records block membership, so the map has to come from you.

The component is only estimable when the blocks say something the model does not
already know. Blocks that coincide with the scales, a single block holding every
item, or one block per item all leave `zeta2` unidentified, and it is dropped
with `details$zeta2_fitted` set to `FALSE` — the same treatment scale
specificity gets on a single-item instrument.

What it costs to *ignore* real blocks depends on their geometry. The general
factor never hands block variance back, so `xi2` is inflated under most layouts
and unchanged under a few — never deflated. The axes variance, the quantity
reliability is actually read from, moves only when block membership carries
information about the *angular distance* between items, beyond what sharing a
scale already tells you.

One case is clean enough to rely on: when **each block draws exactly one item
from every scale**, every within-block pair is a different-scale pair and the
blocks cover every pair of scale positions equally often. Membership then says
nothing about angular distance, and `xi1`, the reliability, and the SEm are
untouched — worth estimating for its own sake, but costless to omit.

Elsewhere, be careful: *"the blocks are spread evenly around the circle" is not
the test.* Blocks pairing diametrically opposite scales are as evenly spread as
blocks get — their angles average to the centre of the circle — and at eight
scales they still pull `xi1` about 9% *below* truth, because every within-block
pair sits exactly half a turn apart, which is very much information about
angular distance. Blocks covering contiguous arcs pull about 12% the other way.
Unless each block takes one item per scale, estimate the component rather than
reasoning about the picture.

Finally, a note on the `SEm`. The standard error of measurement supports a
location interval for a single profile (Strack et al., 2013, use ±1.65·SEm). By
default `axes_reliability()` reports the z-standardized SEm, `sqrt(1 -
reliability)`; pass `sd = "raw"` (or your own axis SDs) to put the SEm on the
raw axis-score scale. Such an interval describes the measurement imprecision of
one profile's axis position; it is not a significance test of that position
against any particular value.

## 6. Wrap-up

`axes_reliability()` gives a compact, per-axis answer to "how reliably does this
instrument measure communion and agency?", isolating the axes variance from the
general and scale-specific components that a simpler reliability formula would
conflate. Use it to characterize a circumplex instrument before leaning on its
axis scores, and read its output with the correlation-as-covariance,
missing-data, and boundary caveats in mind.

The examples above all use the canonical eight octant scales, but nothing in
the model requires them: any **equally spaced** set of angles works, at any
rotation and at any count from four scales upward. Equal spacing is what
matters, and it is required rather than merely preferred — a quasi-circumplex,
whose scales sit at slightly unequal intervals, is refused rather than
approximated, because Strack et al. (2013) excluded such instruments when
validating the model. Three scales are refused for a different reason: at that
count every pair of scales sits the same angular distance apart, so the
general, axes, and scale-specificity variances can no longer be told apart.

Scales may also carry a **single item** each, as Strack et al.'s single-item
circumplex types do. With one item at every position no two items ever share a
scale, so nothing in the data distinguishes scale specificity from item
specificity: the scale-specificity component is not identified, and
`axes_reliability()` drops it from the model rather than estimating a quantity
the data cannot support. The components table then carries three rows instead
of four, and `details$zeta1_fitted` records the drop. An instrument that
*mixes* single-item and multi-item scales still estimates the component — one
multi-item scale supplies the information, and the shared-value restriction
carries it to the remaining scales. One consequence follows for the comparison
figure: coefficient alpha is undefined for a one-item scale, so the
Nunnally–Bernstein reliability is reported as `NA` with its reason whenever any
scale has fewer than two items, exactly as Strack et al. leave that column
blank for such instruments.

One further boundary is worth stating plainly, because the paper crosses it and
this function does not: the model implemented here is **two**-dimensional.
Strack et al. also analyze a *spherical* instrument, SYMLOG, whose items span
three orthogonal dimensions; the figures they publish for it come from that
three-axis model, and no configuration `axes_reliability()` accepts will
reproduce them.

## References

* Cudeck, R. (1989). Analysis of correlation matrices using covariance
  structure models. _Psychological Bulletin, 105_(2), 317–327.

* Satorra, A., & Bentler, P. M. (1994). Corrections to test statistics and
  standard errors in covariance structure analysis. In _Latent variables
  analysis: Applications for developmental research_ (pp. 399–419).

* Strack, S., Jacobs, K. A., & Grosse Holtforth, M. (2013). The reliability of
  circumplex axes. _SAGE Open, 3_(2). https://doi.org/10.1177/2158244013486115
