Sweep the loading-ridge penalty and report a stability path

ridge_path() refits the same model across a grid of loading-ridge penalty scales (tau, i.e. gllvmTMBcontrol(loading_ridge = tau)) and reports, for every trait, how the largest loading and its communality move as the penalty weakens (tau growing toward Inf, i.e. plain maximum likelihood with no ridge).

Usage

ridge_path(
  formula,
  data,
  family,
  tau = c(0.5, 1, 2, 4, 8, Inf),
  unit = NULL,
  trait = "trait",
  weights = NULL,
  missing = miss_control(),
  control = gllvmTMBcontrol(),
  ...
)

Arguments

formula, data, family, unit, trait, weights, missing: Passed to gllvmTMB() at every grid point, exactly as you would call it directly. unit has no default – it must name the sampling-unit column, matching gllvmTMB()'s own required argument.
tau: Numeric vector of ridge scales to sweep. tau is sorted into increasing order before refitting, so the returned rows and the printed table run from the strongest penalty to the weakest whatever order you give. Every entry must be a positive number; Inf is allowed and requests plain ML. Default c(0.5, 1, 2, 4, 8, Inf).
control: A gllvmTMBcontrol() object providing every OTHER setting (n_init, se, warn_runaway, ...). Its own loading_ridge / aghq_ridge value is overridden at every grid point and does not need to be set. The estimator is always "ml": loading_ridge cannot be combined with estimator = "mspl" (see gllvmTMB()).
...: Passed to gllvmTMB() at every grid point.

Value

A data frame of class c("gllvmTMB_ridge_path", "data.frame") with one row per tau x trait combination:

tau: The ridge scale for this row's refit.
trait: Trait name.
max_loading: The largest |loading| for this trait across the model's reduced-rank blocks at this tau.
communality: This trait's between-unit communality at this tau (from extract_communality(), level = "unit"); NA if unavailable.
logLik_at_map: The UNPENALISED log-likelihood evaluated at the penalised (MAP) point – not a maximum for any finite tau, so it is not comparable across tau as a model-selection criterion. At tau = Inf this is the ordinary maximum log-likelihood.
convergence: The optimiser's convergence code (0 = converged); NA if the refit errored.
fit_error: The error message if this grid point's refit failed outright; NA otherwise.

Details

This is a SENSITIVITY diagnostic, never an identification certificate. A ridge that stabilises a fit does not prove the underlying model is identified; it only prevents the optimiser from drifting to an unbounded point. The path lets you tell two situations apart that a single fit cannot: a trait whose loading and communality settle down as the penalty weakens is behaving like an ordinary, data-determined estimate; a trait that keeps moving toward the boundary as the penalty weakens is one whose finite estimate under any single tau is being determined by the penalty, not by the data.

Every grid point is a COLD refit: each tau starts from the model's own default initial values, not from the previous grid point's solution. This trades some speed for comparability across tau – a warm-started path could stay near one local optimum across the whole grid and never reveal that a later tau has a different, better optimum, which would understate exactly the instability this function exists to detect.

print() classifies each trait using a SIGNED statistic – a shrinking or flat loading is always "interior", regardless of magnitude; only a substantial POSITIVE move counts as "penalty-determined". Two comparisons are possible, chosen automatically per trait:

A converged tau = Inf point is available: The largest finite-tau loading is compared directly to the Inf loading (the plain signed relative change). A blow-up there – exactly what the paragraph above calls "moving toward boundary" – is compared against boundary_rel_thresh (default 0.10, i.e. a >10% increase from the last finite tau to Inf). This comparison is never skipped when a converged Inf point exists: it is the diagnostic the previous paragraph promises.
No converged tau = Inf point is available: The last two finite-tau grid points are compared using the log-log slope (elasticity) d(log loading) / d(log tau), not the raw relative change. The raw relative change decays like 1/tau for a fixed absolute grid spacing even when a trait is genuinely still moving – a tau = c(2000, 2100) pair (a 1.05x ratio) gives a still-growing trait only a 5% raw change, an artifact of the grid geometry rather than the trait's behaviour. The elasticity is invariant to that spacing (for loading ~ tau^k it recovers k exactly at any two grid points) and is compared against elasticity_thresh (default 0.10).

Both are heuristics, not formal tests; see print.gllvmTMB_ridge_path's arguments to adjust either threshold.

tau = Inf requests a plain maximum-likelihood refit with no ridge. The boundary comparison uses that refit only when its convergence code is 0. An Inf row that returns a loading whose convergence code is not 0 is not the boundary anchor, and the printed line says that the Inf refit did not converge. The elasticity comparison is then used only when at least two finite-tau rows remain. An outright error leaves that loading missing, so it is not the boundary anchor and does not add the sentence. A converged Inf loading, including an extreme one, stays the boundary comparison. fit_error records the message when a grid point's refit errors outright, and the row's other columns are NA for that point.

