Rotate the loadings of a fitted multivariate model

Applies a varimax or promax rotation to the loading matrix \(\Lambda\) from a fit returned by gllvmTMB(). Use level = "unit" for the between-unit reduced-rank component and level = "unit_obs" for the within-unit component. The latent scores are rotated by the complementary transform so the linear predictor, fitted log-likelihood, and implied covariance are unchanged.

Usage

rotate_loadings(
  fit,
  level = "unit",
  method = c("varimax", "promax", "none"),
  order_axes = TRUE,
  sign_anchor = c("auto", "none"),
  anchor_traits = NULL
)

Arguments

fit: A fitted multivariate model returned by gllvmTMB().
level: "unit" (between-unit) or "unit_obs" (within-unit). Deprecated aliases "B" and "W" are still accepted with a warning.
method: One of "varimax", "promax", or "none".
order_axes: Logical. When TRUE (default for rotated output), reorder rotated axes by decreasing shared variance colSums(Lambda^2). Ignored when method = "none".
sign_anchor: One of "auto" or "none". "auto" (default for rotated output) flips each rotated axis so its anchor trait has a positive loading. Ignored when method = "none".
anchor_traits: Optional character vector of trait names used for sign anchoring. Supply one trait per axis after ordering. Axes without a supplied anchor use the trait with the largest absolute loading.

Value

A list with rotated Lambda (n_traits × d), rotated scores (with rows = units or within-unit observations, columns = factors), and the rotation matrix T such that \(\Lambda_{\text{rotated}} = \Lambda T\). The list also includes axis_variance, axis_order, axis_sign, and anchor_traits metadata after any ordering and sign anchoring.

Details

Rotation is for interpretation of the loading columns, especially figures. It does not make the axes uniquely "right"; it chooses a readable orientation among mathematically equivalent orientations. For quantitative interpretation, start from the model-implied Sigma, correlations, communality, and uniqueness, then use rotated axes for labels and plots.

A defensible plotting workflow is: fit the model, check convergence, reconstruct covariance summaries, then rotate loadings for visual interpretation. Rotate level = "unit" and level = "unit_obs" separately because they represent different covariance structures. method = "varimax" is the usual first choice for a simple loading pattern; use method = "promax" only when correlated latent axes are intended.

The rotation is applied to \(\Lambda\) on the left and to the latent scores on the right using the inverse transform, so \(\Lambda_{\text{rot}} \mathbf{z}_{\text{rot}} = \Lambda \mathbf{z}\) is unchanged. For varimax, T is orthogonal so the rotation preserves the implied covariance \(\Lambda \Lambda^\top\); for promax, T is oblique and the columns of \(\Lambda_{\text{rot}}\) are no longer orthogonal.

Examples

set.seed(1)
# Simulate six traits, each measured twice on 100 units (e.g. sites).
# Traits 2-3 share one latent axis, traits 4-6 share a second axis,
# and trait_1 loads on both.
Lambda <- cbind(c(0.7, 1.0, 0.9, 0.0, 0.0, 0.0),
                c(0.7, 0.0, 0.0, 1.0, 0.9, 0.8))
sim <- simulate_unit_trait(n_units = 100, n_obs_per_unit = 2,
                           n_traits = 6, Lambda_B = Lambda,
                           psi_B = 0.4, sigma2_eps = 0.3)
# Two latent axes plus a trait-specific variance (unique = TRUE).
fit <- gllvmTMB(value ~ 0 + trait +
                  latent(0 + trait | unit, d = 2, unique = TRUE),
                data = sim$data, trait = "trait", unit = "unit")

# Raw loadings: to pin down the axes, the fit fixes trait_1's loading
# on axis 2 at zero, so the other traits are spread over both axes.
raw <- extract_ordination(fit, "unit")
round(raw$loadings, 2)
#>          LV1   LV2
#> trait_1 0.89  0.00
#> trait_2 0.73  0.59
#> trait_3 0.83  0.58
#> trait_4 0.56 -0.79
#> trait_5 0.55 -0.64
#> trait_6 0.63 -0.62

# Varimax rotation: traits 4-6 now load on one axis and traits 2-3 on
# the other, with trait_1 shared between them.
rot <- rotate_loadings(fit, level = "unit", method = "varimax")
round(rot$Lambda, 2)
#>           LV1   LV2
#> trait_1  0.53  0.71
#> trait_2 -0.04  0.94
#> trait_3  0.03  1.01
#> trait_4  0.97 -0.02
#> trait_5  0.84  0.06
#> trait_6  0.87  0.14

# Rotation only turns the axes. The scores turn with them, so each
# unit's latent contribution (scores times loadings) is unchanged.
head(round(rot$scores, 2))
#>     LV1   LV2
#> 1 -0.02  0.31
#> 2  0.89  0.55
#> 3  1.04  0.68
#> 4  1.50 -0.13
#> 5 -0.93  1.27
#> 6 -0.09  0.14
all.equal(rot$scores %*% t(rot$Lambda),
          raw$scores %*% t(raw$loadings))
#> [1] TRUE