Extract ordination scores and loadings from a fitted multivariate model

Usage

extract_ordination(
  fit,
  level = "unit",
  component = c("total", "innovation", "mean")
)

Arguments

fit: A fitted multivariate model returned by gllvmTMB(). Admitted engine = "julia" bridge fits expose raw unit-tier loadings and scores. Gaussian bridge fits with latent(..., lv = ~ x) also expose retained "mean" and "innovation" score components. Within-unit, structured-tier, and rotated ordinations remain gated for Julia bridge extractors. A gllvmTMB_va (integration = "va") fit exposes level = "unit" loadings and scores only, with generic trait/unit row and column names synthesised (a variational fit carries no stored data or column names); level = "unit_obs" returns NULL because the variational route fits no within-unit tier.
level: "unit" (between-unit) or "unit_obs" (within-unit). Deprecated aliases "B" and "W" are still accepted with a warning.
component: Score component to return. "total" returns the latent score entering the linear predictor. "innovation" returns the zero-mean latent innovation. "mean" returns the predictor-informed score mean and is non-zero only for latent(..., lv = ~ x) fits. Temporal ordination supports "total" only; other component requests raise an error.

Value

A list with scores (units or within-unit observations in rows, latent axes in columns), loadings (traits in rows, axes in columns), and row_id (the row identifiers of scores). The temporal branch also returns row_index (the pair table behind the rows of scores) and temporal_sign (the sign convention applied to the axes).

Examples

## Simulate 6 traits recorded on a few species at each of 80 sites.
## Between sites, the traits vary along two latent axes (the loadings
## in Lambda_B): traits 1-2 load mainly on the first, traits 4-5 on
## the second, and traits 3 and 6 on both. Each trait also has its
## own unique variance (psi_B), and sigma2_eps is the variance of the
## remaining noise in each observation. beta = 0 turns off the
## simulator's environmental effects, which this model does not
## include.
set.seed(5)
Lambda_B <- cbind(c(0.9, 0.7, 0.5, 0.0, 0.0, 0.3),
                  c(0.0, 0.3, 0.5, 0.9, 0.7, -0.4))
sim <- simulate_site_trait(
  n_sites = 80, n_species = 10, n_traits = 6,
  mean_species_per_site = 3,
  beta = matrix(0, 6, 2),
  Lambda_B = Lambda_B,
  psi_B = rep(0.3, 6),
  sigma2_eps = 0.3
)
## Two latent axes between sites, plus a unique variance for each
## trait (unique = TRUE)
fit <- gllvmTMB(
  value ~ 0 + trait +
    latent(0 + trait | site, d = 2, unique = TRUE),
  data  = sim$data,
  trait = "trait",
  unit  = "site"
)
## level = "unit" asks for the between-site ordination
ord <- extract_ordination(fit, level = "unit")
## Site scores: each site's predicted position on the two latent
## axes (one row per site, one column per axis)
dim(ord$scores)
#> [1] 80  2
head(ord$scores)
#>           LV1           LV2
#> 1  0.18674829  0.9933514612
#> 2 -0.02002702  0.0002150301
#> 3  0.30803244 -1.1372858332
#> 4 -1.08366032  0.5657920643
#> 5 -1.16119710 -0.3676723199
#> 6 -0.66365842 -0.1036673315
## Trait loadings: how strongly each trait follows each axis. A site
## with a high LV1 score tends to have high values of the traits with
## positive LV1 loadings. Rotating the axes (or flipping the sign of
## one) fits the data equally well; gllvmTMB fixes trait_1's loading
## on LV2 at zero, which settles the rotation but not the signs.
ord$loadings
#>                LV1        LV2
#> trait_1 0.69947489  0.0000000
#> trait_2 0.81301538  0.2602390
#> trait_3 0.64294477  0.2529062
#> trait_4 0.15181676  0.8959594
#> trait_5 0.08368232  0.6840032
#> trait_6 0.12552714 -0.4151906