Skip to contents

Per-unit covariance of ordination (latent) scores

Source: R/ordination-uncertainty.R
ordination_uncertainty.Rd

getLV(se = TRUE) reports a standard error for every latent-score cell, one axis at a time. ordination_uncertainty() reports the full covariance of a unit's score VECTOR across its axes – the object needed to draw an uncertainty ellipse (rather than an axis-aligned error bar) around a site in an ordination biplot.

Usage

ordination_uncertainty(fit, level = "unit")

Arguments

fit

A fitted multivariate model returned by gllvmTMB(). Must be a native TMB fit with gllvmTMBcontrol(se = TRUE) (the default); engine = "julia" bridge fits, integration = "va" fits, and estimator = "mspl" fits are refused (see Details).

level

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

Value

NULL if the fit has no latent() term at the requested level (matching extract_ordination()). Otherwise a list of class "gllvmTMB_ordination_uncertainty":

  • scores – the n x K score matrix (identical to getLV(fit, level), native rotate = "none" orientation).

  • se – an n x K matrix of per-(unit, axis) standard deviations; sqrt(diag(cov[, , s])) for row s, and numerically identical to getLV(fit, level, se = TRUE)$se.

  • cov – a K x K x n array; cov[, , s] is unit s's latent-score covariance matrix, dimnamed by axis on the first two margins and by unit on the third.

  • level – the canonical level name ("unit" / "unit_obs").

se and cov are filled with NA (with a warning) when the fit's Hessian is not positive-definite, matching getLV(se = TRUE).

Estimand (read this before using the output)

Latent scores are random effects, not fixed parameters. What this function returns is the model's best guess at how far each site's score could plausibly have landed given (a) the data actually observed at that site and (b) the package's uncertainty about the fitted parameters (loadings, dispersion, etc.) – not a sampling-distribution standard error of a fixed number, the way std.error on a regression coefficient is. Two model fits refit on two different datasets would not be expected to recover this quantity "on average" the way a coefficient SE is; it is closer to a prediction interval than a confidence interval. Practically: it answers "how sure is the model about THIS site's position on the ordination plot," not "how much would the whole ordination move if I recollected the data."

Concretely, it is the covariance of \(\mathbf u_s\) (site \(s\)'s latent score vector) given the data and the fitted model, read off the fit's joint (fixed effects + random effects) precision matrix (TMB::sdreport(getJointPrecision = TRUE)), so it propagates the package's uncertainty about the loadings and other fixed parameters into the score covariance, not just the Laplace curvature at the fixed-effect point estimate.

Rotation / identifiability

Scores are reported in gllvmTMB's native fitting orientation (rotate = "none") only. That orientation is pinned by a structural lower-triangular constraint on the loading matrix, which fixes the axes up to a per-axis SIGN FLIP (a discrete choice), not a continuous rotation. Per-axis variance (se) is unaffected by that sign choice. The cross-axis covariance in cov is the one that belongs with the returned scores from THIS fit, so an ellipse built from them is correctly oriented for this fit's own biplot – but do not compare the sign of a cross-axis entry, or of an axis, across independent (re)fits; that sign is this fit's own arbitrary convention, exactly as in ordinary PCA or factor analysis. Rotated scores (getLV(..., rotate = "varimax") / "promax") are not supported – the covariance's transform under a post-hoc rotation is not implemented, matching getLV(se = TRUE)'s existing rotate = "none"-only restriction.

What this is not

A Wald quantity read from asymptotic (delta-method / Laplace) theory, not a resampled, profiled, or simulation-based interval. Its repeated -sampling coverage has not been measured for this function; se and cov are reported as-is, with no lower/upper bound, so that no unmeasured coverage claim is implied by the return shape (the same choice getLV(se = TRUE) already makes). A confidence ellipse built from cov at some nominal level (e.g. via stats::qchisq(level, df = 2)) is a standard asymptotic-normal construction, not a certified interval.

See also

getLV() for the marginal per-cell standard error this function generalises to a full per-unit covariance; ordiplot()'s ellipse argument, which draws the covariance this function returns.

Examples

if (FALSE) { # \dontrun{
u <- ordination_uncertainty(fit, level = "unit")
u$se           # same numbers as getLV(fit, "unit", se = TRUE)$se
u$cov[, , 1]   # site 1's 2x2 (or dxd) score covariance
} # }