Skip to contents

Extract latent-variable scores from a fitted multivariate model

Source: R/output-methods.R
getLV.Rd

Returns latent scores from a fit returned by gllvmTMB(). This is a small compatibility wrapper around extract_ordination() for readers familiar with gllvm::getLV().

Usage

getLV(fit, level = "unit", rotate = c("none", "varimax", "promax"), se = FALSE)

Arguments

fit

A fitted multivariate model returned by gllvmTMB(). Admitted engine = "julia" bridge fits expose raw unit-tier loadings and scores; rotated ordinations remain gated for Julia bridge fits.

level

"unit" (between-unit) or "unit_obs" (within-unit). Deprecated aliases "B" and "W" are still accepted with a warning.

rotate

Optional "varimax" or "promax" rotation after fitting. Default "none" returns the engine's native lower-triangular Lambda. For Julia bridge fits only "none" is currently routed.

se

Logical; if TRUE, also return uncertainty for every latent score: a frequentist standard error for an ordinary Laplace fit, or a variational posterior SD for a gllvmTMB_va fit. Default FALSE, which preserves the original behaviour (a bare matrix). See the Score uncertainty section.

Value

When se = FALSE (default): a matrix with one row per unit (level = "unit") or one row per within-unit observation (level = "unit_obs"), and one column per latent factor. When se = TRUE: a list with scores (that same matrix) and se (a matrix of identical shape and dimnames holding the score uncertainty described below).

Score uncertainty

For an ordinary Laplace fit, se = TRUE reads the marginal standard error of every unit-level (or within-unit) latent-score random effect – z_B at level = "unit", z_W at level = "unit_obs" – from the fit's TMB sdreport() (sqrt(sd_report$diag.cov.random)), and reshapes it with the identical matrix(..., nrow = d, ncol = n) then transpose convention that extract_ordination() uses for the point estimates, so scores[i, k] and se[i, k] always refer to the same (unit, axis) cell. This value is mathematically equivalent to inverting the fit's full joint precision matrix (TMB::sdreport(getJointPrecision = TRUE)) and reading the diagonal of the same block; the two routes were verified to agree to machine precision during development. Requirements:

  • The fit must carry a valid sdreport() (gllvmTMBcontrol(se = TRUE), the default) with a positive-definite Hessian; otherwise se = TRUE raises an error (no sdreport) or a warning with NA standard errors (non-positive-definite Hessian).

  • rotate must be "none": rotating scores changes their covariance, which is not currently propagated, so se = TRUE together with rotate != "none" raises an error rather than silently pairing rotated scores with un-rotated standard errors.

  • Predictor-informed latent(..., lv = ~ x) fits at level = "unit" are not yet supported (the score mean's own uncertainty is not yet propagated) and raise an error.

  • engine = "julia" bridge fits are not yet supported and raise an error.

For loading-ridge fits, score SDs require matching retained penalised-curvature metadata. They are local posterior approximations conditional on the ridge scale. Older or mismatched reports raise a warning and return NA SDs, while preserving the score point estimates.

For a gllvmTMB_va (variational) fit, se is instead the per-unit variational posterior SD read from the fit's own variational distribution at its optimum – not a Wald standard error, and not calibrated (the returned matrix carries a "uncertainty_basis" and "calibrated" attribute making this explicit in the object itself, not only here – though those attributes are silently dropped by se[i, ], head(se) and as.data.frame(se)).

Two independent gates restrict it. By tier, an explicitly requested "jj" (Jaakkola-Jordan) fit is refused pending its own measurement. By mechanism, any fit whose per-unit SD turns out to be constant across units is refused outright: the array would carry one row per unit while containing no per-unit information. That degeneracy is provable under "ac" (Albert-Chib) and was measured on a Gaussian fit (coefficient of variation 1.6e-15) which the tier gate alone did not catch.

See also

extract_ordination() for scores and loadings together.

Examples

set.seed(1)
## Five traits measured three times on each of 60 units. Between units,
## the traits share two latent axes (Lambda) plus a trait-specific
## unique variance (psi = 0.2 for every trait).
Lambda <- matrix(c(0.9, 0.7, 0.5, -0.3, 0.2,
                   0.0, 0.4, -0.5, 0.6, 0.8), nrow = 5, ncol = 2)
sim <- simulate_unit_trait(
  n_units = 60, n_obs_per_unit = 3, n_traits = 5,
  Lambda_B = Lambda, psi_B = rep(0.2, 5), sigma2_eps = 0.3
)
fit <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | unit, d = 2, unique = TRUE),
  data  = sim$data,
  trait = "trait",
  unit  = "unit"
)
## Latent-variable scores: each unit's predicted position on the two
## between-unit axes (one row per unit, one column per axis)
lv <- getLV(fit, level = "unit")
dim(lv)
#> [1] 60  2
head(lv)
#>           LV1         LV2
#> 1 -1.05207923  0.96595727
#> 2 -0.08325838 -1.76034400
#> 3  0.70690818  1.71992748
#> 4 -0.10253885  0.08801185
#> 5 -0.23204526  1.92088889
#> 6  0.85219464  0.39951765
## se = TRUE also returns the standard error of each score, from the
## fitted model's Wald (delta-method) approximation
lv_se <- getLV(fit, level = "unit", se = TRUE)
str(lv_se)
#> List of 2
#>  $ scores: num [1:60, 1:2] -1.0521 -0.0833 0.7069 -0.1025 -0.232 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:60] "1" "2" "3" "4" ...
#>   .. ..$ : chr [1:2] "LV1" "LV2"
#>  $ se    : num [1:60, 1:2] 0.494 0.476 0.467 0.453 0.471 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:60] "1" "2" "3" "4" ...
#>   .. ..$ : chr [1:2] "LV1" "LV2"