Select a latent-variable rank by information criterion

Fits the same model at latent rank (number of ordination axes) d = 1, ..., d_max by sweeping the d argument of the formula's single ordinary latent(...) covstruct term, and reports AIC, BIC, and AICc for each. The fit at the criterion-minimising rank is reported as selected_d.

select_lv() is a information-criterion workflow tool, not a hypothesis test: it says nothing about statistical significance and carries no interval on the chosen d. For a likelihood-ratio test of one additional latent dimension (with the appropriate boundary correction, and an explicit refusal when that correction is not justified), see anova.gllvmTMB_multi().

Usage

select_lv(formula, data, ..., d_max, criterion = c("bic", "aic", "aicc"))

# S3 method for class 'gllvmTMB_select_lv'
print(x, ...)

Arguments

formula: A gllvmTMB model formula containing exactly one ordinary latent(...) term (the between-unit reduced-rank ordination; latent(0 + trait | unit, ...)). Any d = argument already present in that term is overwritten for each swept value; write it without d =, or with any placeholder value, since it will be replaced. Structured source-specific terms (phylo_latent(), spatial_latent(), kernel_latent(), animal_latent()) are not swept and are rejected if present, because their own d is a different rank than the one this function selects.
data: A data frame, as passed to gllvmTMB().
...: Currently unused.
d_max: Single positive integer: the largest rank to try. Fitting is only identifiable up to the number of traits p (a p-row loading matrix cannot have rank greater than p); d_max greater than p is rejected before any fitting is attempted, naming p.
criterion: One of "aic", "bic" (default), or "aicc". Selects the row that minimises this criterion among fits that converged with a positive-definite Hessian; a fit that did not is excluded from selection (and reported with a warning) even though its row still appears in the table.
x: A "gllvmTMB_select_lv" object.

Value

An object of class "gllvmTMB_select_lv", a list with:

table: A data.frame with one row per attempted d: d, npar, logLik, aic, bic, aicc, converged (optimizer convergence flag), pd_hessian, seconds, and error (the error message when a fit failed, else NA).
selected_d: The chosen rank under criterion.
criterion: The criterion used for selection.
fits: A named list (names = d) of the fitted gllvmTMB() objects that succeeded; NULL for d that failed.
selected_fit: The fit at selected_d.

Scope

Every fit in the sweep uses ordinary maximum likelihood (REML = TRUE is rejected — AIC/BIC comparisons across different random-effect structures under REML are not meaningful when the induced conditioning changes with d, and gllvmTMB's REML route is Gaussian-only regardless). A single failing k (an error from gllvmTMB(), non-convergence, or a non-PD Hessian) does not abort the sweep; it is recorded in the table with NA criteria and excluded from selection, with a warning naming which d failed and why.

Examples

# \donttest{
set.seed(7)
n_units <- 70L
traits <- paste0("t", 1:4)
units <- paste0("u", seq_len(n_units))
Lambda_true <- matrix(
  c(0.8, 0.5, -0.6, 0.4, 0.3, 0.7, -0.5, 0.6),
  nrow = 4, ncol = 2
)
scores <- matrix(rnorm(n_units * 2), n_units, 2)
eta <- tcrossprod(scores, Lambda_true)
dat <- do.call(rbind, lapply(seq_len(n_units), function(i) {
  data.frame(unit = units[i], trait = traits,
    value = eta[i, ] + rnorm(4, sd = 0.2))
}))
dat$unit <- factor(dat$unit, levels = units)
dat$trait <- factor(dat$trait, levels = traits)

sel <- select_lv(
  value ~ 0 + trait + latent(0 + trait | unit, d = 1),
  data = dat, unit = "unit", trait = "trait", d_max = 3,
  criterion = "bic",
  control = gllvmTMBcontrol(optimizer = "optim", optArgs = list(method = "BFGS"))
)
#> ℹ Auto-suppressing `sigma_eps`: `indep(0 + trait | unit)` is at the per-row
#>   level, so it already absorbs the observation residual.
#> • Fixed at 0.000848 (~1/1000 of sd(y)) to keep the Gaussian density
#>   well-defined; the row-level residual variance is fully captured by the
#>   per-row diagonal term.
#> ℹ Auto-suppressing `sigma_eps`: `indep(0 + trait | unit)` is at the per-row
#>   level, so it already absorbs the observation residual.
#> • Fixed at 0.000848 (~1/1000 of sd(y)) to keep the Gaussian density
#>   well-defined; the row-level residual variance is fully captured by the
#>   per-row diagonal term.
#> Warning: NaNs produced
#> ℹ Auto-suppressing `sigma_eps`: `indep(0 + trait | unit)` is at the per-row
#>   level, so it already absorbs the observation residual.
#> • Fixed at 0.000848 (~1/1000 of sd(y)) to keep the Gaussian density
#>   well-defined; the row-level residual variance is fully captured by the
#>   per-row diagonal term.
#> Warning: NaNs produced
#> Warning: 2 of 3 fit(s) were excluded from selection (still shown in the table, criteria
#> NA where unavailable).
#> ℹ d = 2: Hessian is not positive-definite
#> ℹ d = 3: Hessian is not positive-definite
sel
#> gllvmTMB latent-rank selection (criterion = bic, selected d = 1)
#>    d npar   logLik     AIC     BIC    AICc conv pdHess
#>  * 1   12 -155.733 335.465 379.083 336.634 TRUE   TRUE
#>    2   15 -141.104 312.208 366.730 314.026 TRUE  FALSE
#>    3   17 -141.104 316.208 378.000 318.544 TRUE  FALSE
sel$selected_d
#> [1] 1
# }