Examples

set.seed(1)
# Simulate presence/absence of five species at 200 sites. One shared
# latent gradient (say, habitat quality) drives all five; species 4
# responds negatively.
n_sites <- 200
gradient <- rnorm(n_sites)
loadings <- c(1.5, 1.2, 1.0, -1.0, 0.8)
intercepts <- c(0, 0.5, -0.5, 0.3, -0.2)
dat <- expand.grid(site = factor(seq_len(n_sites)),
                   trait = factor(paste0("sp", 1:5)))
eta <- intercepts[dat$trait] + loadings[dat$trait] * gradient[dat$site]
dat$value <- rbinom(nrow(dat), size = 1, prob = plogis(eta))

# Refit a one-factor model while a ridge penalty on the loadings
# weakens: tau = 0.5 shrinks hard, tau = Inf is plain maximum
# likelihood with no penalty. Each tau is a separate fit. With one
# presence/absence record per species and site, the model uses
# loadings only (unique = FALSE): extra per-species variance could not
# be told apart from the coin-flip noise of 0/1 data.
path <- ridge_path(
  value ~ 0 + trait + latent(0 + trait | site, d = 1, unique = FALSE),
  data = dat, family = binomial(), unit = "site", trait = "trait",
  tau = c(0.5, 2, 8, Inf)
)

# max_loading is the size of each species' loading. It levels off as
# tau grows, so every species is classed "interior": the data, not the
# penalty, set its size. A loading that kept growing as the penalty
# weakened would be flagged "penalty-determined".
path
#> gllvmTMB ridge-path sensitivity table
#>  tau trait max_loading communality logLik_at_map convergence fit_error
#>  0.5   sp1   0.8010259  0.16320507     -671.9534           0      <NA>
#>  0.5   sp2   0.7757318  0.15462933     -671.9534           0      <NA>
#>  0.5   sp3   0.5290655  0.07841111     -671.9534           0      <NA>
#>  0.5   sp4   0.6033933  0.09964101     -671.9534           0      <NA>
#>  0.5   sp5   0.7266561  0.13830364     -671.9534           0      <NA>
#>  2.0   sp1   1.0380132  0.24671116     -670.6290           0      <NA>
#>  2.0   sp2   1.0036789  0.23442290     -670.6290           0      <NA>
#>  2.0   sp3   0.6227023  0.10543707     -670.6290           0      <NA>
#>  2.0   sp4   0.7247513  0.13767922     -670.6290           0      <NA>
#>  2.0   sp5   0.9119916  0.20179769     -670.6290           0      <NA>
#>  8.0   sp1   1.0602297  0.25466661     -670.6197           0      <NA>
#>  8.0   sp2   1.0248824  0.24200990     -670.6197           0      <NA>
#>  8.0   sp3   0.6293424  0.10745479     -670.6197           0      <NA>
#>  8.0   sp4   0.7339021  0.14068565     -670.6197           0      <NA>
#>  8.0   sp5   0.9277157  0.20736089     -670.6197           0      <NA>
#>  Inf   sp1   1.0617841  0.25522316     -670.6197           0      <NA>
#>  Inf   sp2   1.0263334  0.24252931     -670.6197           0      <NA>
#>  Inf   sp3   0.6297877  0.10759053     -670.6197           0      <NA>
#>  Inf   sp4   0.7345361  0.14089458     -670.6197           0      <NA>
#>  Inf   sp5   0.9287953  0.20774345     -670.6197           0      <NA>
#> 
#> Per-trait classification (see ?ridge_path for the boundary/elasticity rule):
#>   sp1: interior (stabilises as penalty weakens) [tau 8 -> Inf: max|loading| 1.060 -> 1.062, relative change to Inf 0.1%]
#>   sp2: interior (stabilises as penalty weakens) [tau 8 -> Inf: max|loading| 1.025 -> 1.026, relative change to Inf 0.1%]
#>   sp3: interior (stabilises as penalty weakens) [tau 8 -> Inf: max|loading| 0.629 -> 0.630, relative change to Inf 0.1%]
#>   sp4: interior (stabilises as penalty weakens) [tau 8 -> Inf: max|loading| 0.734 -> 0.735, relative change to Inf 0.1%]
#>   sp5: interior (stabilises as penalty weakens) [tau 8 -> Inf: max|loading| 0.928 -> 0.929, relative change to Inf 0.1%]