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. 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 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 restrictions apply. By method, 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 method restriction alone did not catch.

See also

extract_ordination() for scores and loadings together.

Examples

if (FALSE) { # \dontrun{
getLV(fit, level = "unit")
getLV(fit, level = "unit", se = TRUE)
} # }