Package 'gllvmTMB'

Description

Transform longitude/latitude columns to an equal-distance UTM coordinate reference system with sf::st_transform().

Usage

add_utm_columns(
  dat,
  ll_names = c("longitude", "latitude"),
  ll_crs = 4326,
  utm_names = c("X", "Y"),
  utm_crs = get_crs(dat, ll_names),
  units = c("km", "m")
)

Arguments

Value

A copy of dat with the requested UTM coordinate columns.

Examples

d <- data.frame(lat = c(52.1, 53.4), lon = c(-130, -131.4))
add_utm_columns(d, c("lon", "lat"))

Description

Full-rank unstructured $T \times T$ additive-genetic covariance matrix $\boldsymbol G$, parameterised as the Cholesky factor of a rank-$T$ $\boldsymbol\Lambda$ (i.e. the animal_latent() case with $K = T$). Mathematical parallel to phylo_dep().

Usage

animal_dep(formula, pedigree = NULL, A = NULL, Ainv = NULL)

Arguments

Value

See animal_scalar().

See Also

animal_latent(), phylo_dep().

Examples

# Full unstructured additive-genetic trait covariance (rank = n_traits).
# Grounded in test-animal-keyword.R.
ped <- data.frame(
  id   = paste0("i", 1:12),
  sire = c(rep(NA, 4), rep(c("i1", "i2"), length.out = 8)),
  dam  = c(rep(NA, 4), rep(c("i3", "i4"), length.out = 8))
)
A <- pedigree_to_A(ped)
yvec <- as.numeric(MASS::mvrnorm(
  1, mu = rep(0, 2 * 12),
  Sigma = kronecker(diag(2), A) * 0.5 + diag(2 * 12) * 0.5
))
df <- data.frame(
  species = factor(rep(ped$id, each = 2), levels = ped$id),
  trait   = factor(rep(c("t1", "t2"), times = 12), levels = c("t1", "t2")),
  value   = yvec
)
fit <- gllvmTMB(
  value ~ 0 + trait + animal_dep(0 + trait | species, A = A),
  data = df, family = gaussian(), unit = "species", cluster = "species"
)

Description

Per-trait animal-model random intercepts with no cross-trait covariance, using the bar-form syntax; the .indep marker disambiguates printing. Mathematical parallel to phylo_indep().

Usage

animal_indep(formula, pedigree = NULL, A = NULL, Ainv = NULL, common = FALSE)

Arguments

Value

See animal_scalar().

See Also

phylo_indep().

Examples

# Independent per-trait animal-model intercepts via the bar form,
# passing the dense relatedness matrix A directly.
# Grounded in test-animal-keyword.R.
ped <- data.frame(
  id   = paste0("i", 1:12),
  sire = c(rep(NA, 4), rep(c("i1", "i2"), length.out = 8)),
  dam  = c(rep(NA, 4), rep(c("i3", "i4"), length.out = 8))
)
A <- pedigree_to_A(ped)
yvec <- as.numeric(MASS::mvrnorm(
  1, mu = rep(0, 2 * 12),
  Sigma = kronecker(diag(2), A) * 0.5 + diag(2 * 12) * 0.5
))
df <- data.frame(
  species = factor(rep(ped$id, each = 2), levels = ped$id),
  trait   = factor(rep(c("t1", "t2"), times = 12), levels = c("t1", "t2")),
  value   = yvec
)
fit <- gllvmTMB(
  value ~ 0 + trait + animal_indep(0 + trait | species, A = A),
  data = df, family = gaussian(), unit = "species", cluster = "species"
)

Description

Reduced-rank decomposition of the additive-genetic covariance matrix using $K$ latent factors. Set unique = TRUE to add a per-trait diagonal $\boldsymbol\Psi$ companion: $\boldsymbol G = \boldsymbol\Lambda \boldsymbol\Lambda^\top + \boldsymbol\Psi$ with $\boldsymbol\Lambda$ a $T \times K$ loadings matrix ($T$ = number of traits, $K \le T$). The latent factors and diagonal companion both carry the same $\mathbf A$ structure across individuals.

Usage

animal_latent(
  id,
  d = 1,
  pedigree = NULL,
  A = NULL,
  Ainv = NULL,
  unique = FALSE
)

Arguments

Details

This is the canonical "factor-analytic G-matrix" model from quantitative genetics (Kirkpatrick & Meyer 2004; Meyer 2009; the WOMBAT method). Mathematical parallel to phylo_latent() – same engine path with a pedigree-derived relatedness matrix instead of phylogenetic VCV.

Value

See animal_scalar().

See Also

animal_scalar(), animal_indep(), phylo_latent().

Examples

# Factor-analytic G-matrix (d latent factors plus diagonal Psi).
# Grounded in test-animal-keyword.R.
ped <- data.frame(
  id   = paste0("i", 1:12),
  sire = c(rep(NA, 4), rep(c("i1", "i2"), length.out = 8)),
  dam  = c(rep(NA, 4), rep(c("i3", "i4"), length.out = 8))
)
A <- pedigree_to_A(ped)
yvec <- as.numeric(MASS::mvrnorm(
  1, mu = rep(0, 2 * 12),
  Sigma = kronecker(diag(2), A) * 0.5 + diag(2 * 12) * 0.5
))
df <- data.frame(
  species = factor(rep(ped$id, each = 2), levels = ped$id),
  trait   = factor(rep(c("t1", "t2"), times = 12), levels = c("t1", "t2")),
  value   = yvec
)
fit <- gllvmTMB(
  value ~ 0 + trait + animal_latent(species, d = 1, pedigree = ped,
                                   unique = TRUE),
  data = df, family = gaussian(), unit = "species", cluster = "species"
)

Description

Canonical name for the single-scalar (per-trait) animal-model random intercept. Each trait carries an independent length-n_ids draw $\mathbf p_t \sim \mathcal{N}(\mathbf{0}, \sigma^{2}_{\text{a}}\,\mathbf{A})$ where $\mathbf A$ is the additive genetic relatedness matrix (the pedigree-derived numerator-relationship matrix). The variance $\sigma^{2}_{\text{a}}$ is shared across traits.

Usage

animal_scalar(id, pedigree = NULL, A = NULL, Ainv = NULL)

Arguments

Details

Mathematical parallel to phylo_scalar() – same engine path; the only difference is that A is supplied via ⁠pedigree =⁠ (a 3-column data frame: id, sire, dam), ⁠A =⁠ (dense $n \times n$ relatedness matrix), or ⁠Ainv =⁠ (precision matrix). Sparse Ainv inputs use the sparse precision route.

Value

Used inside a gllvmTMB() formula. Returns invisible(NULL) when called outside a formula – the keyword is a syntactic marker that the parser rewrites internally to the canonical phylo_scalar() path (which is family-agnostic at the math level).

See Also

animal_indep(), animal_latent(), animal_dep(), animal_slope(), phylo_scalar(), meta_V() (sampling variance, distinct from relatedness).

Examples

# Pedigree-derived additive-genetic relatedness; single shared
# variance across traits. Grounded in test-animal-keyword.R.
ped <- data.frame(
  id   = paste0("i", 1:12),
  sire = c(rep(NA, 4), rep(c("i1", "i2"), length.out = 8)),
  dam  = c(rep(NA, 4), rep(c("i3", "i4"), length.out = 8))
)
A <- pedigree_to_A(ped)
yvec <- as.numeric(MASS::mvrnorm(
  1, mu = rep(0, 2 * 12),
  Sigma = kronecker(diag(2), A) * 0.5 + diag(2 * 12) * 0.5
))
df <- data.frame(
  species = factor(rep(ped$id, each = 2), levels = ped$id),
  trait   = factor(rep(c("t1", "t2"), times = 12), levels = c("t1", "t2")),
  value   = yvec
)
fit <- gllvmTMB(
  value ~ 0 + trait + animal_scalar(species, pedigree = ped),
  data = df, family = gaussian(), unit = "species", cluster = "species"
)

Description

Adds a per-individual random slope on a continuous covariate x, with individual-level slopes correlated according to the additive- genetic relatedness matrix $\mathbf A$. Mathematically:

$\eta_{io} \;\mathrel{{+}{=}}\; \beta_{\text{a}}(i)\, x_{io}, \qquad \boldsymbol\beta_{\text{a}} \sim \mathcal{N}\bigl(\mathbf 0,\, \sigma^2_{\text{slope}}\,\mathbf A\bigr).$

Usage

animal_slope(formula, pedigree = NULL, A = NULL, Ainv = NULL)

Arguments

Details

This is the canonical quantitative-genetic random-regression (reaction-norm) model – heritable variation in the slope on an environmental gradient. Mathematical parallel to phylo_slope(); the engine recycles the same $\mathbf A^{-1}$ prepared by animal_scalar() / animal_latent().

Value

See animal_scalar().

See Also

animal_scalar(), animal_latent(), phylo_slope().

Examples

# Additive-genetic random regression: one heritable slope variance on a
# continuous covariate x, slopes correlated by the relatedness matrix A.
# Grounded in test-animal-slope-recovery.R.
ped <- data.frame(
  id   = paste0("i", 1:12),
  sire = c(rep(NA, 4), rep(c("i1", "i2"), length.out = 8)),
  dam  = c(rep(NA, 4), rep(c("i3", "i4"), length.out = 8))
)
A <- pedigree_to_A(ped)
df <- expand.grid(
  species = factor(rownames(A), levels = rownames(A)),
  trait   = factor(c("t1", "t2", "t3")),
  rep     = 1:5
)
df$x <- rnorm(nrow(df))
df$value <- rnorm(nrow(df))
fit <- gllvmTMB(
  value ~ 0 + trait + animal_slope(x | species, A = A),
  data = df, unit = "species", cluster = "species"
)

Description

Usage

animal_unique(id, pedigree = NULL, A = NULL, Ainv = NULL)

Arguments

Details

animal_unique() is soft-deprecated as compatibility syntax in gllvmTMB 0.2.0. Use animal_indep() for standalone marginal diagonal animal-model terms. Paired explicit-Psi use remains accepted; use animal_latent(..., unique = TRUE) when the folded latent term itself should carry this diagonal companion.

Canonical name for the D independent animal-model random intercept on a shared relatedness matrix $\mathbf A$. Each trait gets its own variance $\sigma^{2}_{\text{a},t}$; traits share the same per-individual draws (rows of A act on each trait independently). Mathematical parallel to phylo_unique().

Value

See animal_scalar().

Correlated intercept + slope reaction norm

Used with an intercept+slope bar, animal_unique(1 + x | id, pedigree = ped) fits a correlated additive-genetic random-regression (reaction norm) $\mathrm{vec}(\mathbf B) \sim N(0, \Sigma_b \otimes \mathbf A)$, where $\Sigma_b$ is a $2 \times 2$ (intercept, slope) covariance with a free intercept-slope correlation. This routes through the same augmented engine as phylo_unique() with a 1 + x | sp bar; no separate likelihood. For the diagonal ($\rho = 0$) special case use animal_indep() with the same 1 + x | id bar. Recover $\Sigma_b$ with extract_Sigma(fit, level = "phy").

See Also

animal_scalar(), animal_indep(), animal_latent(), animal_dep(), phylo_indep().

Examples

# Per-trait independent additive-genetic variances on a shared A.
# Grounded in test-animal-keyword.R.
ped <- data.frame(
  id   = paste0("i", 1:12),
  sire = c(rep(NA, 4), rep(c("i1", "i2"), length.out = 8)),
  dam  = c(rep(NA, 4), rep(c("i3", "i4"), length.out = 8))
)
A <- pedigree_to_A(ped)
yvec <- as.numeric(MASS::mvrnorm(
  1, mu = rep(0, 2 * 12),
  Sigma = kronecker(diag(2), A) * 0.5 + diag(2 * 12) * 0.5
))
df <- data.frame(
  species = factor(rep(ped$id, each = 2), levels = ped$id),
  trait   = factor(rep(c("t1", "t2"), times = 12), levels = c("t1", "t2")),
  value   = yvec
)
fit <- gllvmTMB(
  value ~ 0 + trait + animal_unique(species, pedigree = ped),
  data = df, family = gaussian(), unit = "species", cluster = "species"
)

Description

Constructs a block-diagonal ⁠n x n⁠ matrix V suitable for the meta_V() / equalto() covstruct in gllvmTMB(), where rows belonging to the same study_id are correlated and rows from different studies are independent. Each block is a compound- symmetric covariance matrix $\mathbf{V}_s = \mathbf{D}_s^{1/2} \mathbf{R}_s \mathbf{D}_s^{1/2}$ with $\mathbf{D}_s = \mathrm{diag}(\sigma^2_{s1}, \ldots)$ the per-row sampling variances within study $s$ and $\mathbf{R}_s$ the within-study correlation matrix (1 on the diagonal, rho_within on the off-diagonal).

Usage

block_V(study_id, sampling_var, rho_within = 0.5)

Arguments

Details

Use this when you have multiple effect sizes per study (multiple outcomes, multiple traits per individual, etc.) and the within-study sampling errors share a common cause. The independent-rows / single- effect-per-study case is block_V(..., rho_within = 0), which equals diag(sampling_var).

Value

A symmetric positive-definite ⁠n x n⁠ matrix.

See Also

meta_V() for the canonical keyword that consumes the matrix produced here.

Examples

set.seed(1)
df <- data.frame(
  study = factor(rep(paste0("s", 1:3), each = 2)),
  var   = c(0.04, 0.05, 0.06, 0.04, 0.05, 0.07)
)
V <- block_V(df$study, df$var, rho_within = 0.5)
round(V, 3)

Description

Use bootstrap_Sigma() when Hessian, Wald, or profile intervals are unavailable or unsafe but the fitted model still has useful point estimates. A pdHess = FALSE or skipped sdreport() is an inference warning, not automatic proof that the fitted mean or rotation-invariant covariance summaries are unusable. Inspect check_gllvmTMB() / gllvmTMB_diagnose() first, then use this helper for parametric simulate-refit uncertainty.

Usage

bootstrap_Sigma(
  fit,
  n_boot = 999,
  level = c("unit", "unit_obs", "phy", "B", "W"),
  what = c("Sigma", "R", "communality", "ICC", "cross_corr"),
  conf = 0.95,
  seed = NULL,
  n_cores = 1,
  progress = TRUE,
  keep_draws = FALSE,
  link_residual = c("auto", "none")
)

Arguments

Details

The function generates n_boot parametric bootstrap replicates of a fitted model and returns percentile confidence intervals for the canonical biological summaries: trait covariance matrices $\hat\Sigma_\mathrm{unit}$, $\hat\Sigma_\mathrm{unit\_obs}$; the corresponding correlation matrices; per-trait communalities $c_t^2 = (\Lambda \Lambda^\top)_{tt} / \Sigma_{tt}$; and per-trait site-level ICCs $R_t = (\Sigma_\mathrm{unit})_{tt} / [(\Sigma_\mathrm{unit})_{tt} + (\Sigma_\mathrm{unit\_obs})_{tt}]$.

Scope: Gaussian bootstrap summaries and mixed-family refit plumbing are covered by current tests. Non-Gaussian bootstrap calibration is experimental: its repeated-sampling behaviour has not been certified, and there is no production calibration evidence for mixed-family intervals. Treat those intervals as indicative.

Each bootstrap replicate (1) draws a new response vector from simulate(fit, nsim = 1), (2) refits the model with the same formula on the simulated data, (3) extracts the requested summaries via extract_Sigma(), extract_communality(), and extract_ICC_site(). Replicates whose refit fails to converge are recorded but excluded from CI calculation.

Parallel refits use future + future.apply and are limited to two workers. Set n_cores = 2 to enable parallel refits. When parallel, replicates use future.apply's L'Ecuyer-CMRG seed stream so the answers are reproducible given a fixed seed, but they are NOT bit-identical to an n_cores = 1 run with the same seed (different RNG streams).

Value

A list with components:

point_est

Named list of point estimates for each requested summary at each requested level (e.g. Sigma_B, R_B, communality_B, ICC_site).

ci_lower, ci_upper

Named lists of percentile CI bounds, element-wise the same shape as the corresponding point_est.

n_effective

Named list (vector summaries only) giving the per-element count of finite bootstrap draws. For ⁠multiple_r_*⁠ entries, any element whose effective count is below 80% of n_boot has its CI bounds set to NA (minimum-effective-B floor).

boot_median

Named list (⁠multiple_r_*⁠ entries only) of the per-element bootstrap median – a cheap sanity comparator against point_est; a large gap flags gross bootstrap corruption and is not a coverage certifier.

ci_method

Character; currently "percentile".

link_residual

Character; the link-residual convention used in point estimates and bootstrap refit summaries.

conf, n_boot, n_failed

Configuration metadata.

draws

NULL unless keep_draws = TRUE; otherwise a named list of bootstrap draw arrays.

Caveats

• Uses the existing simulate.gllvmTMB_multi() method. By default that method REDRAWS the random effects from the fitted covariance, which is what makes these intervals span between-unit variability. CIs reflect parametric simulate-refit variability, not a Bayesian posterior distribution for variance components.

• Intervals are too narrow when the simulator cannot redraw a tier. Redraw is not implemented for every random-effect tier, notably the SPDE spatial tier and the diagonal phylogenetic tier. For a fit using one of those, simulate.gllvmTMB_multi() falls back to reusing the fitted random-effect modes and emits a one-shot warning. That fallback understates between-unit variability, so the intervals returned here are not calibrated for such fits. Treat the warning as a signal that these intervals cannot be trusted. Unsupported families fall back through the simulator's own warning path.

• Refits use the same formula reconstructed from fit$formula and fit$covstructs, and forward the fit's auxiliary structure (phylo_vcv, phylo_tree, mesh, lambda_constraint) so each refit matches the original model.

• Convergence: replicates whose refit fails or whose optimiser does not return convergence == 0 are counted in n_failed and excluded from CIs.

Examples

set.seed(1)
s <- simulate_site_trait(n_sites = 30, n_species = 1, n_traits = 3,
                         mean_species_per_site = 1,
                         Lambda_B = matrix(c(1, .5, -.4), 3, 1),
                         psi_B = c(.2, .15, .1))
fit <- gllvmTMB(value ~ 0 + trait + latent(0 + trait | site, d = 1),
                data  = s$data,
                trait = "trait",
                unit  = "site")
# Fifty refits illustrate the API; they do not establish interval coverage.
boot <- bootstrap_Sigma(fit, n_boot = 50, level = "unit",
                        what = c("Sigma", "R"), seed = 42)
boot$point_est$Sigma_B
boot$ci_lower$Sigma_B
boot$ci_upper$Sigma_B

Description

categorical() declares the predictor-model family for an unordered categorical missing predictor used inside mi() with impute_model(family = categorical()). The missing unordered predictor is marginalised EXACTLY by a finite-state baseline-category softmax sum over its K categories (K >= 3); there is no latent variable. The predictor link is a baseline-category softmax with the FIRST level as baseline (eta_state(1) = 0); for state ⁠k = 2, ..., K⁠, ⁠eta_state(k) = X_x beta_x[block (k-1)]⁠, and ⁠Pr(x = k | z) = exp(eta_state(k)) / sum_j exp(eta_state(j))⁠. The coefficient vector beta_x is packed as K - 1 blocks of the covariate design width. The missing predictor must be an unordered factor, character, or integer category scores with at least three levels, and every category must appear among the observed values; a two-level predictor uses binomial() (the binary route) and an ordered predictor uses cumulative_logit().

Usage

categorical()

Details

This is a predictor-model family tag consumed by impute_model(family = ), NOT a response family for gllvmTMB().

Value

A gllvmTMB_impute_family object for the family argument of impute_model().

See Also

impute_model(), cumulative_logit() for an ordered missing predictor, and binomial() for a binary missing predictor.

Examples

impute_model(habitat ~ z, family = categorical())

Description

Inspects the fit's per-row family vector and flags two configurations that make the link_residual = "auto" path incoherent: within-trait family mixing (errors) and ordinal-probit traits (warns; the latent residual is already standardised at 1, so adding 1 again over-counts).

Usage

check_auto_residual(fit)

Arguments

Value

Invisibly, a list with components status (one of "ok", "warn", or "err" – though the "err" path is unreachable from a successful return because cli::cli_abort halts) and messages (a character vector of human-readable findings; empty when status = "ok").

Side effects

• "err" path: cli::cli_abort is fired and the function does not return.

• "warn" path: cli::cli_warn is fired before the function returns the list invisibly.

• "ok" path: silent.

See Also

extract_correlations(), extract_Sigma(), extract_Omega().

Description

Run check_gllvmTMB() right after fitting, before interpreting confidence intervals or covariance summaries. It returns a stable table of optimiser, gradient, Hessian, sdreport(), restart, boundary, latent-identifiability, and binomial prevalence/loading diagnostics. It is the machine-readable companion to gllvmTMB_diagnose(): use this in simulations, tests, and reports where parsing printed messages would be brittle.

Usage

check_gllvmTMB(
  object,
  gradient_thresh = 0.01,
  se_thresh = 100,
  weak_axis_thresh = 0.05,
  psi_thresh = 1e-04,
  psi_rel_thresh = 0.01,
  sigma_eps_thresh = 1e-04,
  cross_loading_thresh = 0.6,
  binary_prevalence_thresh = 0.9,
  binary_saturation_prob_thresh = 0.99,
  binary_saturation_share_thresh = 0.5,
  loading_relative_thresh = 8,
  loading_runaway_thresh = 25,
  loading_absolute_thresh = 8,
  multinomial_collapse_floor = 1e-10,
  multinomial_collapse_rel_thresh = Inf,
  multinomial_rail_thresh = 0.99,
  multinomial_range_collapse_thresh = 0.02,
  ordinal_loading_runaway_thresh = Inf,
  ordinal_loading_absolute_thresh = Inf
)

Arguments

Details

Scope: optimisation and inference-risk signals for fitted models, including latent-axis rotation, weak-axis, near-zero psi, residual-scale boundary flags, a binomial near-constant/loading/saturation screen, and the intentional gllvmTMBcontrol(se = FALSE) point-estimate path. The table does not calibrate interval coverage, prove formal separation, or prove the selected latent rank by itself. Target-explicit known-DGP simulations will decide when broader interval or rank-selection claims can move beyond diagnostic status.

A WARN row, including pdHess = FALSE, means that Wald standard errors or curvature-based inference need more care; it is not by itself proof that the fitted mean, likelihood, or rotation-invariant covariance summaries are unusable.

Template for an existing fitted model or data set; replace the named inputs with your own objects before running this code.

fit <- gllvmTMB(value ~ 0 + trait + latent(0 + trait | site, d = 2),
                data = dat, trait = "trait", unit = "site")
check_gllvmTMB(fit)

Value

A data frame with columns component, status, value, threshold, message, and action. Status values are "PASS", "WARN", or "FAIL".

Description

Refits the user's data with full unstructured T x T phylogenetic and non-phylogenetic trait covariances (phylo_dep + dep) using the same engine, family, link, and unit / cluster grouping as the supplied two-psi fit. Compares the joint two-psi fit's implied $\boldsymbol\Sigma_{\mathrm{phy}} = \boldsymbol\Lambda_{\mathrm{phy}}\boldsymbol\Lambda_{\mathrm{phy}}^\top + \boldsymbol\Psi_{\mathrm{phy}}$ (and the analogous $\boldsymbol\Sigma_{\mathrm{non}}$) against the unstructured baseline. Per-component RMSE plus a flag field (TRUE when any component disagrees beyond threshold) identifies two-psi identifiability failures.

Usage

compare_dep_vs_two_psi(
  fit_two_psi,
  threshold = 0.1,
  phylo_vcv = NULL,
  phylo_tree = NULL
)

Arguments

Details

Conceptual basis: Williams et al. (2025) bioRxiv 2025.12.20.695312 Eq. 3, generalised across traits. The unstructured fit has $T(T+1)/2$ free parameters per tier and is the strongest available benchmark for the implied total covariance. When the joint two-psi fit is well-identified, its implied total Sigma matches the unstructured baseline; when it isn't, the diagnostic flags the disagreement and the user knows to relax the rank or revisit identifiability.

The cross-check separates three questions:

Level 1: can we fit?

Both fits converged.

Level 2: total Sigma_phy and Sigma_non agree?

This diagnostic answers Level 2 directly. The two estimators target the SAME total covariance with different parameterisations, so agreement is the identifiability check.

Level 3: split into Lambda Lambda^T + Psi?

Diagnostic does not test Level 3 directly; Level 2 agreement is the pre-requisite for Level 3 to be meaningful.

Computational scope: tractable for T <= ~30. For larger T, prefer compare_indep_vs_two_psi().

Template for an existing fitted model or data set; replace the named inputs with your own objects before running this code.

library(ape)
tree <- ape::rcoal(200)
tree$tip.label <- paste0("sp", seq_len(200))
Cphy <- ape::vcv(tree, corr = TRUE)
fit  <- gllvmTMB(
  value ~ 0 + trait + phylo_latent(species, d = 1, unique = TRUE) +
          indep(0 + trait | species),
  data      = df,
  trait     = "trait",
  unit      = "species",
  cluster   = "species",
  phylo_vcv = Cphy
)
diag <- compare_dep_vs_two_psi(fit)
diag$flag
diag$agreement

Value

A list with components:

joint

The two-psi fit's implied Sigma_phy and Sigma_non (T x T matrices).

dep

The unstructured phylo_dep + dep baseline's implied Sigma_phy and Sigma_non (T x T matrices). May contain NULL entries if the alt fit failed.

agreement

Data frame with rows for Sigma_phy and Sigma_non, columns rmse (Frobenius RMSE between joint and dep), dep_mag (Frobenius magnitude of the dep estimate), rel_disagreement (rmse / dep_mag), and flag (rel_disagreement > threshold, or NA when the alt refit failed so the component could not be tested).

flag

Logical: TRUE if any component is flagged, FALSE if all tested components agreed, NA if nothing could be tested (distinguishes agreement from a failed refit).

threshold

The threshold used.

alt_fit

The refit gllvmTMB_multi object (or NULL on failure), retained so users can inspect convergence and pull other extractor outputs.

References

Williams, M. J., McGillycuddy, M., Drobniak, S. M., Bolker, B. M., Warton, D. I., & Nakagawa, S. (2025). Fast phylogenetic generalised linear mixed-effects modelling using the glmmTMB R package. bioRxiv 2025.12.20.695312. doi:10.64898/2025.12.20.695312

Hadfield, J. D. & Nakagawa, S. (2010). General quantitative genetic methods for comparative biology: phylogenies, taxonomies and multi- trait models for continuous and categorical characters. Journal of Evolutionary Biology 23, 494-508. doi:10.1111/j.1420-9101.2009.01915.x

Meyer, K. & Kirkpatrick, M. (2008). Perils of parsimony: properties of reduced-rank estimates of genetic covariance matrices. Genetics 180, 1153-1166. doi:10.1534/genetics.108.090159

Felsenstein, J. (2005). Using the quantitative genetic threshold model for inferences between and within species. Philosophical Transactions of the Royal Society B 360, 1427-1434. doi:10.1098/rstb.2005.1669

Felsenstein, J. (2012). A comparative method for both discrete and continuous characters using the threshold model. American Naturalist 179, 145-156. doi:10.1086/663681

See Also

compare_indep_vs_two_psi() for the cheap diagonal fallback when T is large; extract_Sigma() for the underlying covariance extractor.

Description

Refits the user's data with phylo_indep + indep – a stacked multivariate generalisation of T univariate phylogenetic GLMMs in the spirit of Williams et al. (2025) Eq. 3 – using the same engine, family, link, and grouping as the supplied two-psi fit. Compares the per-trait diagonals of the joint two-psi fit's implied $\boldsymbol\Sigma_{\mathrm{phy}}$ (and $\boldsymbol\Sigma_{\mathrm{non}}$) against the marginal-only baseline. Tests only the diagonals (cheap); does not test the off-diagonal cross-trait covariances.

Usage

compare_indep_vs_two_psi(
  fit_two_psi,
  threshold = 0.1,
  phylo_vcv = NULL,
  phylo_tree = NULL
)

Arguments

Details

Use this diagnostic when T is large (T $\geq$ 30) and the unstructured fit needed by compare_dep_vs_two_psi() is intractable. Otherwise prefer compare_dep_vs_two_psi() (gold standard).

Cross-check intent: the per-trait diagonals are the simplest summary the joint two-psi fit and the marginal baseline both target. If the diagonals already disagree beyond the threshold, no inspection of the cross-trait off-diagonals is warranted.

Template for an existing fitted model or data set; replace the named inputs with your own objects before running this code.

diag <- compare_indep_vs_two_psi(fit)
diag$flag
diag$agreement

Value

A list with components:

joint

Joint two-psi fit's implied diag(Sigma_phy) and diag(Sigma_non) (named numeric vectors of length T).

indep

The marginal-only baseline's diag(Sigma_phy) and diag(Sigma_non) (same shape; may contain NULL if alt fit failed).

agreement

Data frame with rows for Sigma_phy_diag and Sigma_non_diag, columns rmse, indep_mag, rel_disagreement, and flag.

flag, threshold, alt_fit

As in compare_dep_vs_two_psi().

References

Williams, M. J., McGillycuddy, M., Drobniak, S. M., Bolker, B. M., Warton, D. I., & Nakagawa, S. (2025). Fast phylogenetic generalised linear mixed-effects modelling using the glmmTMB R package. bioRxiv 2025.12.20.695312. doi:10.64898/2025.12.20.695312

See Also

compare_dep_vs_two_psi() (canonical, full Sigma); extract_Sigma().

Description

When two latent-variable fits use the same traits, the loading matrices are only identified up to rotation and sign. Procrustes alignment finds the orthogonal transform that brings one matrix as close as possible to the other, then reports the residual disagreement. This is mainly a validation helper for comparing gllvmTMB() with another implementation.

Usage

compare_loadings(Lambda_a, Lambda_b)

Arguments

Details

Template for an existing fitted model or data set; replace the named inputs with your own objects before running this code.

fit_g <- gllvmTMB(value ~ 0 + trait + latent(0 + trait | site, d = 2),
                  data = df, trait = "trait", unit = "site")
fit_t <- glmmTMB::glmmTMB(value ~ 0 + trait + rr(0 + trait | site, d = 2),
                         data = df, REML = FALSE)
L_g <- extract_ordination(fit_g, "unit")$loadings
L_t <- attr(glmmTMB::ranef(fit_t)$cond$site, "loadings")
compare_loadings(L_g, L_t)

Value

A list with the optimal rotation R, the rotated Lambda_a_rot, and the Frobenius distance after alignment.

Description

compare_Sigma_table() joins report-ready extract_Sigma_table() rows to a known covariance or correlation matrix. It is designed for simulation and teaching articles that need estimate-vs-truth tables without hand-indexing matrices inside the article.

Usage

compare_Sigma_table(
  x,
  truth,
  level = "unit",
  part = c("total", "shared", "unique"),
  measure = c("covariance", "correlation"),
  entries = c("unique", "all", "upper", "lower", "offdiag", "diag"),
  link_residual = c("auto", "none")
)

Arguments

Details

Scope: the helper compares fitted or precomputed extract_Sigma_table() rows against one supplied truth matrix (covered by the package's Sigma-table validation tests). It is a table helper only: it does not compute uncertainty, simulate data, or validate calibration. For the first visual comparison layer, use plot_Sigma_comparison(); richer article-specific calibration summaries remain future visualization work.

Value

A data frame with the columns from extract_Sigma_table() plus truth, error, abs_error, and comparison_status.

See Also

extract_Sigma_table(), plot_Sigma_comparison(), plot_Sigma_table().

Examples

rows <- data.frame(
  estimand = "R_unit[length,mass]",
  trait_i = "length",
  trait_j = "mass",
  i = 1L,
  j = 2L,
  level = "unit",
  component = "total",
  matrix = "R",
  estimate = 0.62,
  lower = NA_real_,
  upper = NA_real_,
  interval_method = "none",
  interval_status = "none",
  scale = "correlation",
  diagonal = FALSE,
  triangle = "upper"
)
truth_R <- matrix(c(1, 0.6, 0.6, 1), 2,
  dimnames = list(c("length", "mass"), c("length", "mass"))
)
compare_Sigma_table(rows, truth_R, measure = "correlation")

Description

Use confint_inspect() when a profile confidence interval is missing, one-sided, very wide, or hard to trust. Most users should start with confint.gllvmTMB_multi() or the matching ⁠extract_*()⁠ interval helper; this function is the advanced visual check for a direct profile target.

Usage

confint_inspect(
  fit,
  parm,
  level = 0.95,
  ystep = 0.5,
  ytol = NULL,
  parm.range = c(-Inf, Inf)
)

Arguments

Details

It returns the full profile-likelihood curve, the deviance bounds, and (when ggplot2 is available) a ggplot visualisation showing the curve, the MLE, the chi-squared threshold, and the resulting confidence-interval bounds. It surfaces a small set of numerical shape warnings before the user interprets the bound; the returned label is not a goodness or calibration certificate.

This helper inspects direct profile-ready targets from profile_targets(). Nonlinear derived targets such as communality, repeatability, and trait correlation do not currently have a public profile-inspection route. A working direct curve does not establish empirical interval coverage for every model class.

Common patterns to look for in the returned plot (see the profile-likelihood article):

• Quadratic and approximately symmetric near the MLE: Wald and profile may agree closely for this target.

• Asymmetric or skewed: inspect the endpoint status and constrained-fit stability before preferring either interval.

• Flat near the MLE: the target may be weakly informed, or the numerical profile may need closer inspection.

• Hits a parameter-range bound: interpret the endpoint on the target's transformed scale and distinguish a natural boundary from a truncated search.

• No usable crossing: treat the endpoint as unavailable; do not replace it automatically with bootstrap output unless simulation, refits, failed replicate counts, and Monte Carlo resolution are credible for that target.

Template for an existing fitted model or data set; replace the named inputs with your own objects before running this code.

fit <- gllvmTMB(value ~ 0 + trait +
                latent(0 + trait | site, d = 1),
                data  = sim$data,
                trait = "trait",
                unit  = "site")
inspect <- confint_inspect(fit, parm = "sigma_eps")
inspect$bounds       # the CI + Wald comparison
inspect$diagnostics  # any shape flags
inspect$plot         # ggplot of the profile curve

Value

An object of class gllvmTMB_confint_inspect with components:

⁠$curve⁠

Data frame with one row per profile-grid point. Columns: parm (user-facing label), parm_value_natural (parameter value on the natural scale after the registered transformation), parm_value_link (parameter value on the TMB optimisation scale), nll (negative log-likelihood), deviance_drop (2 * (nll - mle_nll), the chi-squared statistic), excess_over_threshold (deviance_drop - chi2_threshold; negative is inside CI, positive is outside), in_ci (logical).

⁠$bounds⁠

Data frame, 1 row: parm, estimate_natural, lower_natural, upper_natural, wald_lower_natural, wald_upper_natural, wald_profile_disagree_lower, wald_profile_disagree_upper. Wald bounds come from fit$sd_report$cov.fixed via the delta method; disagreement flags are TRUE when the absolute difference exceeds 25% of the Wald half-width (a heuristic for "look at the plot").

⁠$plot⁠

A ggplot object visualising the curve, MLE, chi-squared threshold, profile bounds (vertical lines), and Wald bounds (dashed vertical lines for comparison). NULL if ggplot2 is not installed.

⁠$diagnostics⁠

Named character vector flagging any of: "no_heuristic_warning" (none of this helper's shape warnings fired; this does not prove a quadratic profile), "asymmetric" (skewed, Wald-vs-profile disagrees), "flat_at_mle", "hits_lower_bound", "hits_upper_bound", "no_lower_crossing", "no_upper_crossing", "profile_failed".

⁠$call⁠

The match.call() of the invocation.

See Also

confint.gllvmTMB_multi() (the corresponding CI extractor; same parm vocabulary), profile_targets() (the target inventory), tmbprofile_wrapper() (the lower-level profile API), and the profile-likelihood article.

Description

Returns 95% (or other-level) confidence intervals for fixed effects, variance components, and trait covariance matrices, with three method choices:

Usage

## S3 method for class 'gllvmTMB_multi'
confint(
  object,
  parm,
  level = 0.95,
  method = c("profile", "wald", "bootstrap"),
  nsim = 500L,
  seed = NULL,
  ...
)

Arguments

Details

• method = "profile" (default): profile-likelihood bounds via TMB::tmbprofile() + stats::uniroot(). A successful numerical route can represent asymmetric likelihood shape, but its reference distribution and coverage remain target-specific.

• method = "wald": Gaussian-approximation CIs from sd_report. Fastest; poor near boundaries.

• method = "bootstrap": fitted-model parametric bootstrap via bootstrap_Sigma(). Slowest; inspect failed refits, Monte Carlo resolution, and whether simulation covers every active model component.

Fixed-effect and direct-parameter intervals use the requested method where supported, and Sigma-matrix intervals accept canonical parm = "Sigma_unit" / "Sigma_unit_obs" names, diagonal extra-grouping "Sigma_cluster" / "Sigma_cluster2" names, plus legacy "Sigma_B" / "Sigma_W" aliases. Profile intervals for full decomposed Sigma entries fall back to bootstrap because those entries are nonlinear functions of rotation-equivalent loadings and diagonal $\Psi$. Non-Gaussian bootstrap and broader derived-target coverage are not universally calibrated; inspect the returned method and the target-specific article before reporting bounds. Fits using non-unit likelihood weights target a weighted estimating objective, not the ordinary likelihood used by these interval routes. confint() therefore refuses them; no sandwich interval is currently certified.

Main parm-class dispatch paths:

• Sigma matrices – when parm is one of "Sigma_unit", "Sigma_unit_obs", "Sigma_cluster", "Sigma_cluster2", or "sigma_phy" (legacy aliases "Sigma_B" and "Sigma_W" still work), returns a tidy data.frame with columns parameter, estimate, lower, upper, method. Profile is computed element-wise via TMB::tmbprofile() for the diagonal entries; structural-zero off-diagonal entries in pure-diagonal tiers are labelled "structural_zero"; off-diagonals for reduced-rank tiers fall back to bootstrap (full Sigma sampling) since they mix two parameters in a non-linear way.

• Fixed effects / variance components – when parm is missing, an integer index, or a character vector of fixed-effect term names, returns a numeric matrix with rows = parameters and columns = lower / upper bounds (same shape as stats::confint()). Method choice applies.

• Lambda and derived summaries – "Lambda...", "icc...", "phylo_signal...", "communality...", "rho...", and "proportion..." tokens dispatch to their corresponding direct-profile, Wald, or bootstrap helper where that route is implemented. Nonlinear penalty-profile prototypes are withheld.

Value

• Sigma path – a data.frame with columns parameter (character, e.g. "Sigma_unit[t1,t1]"), estimate (point estimate), lower, upper, and method (the method used for that row, or "structural_zero" for off-diagonal entries that are fixed at zero by a diagonal-only tier). The parameter prefix follows the requested parm, so legacy calls still return legacy "Sigma_B[...]" or "Sigma_W[...]" labels.

• Lambda path – a data.frame with columns parameter (e.g. "Lambda[trait_1,LV1]" for raw rows or "rho[trait_1,LV1]" for standardised rows), estimate, lower, upper, method, loading_scale, pd_hessian, and ci_status. When pdHess = FALSE on a Wald path, lower/upper are NA and ci_status flags the reason).

• Derived-quantity path ("icc" / "phylo_signal" / "communality" / "rho" / "proportion") – a numeric matrix with two columns named after the requested level (e.g. "2.5 %" / "97.5 %") and rownames identifying the entry, e.g. "icc:trait_1", "communality:unit:trait_1", "rho:unit:1,2", or "proportion:shared_unit:trait_1".

• Fixed-effects / variance-component path – a numeric matrix with rownames = parameter names and two columns named "2.5 %" / "97.5 %" (or the analogous quantiles for the requested level).

References

Pawitan, Y. (2001). In All Likelihood: Statistical Modelling and Inference Using Likelihood, Oxford University Press, ch. 9.

Venzon, D. J. & Moolgavkar, S. H. (1988). A method for computing profile-likelihood-based confidence intervals. Applied Statistics 37, 87-94. doi:10.2307/2347496

McCune, K. B., et al. (2024) coxme_icc_ci() – the Nakagawa-authored coxme-based profile-CI helper that inspired this work, in https://github.com/kelseybmccune/Time-to-Event_Repeatability/blob/main/R/rptRsurv.R.

See Also

bootstrap_Sigma(), extract_Sigma(), extract_correlations(), extract_repeatability(), extract_communality(), tmbprofile_wrapper(), loading_ci(), loading_profile(), gllvmTMB().

Examples

## Fit a tiny example
set.seed(1)
s <- simulate_site_trait(
  n_sites = 30, n_species = 4, n_traits = 3,
  mean_species_per_site = 4,
  Lambda_B = matrix(c(0.9, 0.4, -0.3), 3, 1),
  psi_B = c(0.20, 0.15, 0.10),
  beta = matrix(0, 3, 2), seed = 1
)
fit <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | site, d = 1),
  data  = s$data,
  trait = "trait",
  unit  = "site"
)

## Bootstrap CIs for the between-site covariance matrix.
## Fifty refits illustrate the API; they do not establish interval coverage.
## Inspect failed refits and use more replicates for a scientific analysis.
ci_unit <- confint(fit, parm = "Sigma_unit", method = "bootstrap",
                   nsim = 50, seed = 42)
ci_unit

## Wald CIs for fixed effects
confint(fit)

Description

Construct the n_species by d confirmatory loading-constraint matrix M from a discrete functional-group hypothesis. Each named group is pinned to load on a specific latent axis (zero on all other axes), and one anchor species per axis is pinned at +1 to fix the scale. Groups present in group but absent from loads_on remain free on all axes.

Usage

confirmatory_lambda(
  species,
  group,
  d,
  loads_on,
  anchors = NULL,
  axis_labels = NULL
)

Arguments

Details

This is the recommended starting point for confirmatory JSDMs where prior knowledge takes the form "species in group A respond to gradient 1, species in group B respond to gradient 2, ...". For free-form constraint patterns (psychometric CFA, idiosyncratic species-by-species pins), build the matrix directly; for purely statistical identification scaffolding when no biology is in play, see suggest_lambda_constraint().

Value

An length(species) by d numeric matrix with NA entries (estimated), 0 entries (pinned to zero), and 1 entries (pinned at the anchor), ready to pass to gllvmTMB() as lambda_constraint = list(unit = M).

See Also

suggest_lambda_constraint() for the lower-level statistical identification scaffold, gllvmTMB() for the lambda_constraint argument.

Examples

species <- c(paste0("A_", 1:3), paste0("B_", 1:3), paste0("C_", 1:4))
group   <- c(rep("A", 3), rep("B", 3), rep("C", 4))

# Group A loads on axis 1; group B loads on axis 2; group C is free.
M <- confirmatory_lambda(
  species  = species,
  group    = group,
  d        = 2L,
  loads_on = list(A = 1L, B = 2L)
)
M

Description

cumulative_logit() declares the predictor-model family for an ordered categorical missing predictor used inside mi() with impute_model(family = cumulative_logit()). The missing ordered predictor is marginalised EXACTLY by a finite-state cumulative-logit sum over its K categories (K >= 3); there is no latent variable. The predictor link is cumulative logit, with K - 1 free cutpoints ⁠c_1 < ... < c_{K-1}⁠ and a cell probability ⁠Pr(x = k | z) = F(c_k - eta_x) - F(c_{k-1} - eta_x)⁠, where ⁠eta_x = X_x beta_x⁠ (no separate predictor intercept; the cutpoints carry the location). The missing predictor must be an ordered factor (or integer category scores) with at least three levels; a two-level predictor uses binomial() (the binary route).

Usage

cumulative_logit()

Details

This is a predictor-model family tag consumed by impute_model(family = ), NOT a response family for gllvmTMB(). The predictor cumulative-logit link is DISTINCT from the ordinal_probit() RESPONSE family (a probit threshold model).

Value

A gllvmTMB_impute_family object for the family argument of impute_model().

See Also

impute_model(), binomial() for a binary missing predictor, and ordinal_probit() for the ordinal RESPONSE family.

Examples

impute_model(score ~ z, family = cumulative_logit())

Description

Canonical name for the always-alone full-unstructured trait covariance term. Fits a full $T \times T$ $\boldsymbol\Sigma$ with $T(T+1)/2$ free parameters via a Cholesky parameterisation $\boldsymbol\Sigma = \mathbf{L}\mathbf{L}^\top$. Mathematically identical to latent(0 + trait | g, d = T) standalone (where $T$ is the number of traits) – the engine's existing packed-triangular $\boldsymbol\Lambda$ at full rank IS the Cholesky factor of an unstructured $\boldsymbol\Sigma$. The keyword choice is documentary: dep declares user intent that the full unstructured trait covariance is the model.

Usage

dep(formula)

Arguments

Details

This completes the structural-mode quartet:

Decomposition (latent by default)

Shared low-rank loadings plus the default trait-specific Psi companion: $\boldsymbol\Sigma = \boldsymbol\Lambda \boldsymbol\Lambda^\top + \boldsymbol\Psi$. $TK + T$ parameters (rotation-removed).

Marginal (indep standalone)

Per-trait variance with identity correlation: $\boldsymbol\Sigma = \mathrm{diag}(\sigma^2_t)$. $T$ parameters.

Full / unstructured (dep standalone)

Free unstructured $\boldsymbol\Sigma$, PSD via Cholesky. $T(T+1)/2$ parameters.

Use dep() when you want to commit to the full unstructured interpretation explicitly. This is a standard multivariate response covariance model. Computational scope: tractable for $T \le \sim 30$ with the default optimiser; for larger $T$, prefer latent(d = K) with $K \ll T$ as the rank-reduced approximation.

Mutual exclusion

Combining dep(0 + trait | g) with latent(0 + trait | g, d = K) on the same grouping is over-parameterised. Combining dep with unique or indep on the same grouping is redundant (dep standalone already includes the trait-level diagonal). The parser raises cli::cli_abort() in any of these cases.

Template for an existing fitted model or data set; replace the named inputs with your own objects before running this code.

# Full unstructured trait covariance fit:
fit <- gllvmTMB(value ~ 0 + trait + dep(0 + trait | site),
                data = df, trait = "trait", unit = "site")

# The mathematically-equivalent decomposition form (full rank):
fit <- gllvmTMB(value ~ 0 + trait + latent(0 + trait | site, d = T),
                data = df, trait = "trait", unit = "site")

# ERROR: dep + latent on the same grouping is over-parameterised.
# gllvmTMB(value ~ 0 + trait +
#           dep(0 + trait | site) +
#           latent(0 + trait | site, d = 2), data = df)

Value

A formula marker; never evaluated.

See Also

latent(), indep(), phylo_dep(), spatial_dep(), extract_Sigma().

Description

-2 * logLik(object), delegating through logLik.gllvmTMB_multi() rather than recomputing anything independently. This is deliberate: on a penalised (ridged) fit – aghq_ridge set to a finite value – logLik() already discloses, with a warning, that the returned value is the unpenalised log-likelihood evaluated AT a penalised MAP point, not at its own maximum. Because deviance() is defined purely in terms of logLik() here, that same disclosure fires exactly once and is not repeated or re-derived.

Usage

## S3 method for class 'gllvmTMB_multi'
deviance(object, ...)

Arguments

Value

A single numeric value, -2 * as.numeric(logLik(object)). On a ridged fit this carries the same MAP-point caveat as logLik(): the value is -2 times the unpenalised likelihood at the penalised estimate, not at its own maximum, and the warning issued by logLik() is inherited unchanged.

Description

diagnose_kernel_separability() compares two or more dense fixed kernels on the same levels before they are used in a multi-kernel model. It is a pre-fit diagnostic rather than evidence that the components are separately identifiable: when candidate kernels are highly overlapping, component-specific Gamma blocks are weak evidence and should be treated as descriptive unless simulations justify the split.

Use this helper to screen candidate K_phy and K_tip definitions, including raw-network and residualized-network choices, before fitting kernel_latent(..., name = ...) tiers. Note its limits: this is a diagnostic, not recovery evidence, interval calibration, or an in-engine identifiability proof.

Usage

diagnose_kernel_separability(
  ...,
  thresholds = c(near_orthogonal = 0.25, high = 0.7)
)

Arguments

Value

A list of class gllvmTMB_kernel_separability with:

similarity

A symmetric matrix of off-diagonal Frobenius-style similarities between kernels.

pairs

A data frame with pair labels, similarity, overlap class, and a conservative recommendation.

thresholds

The thresholds used for the classes.

note

A short claim-boundary note.

Examples

A_H <- diag(2)
A_P <- diag(2)
rownames(A_H) <- colnames(A_H) <- c("H1", "H2")
rownames(A_P) <- colnames(A_P) <- c("P1", "P2")
W_phy <- matrix(c(1, 0.2, 0.2, 1), 2, 2,
  dimnames = list(rownames(A_H), rownames(A_P)))
W_tip <- matrix(c(0.2, 1, 1, 0.2), 2, 2,
  dimnames = list(rownames(A_H), rownames(A_P)))
K_phy <- make_cross_kernel(A_H, A_P, W_phy, rho = 0.5)
K_tip <- make_cross_kernel(A_H, A_P, W_tip, rho = 0.5)
diagnose_kernel_separability(phy = K_phy, tip = K_tip)

Description

diagnostic_table() turns the metadata attached by predictive_check() and residuals.gllvmTMB_multi() into plain data frames. It is the public path for articles, reports, and tests that need the plotted diagnostic rows, residual rows, row-status counts, or fit-health table without inspecting object attributes directly.

Usage

diagnostic_table(
  x,
  table = c("data", "row_status", "fit_health_status", "check_gllvmTMB")
)

Arguments

Details

Scope: table extraction from existing attr(x, "gllvmTMB_diagnostic") metadata on predictive_check() plots and diagnostic residual data frames. This helper does not compute new diagnostics, refit models, run formal residual tests, or calibrate uncertainty. Richer article tables can build on this stable metadata path once the diagnostic examples are restored.

Value

A plain data frame.

Examples

set.seed(3)
n <- 24
df <- data.frame(
  unit = factor(rep(seq_len(n), each = 2)),
  trait = factor(rep(c("a", "b"), n)),
  value = rpois(2 * n, lambda = 2)
)
fit <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | unit, d = 1),
  data = df,
  trait = "trait",
  unit = "unit",
  family = poisson()
)
p <- predictive_check(fit, type = "rq_qq", seed = 1)
diagnostic_table(p, table = "row_status")
diagnostic_table(p, table = "fit_health_status")
diagnostic_table(p, table = "data")

r <- residuals(fit, type = "randomized_quantile", seed = 1)
diagnostic_table(r, table = "data")

Description

extract_coevolution_modules() standardizes a component-specific cross-lineage covariance block and decomposes it into coupled trait axes. For a named coevolution kernel tier, the helper computes ⁠R = Sigma_row^{-1/2} Gamma Sigma_col^{-1/2}⁠ from the shared covariance returned by extract_Sigma(), then applies a singular-value decomposition.

This is a point-estimate derived-output helper, covered by the package's validation tests, for fitted dense-kernel coevolution models whose component-specific Gamma blocks are already extractable with extract_Gamma(). It does not report uncertainty, choose the biological rank, estimate rho, or calibrate null thresholds. Use simulation, bootstrap, and kernel-separability checks before treating the returned axes as scientific evidence.

Usage

extract_coevolution_modules(
  fit,
  level,
  row_traits,
  col_traits,
  scale = c("shape", "effect"),
  n_modules = NULL,
  tol = sqrt(.Machine$double.eps)
)

Arguments

Details

Template for an existing fitted model or data set; replace the named inputs with your own objects before running this code.

mods <- extract_coevolution_modules(
  fit,
  level = "phy",
  row_traits = c("host_size", "host_defence"),
  col_traits = c("partner_size", "partner_attack")
)
mods$modules
mods$row_axes
mods$col_axes

Value

A list with R, modules, row_axes, and col_axes. R is the standardized cross-lineage correlation-like block. modules has one row per singular axis with its singular value and squared-value share. row_axes and col_axes contain trait loadings for each coupled axis.

Description

$c_t^2 \;=\; \frac{(\boldsymbol\Lambda \boldsymbol\Lambda^{\!\top})_{tt}}{(\boldsymbol\Lambda \boldsymbol\Lambda^{\!\top})_{tt} + \psi_{tt}}.$

Usage

extract_communality(
  fit,
  level = c("unit", "unit_obs", "phy", "B", "W"),
  link_residual = c("auto", "none"),
  ci = FALSE,
  conf_level = 0.95,
  method = c("wald", "bootstrap", "profile"),
  nsim = 500L,
  seed = NULL
)

Arguments

Details

The proportion of trait $t$'s variance that is shared with the other traits via the latent factors. Bounded between 0 and 1. Calls extract_Sigma() internally for the chosen level, so the diagonal uses the full $\boldsymbol\Sigma = \boldsymbol\Lambda \boldsymbol\Lambda^{\!\top} + \boldsymbol\Psi$ decomposition. Ordinary latent() includes $\boldsymbol\Psi$ by default.

Caveat: communality with no-Psi fits

If the fit uses latent(..., unique = FALSE) at the requested level, then $\boldsymbol\Psi = \mathbf 0$ and c_t^2 = 1 for every trait. This is mathematically correct for the no-residual subset but tells you nothing about trait integration. The extract_Sigma() advisory message will fire to flag this. To get meaningful communalities, use ordinary latent() with the default unique = TRUE.

For binomial fits the link-specific implicit residual ($\pi^2/3$ for logit, 1 for probit, $\pi^2/6$ for cloglog) is added to the denominator by default; pass link_residual = "none" to suppress.

Value

When ci = FALSE: a numeric vector indexed by trait. When ci = TRUE: a data frame with columns trait, tier, c2, lower, upper, method. For a bootstrap_Sigma input, the interval columns are copied from the bootstrap object.

See Also

extract_Sigma(); extract_ICC_site(); extract_correlations(); extract_repeatability(); confint.gllvmTMB_multi().

Examples

sim <- simulate_site_trait(
    n_sites = 20, n_species = 6, n_traits = 4,
    mean_species_per_site = 4, seed = 1
  )
  fit <- gllvmTMB(
    value ~ 0 + trait +
            latent(0 + trait | site, d = 2),
    data  = sim$data,
    trait = "trait",
    unit  = "site"
  )
  ## Per-trait between-unit communality.
  extract_communality(fit, level = "unit")
  ## With delta-method Wald CIs.
  extract_communality(fit, level = "unit", ci = TRUE, method = "wald")
  # Fifty refits illustrate the API; they do not establish interval coverage.
  boot <- bootstrap_Sigma(fit, n_boot = 50, level = "unit",
                          what = "communality", progress = FALSE)
  extract_communality(boot, level = "unit", ci = TRUE)

Description

Returns the implied cross-trait correlations from a fit returned by gllvmTMB() at one or more covariance levels. Use the canonical input names "unit", "unit_obs", "phy", and "spatial"; legacy aliases "B", "W", and "spde" still work. The helper returns 95% (or other-level) point estimates by default. Fisher-z and bootstrap intervals can be requested explicitly with the method argument; the former profile token is retained only to return a clear withdrawal error:

Usage

extract_correlations(
  fit,
  tier = "all",
  pair = NULL,
  level = 0.95,
  method = c("none", "fisher-z", "profile", "wald", "bootstrap"),
  n_eff = NULL,
  nsim = 500L,
  seed = NULL,
  link_residual = c("auto", "none")
)

Arguments

Details

• "none" (default): point correlations only. No universal interval calibration has been established across covariance tiers, families, and targets.

• "fisher-z": Fisher's z-transform heuristic interval. Computes $\hat z = \mathrm{atanh}(\hat\rho)$, $\widehat{\mathrm{SE}}(\hat z) = 1/\sqrt{n_{\text{eff}} - 3}$, constructs the CI on z, then back-transforms via $\tanh(\cdot)$ (so bounds are guaranteed inside $[-1, 1]$). Fast (seconds for any T), but the classical iid-correlation variance is not a calibrated mixed-model standard error. Treat these bounds as a sensitivity display and see n_eff.

• "profile": withdrawn. The nonlinear penalty-profile prototype did not yet provide a sufficiently strict constraint and constrained-optimizer diagnostic contract. This token stops with an explanation rather than returning bounds.

• "wald": backward-compat alias of "fisher-z" with the same numerics. Emits a one-shot inform pointing at the canonical name. Kept for scripts that filter on method == "wald".

• "bootstrap": parametric bootstrap via bootstrap_Sigma. Slowest (full sampling distribution); use when a fitted model gives useful point estimates but Hessian- or profile-based intervals are not the right uncertainty summary. Structured tiers not yet resampled by bootstrap_Sigma, currently the SPDE spatial tier, return an explicit Wald/Fisher-z fallback with a message rather than fake bootstrap support.

For T traits at one tier, there are T(T-1)/2 unique correlations. A fit with T = 6 and four covariance levels present has up to 60 cross-trait correlations to report.

Value

A data frame (tibble-like) with columns:

tier

Character level label. The current output stores internal labels "B", "W", "phy", and "spde"; use "unit" and "unit_obs" as input names in new calls.

trait_i, trait_j

Trait names with i < j.

correlation

Point estimate.

lower, upper

Confidence-interval bounds.

method

Method used to compute the CI.

interval_status

Claim-boundary marker: "none" for point-only output, "heuristic_unvalidated" for Fisher-z/Wald bounds, and "target_specific_uncalibrated" for bootstrap bounds. The bootstrap route computes intervals, but this function does not certify their frequentist coverage for the fitted target.

For an engine = "julia" bridge fit, lower/upper are NA, method = "none", and interval_status = "none".

Caveats

• The former nonlinear penalty-profile route is withheld pending an exact constraint solver and explicit constrained-fit diagnostics.

• Bootstrap uses bootstrap_Sigma refits and is the practical fallback when point estimates are useful but Hessian- or profile-based intervals are unavailable. Inspect bootstrap warnings, failed replicates, and interval width before treating the intervals as final.

See Also

extract_Sigma, bootstrap_Sigma, confint.gllvmTMB_multi, extract_communality.

Examples

set.seed(1)
s <- simulate_site_trait(
  n_sites = 80, n_species = 6, n_traits = 4,
  mean_species_per_site = 4,
  Lambda_B = matrix(c(0.9, 0.4, -0.3, 0.5), 4, 1),
  psi_B = c(0.4, 0.3, 0.5, 0.2)
)
fit <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | site, d = 1),
  data  = s$data,
  trait = "trait",
  unit  = "site"
)
## Default: point correlations only.
cors <- extract_correlations(fit, tier = "unit")
## Opt-in Fisher-z sensitivity bounds with an explicitly justified n_eff.
cors2 <- extract_correlations(fit, tier = "unit", method = "fisher-z",
                               n_eff = 60L)
## Fifty refits illustrate the API, rather than establish coverage.
## Inspect failed refits and increase B for a scientific analysis.
cors_b <- extract_correlations(fit, tier = "unit", method = "bootstrap",
                               nsim = 50, seed = 42)

Description

For a fit where a multinomial() trait shares a latent factor with other-family traits, report the association between the nominal trait and each partner trait. A nominal trait spans K-1 baseline-category contrasts, so its association with a single-scale partner is a vector, summarized two ways:

• multiple_r: the reference-invariant multiple correlation $R = \sqrt{\Sigma_{pc}\,\Sigma_{cc}^{-1}\,\Sigma_{cp}/\sigma_{pp}}$ between partner p and the whole K-1 contrast block c. Invariant to the baseline category; magnitude in [0, 1].

• contrast_r (when contrasts = TRUE): the (K-1)-vector of individual contrast correlations, labelled by category-vs-baseline.

Usage

extract_cross_correlations(
  fit,
  level = "unit",
  contrasts = FALSE,
  link_residual = c("auto", "none"),
  method = c("point", "bootstrap", "profile", "wald"),
  conf = 0.95,
  nsim = 500L,
  seed = NULL
)

Arguments

Value

a data.frame, one row per (nominal, partner) pair. With method = "point" the columns are nominal, partner, multiple_r (and contrast_r when contrasts = TRUE). With method = "wald", "bootstrap", per-estimand interval columns are added: multiple_r_lower/multiple_r_upper/multiple_r_method/multiple_r_interval_status (scalar), and when contrasts = TRUE contrast_r_lower/contrast_r_upper (list columns) plus contrast_r_method/contrast_r_interval_status. The estimand NOT served by the requested method keeps point-only interval columns (NA bounds, method = "point"). Computed intervals carry an interval_status flag ("heuristic_unvalidated" for "wald", "target_specific_uncalibrated" for "bootstrap") – coverage is not yet certified for any route (recovery evidence for this route is not yet complete).

Examples

set.seed(12)
n_unit <- 120L
n_rep <- 3L
unit_id <- rep(seq_len(n_unit), each = n_rep)
score <- stats::rnorm(n_unit)
score_obs <- score[unit_id]
gaussian_y <- 0.8 * score_obs + stats::rnorm(length(unit_id), sd = 0.4)
category_y <- vapply(seq_along(unit_id), function(i) {
  eta <- c(0, 0.9 * score_obs[i], -0.6 * score_obs[i])
  prob <- exp(eta - max(eta))
  sample.int(3L, 1L, prob = prob / sum(prob))
}, integer(1L))
dat <- data.frame(
  unit = factor(c(unit_id, unit_id)),
  trait = factor(rep(c("g", "cat"), each = length(unit_id)),
                 levels = c("g", "cat")),
  family = factor(rep(c("g", "m"), each = length(unit_id)),
                  levels = c("g", "m")),
  value = c(gaussian_y, category_y)
)
families <- list(g = gaussian(), m = multinomial())
attr(families, "family_var") <- "family"
fit <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | unit, d = 1, unique = TRUE),
  data = dat, family = families, trait = "trait", unit = "unit"
)
extract_cross_correlations(fit, contrasts = TRUE)

Description

For traits fitted with ordinal_probit(), returns a tidy data frame with the K - 2 estimated cutpoints $\tau_2, \ldots, \tau_{K-1}$ per trait, with optional standard errors from the joint sdreport.

Usage

extract_cutpoints(fit, quiet = FALSE)

Arguments

Details

Convention: gllvmTMB follows Hadfield (2015) — $\tau_1 = 0$ is fixed for identifiability and the K - 2 free cutpoints are reported as cutpoint_2, cutpoint_3, etc. This differs from brms, which reports K - 1 cutpoints as ⁠Intercept[1..K-1]⁠.

This extractor includes Wald standard errors, so it refuses fits made with non-unit likelihood weights. The fitted cutpoint parameters remain point estimates, but weighted-objective uncertainty is not calibrated.

Template for an existing fitted model or data set; replace the named inputs with your own objects before running this code.

fit <- gllvmTMB(
  value ~ 0 + trait + phylo_indep(species, tree = tree),
  data    = ordinal_dat,
  trait   = "trait",
  unit    = "individual",
  cluster = "species",
  family  = ordinal_probit()
)
extract_cutpoints(fit)

Value

A data frame with columns

trait

Trait label (factor level from data[[trait]]).

cutpoint_index

Integer index $k \in \{2, \ldots, K-1\}$.

cutpoint_label

Character label "cutpoint_<k>".

tau_estimate

Estimated $\tau_k$ on the latent (probit) scale.

tau_se

Standard error from the joint sdreport, or NA if the report is unavailable.

If the fit contains no ordinal_probit() traits, returns a zero-row data frame with the same columns.

References

Hadfield, J. D. (2015). Increasing the efficiency of MCMC for hierarchical phylogenetic models of categorical traits using reduced mixed models. Methods Ecol. Evol. 6:706-714. doi:10.1111/2041-210X.12354

See Also

ordinal_probit() for the family constructor and the threshold-trait theory reference list.

Description

extract_Gamma() slices the shared covariance matrix returned by extract_Sigma() to give the row-lineage by column-lineage block ⁠Gamma_shape = Lambda_row Lambda_col^T⁠. In a fitted cross-lineage kernel model, level is the named ⁠kernel_*()⁠ tier, row_traits are the host traits, and col_traits are the partner traits. This is a point-estimate extractor for the shared covariance block of a fitted dense-kernel model, covered by the package's validation tests at the fixed-rho scale: when the fitted kernel tier was built from make_cross_kernel(), scale = "effect" returns Gamma_effect = rho * Gamma_shape using the fixed rho recorded on that kernel. Uncertainty for Gamma is not yet reported by this helper; use bootstrap workflows when interval estimates are required. Use profile_cross_rho() to compare a fixed rho grid; in-engine rho estimation and profile intervals remain planned work.

Usage

extract_Gamma(fit, level, row_traits, col_traits, scale = c("shape", "effect"))

Arguments

Details

rho is part of the supplied K matrix, not a fitted parameter in the current engine. To profile it, use profile_cross_rho() to rebuild K_star over a small rho grid, refit the same formula, and compare logLik() values. scale = "effect" is therefore a fixed-kernel transformation, not an estimate of rho. Treat Gamma from a single association matrix W as data-condition sensitive: the C2 recovery test includes a sparse-versus-dense W check, and sparse or poorly replicated host-partner links should be reported as weaker evidence rather than as a precise coevolution estimate.

For fixed multi-kernel fits, extract_Gamma() also inspects the fitted kernel-similarity diagnostics. If the requested component participates in a high-overlap kernel pair, the returned block is still available, but a warning reminds the caller that component-specific separation is weak evidence.

Template for an existing fitted model or data set; replace the named inputs with your own objects before running this code.

Gamma_HP <- extract_Gamma(
  fit,
  level = "cross",
  row_traits = c("host_size", "host_defence"),
  col_traits = c("partner_size", "partner_attack")
)

Value

A numeric matrix with rows row_traits and columns col_traits.

Description

Returns the trait loading matrix from a fit returned by gllvmTMB(). This is the canonical snake_case accessor in the ⁠extract_*()⁠ family; getLoadings() is an accepted compatibility spelling for readers coming from gllvm. Both return the same matrix.

Usage

extract_loadings(fit, level = "unit", rotate = c("none", "varimax", "promax"))

Arguments

Details

Template for an existing fitted model or data set; replace the named inputs with your own objects before running this code.

extract_loadings(fit, level = "unit", rotate = "varimax")

Value

An ⁠n_traits x d⁠ numeric matrix.

See Also

extract_ordination() for scores and loadings together; getLoadings() for the gllvm-style spelling.

Description

For a predictor-informed ordinary latent fit, latent(..., lv = ~ x) uses the unit-level score model

$\mathbf u_s = \mathbf z_s + \mathbf X_s \boldsymbol\alpha,$

where the innovation $\mathbf z_s$ keeps the usual standard-normal prior. extract_lv_effects() reports either the raw axis-scale $\boldsymbol\alpha$ coefficients or the induced trait-scale contribution $\mathbf B_{\mathrm{lv}} = \boldsymbol\Lambda \boldsymbol\alpha^\top$. The axis-scale table is the default because it matches the usual constrained latent-variable / ordination coefficient. It is conditional on the fitted loading constraint and axis orientation: like the loadings themselves (extract_loadings()), $\boldsymbol\alpha$ is identified only up to the same orthogonal (or oblique) rotation, so raw axis-scale coefficients from two fits – even of the same model on the same data – need not sit on a common axis and are not directly comparable without alignment. This is intrinsic to latent-variable identifiability, not a defect. For cross-fit comparison, prefer type = "trait_effect", the rotation-invariant induced slope surface on the trait linear-predictor scale, described next.

Usage

extract_lv_effects(
  fit,
  level = "unit",
  type = c("axis_effect", "trait_effect"),
  conf.level = 0.95,
  method = c("wald", "profile", "bootstrap"),
  ...
)

Arguments

Details

If the axis-scale representation is specifically required and cross-fit alignment matters, apply the same rotation used for the loadings: the mean latent score is $\mathbf X_s \boldsymbol\alpha$, so it rotates exactly as the scores in rotate_loadings() do – $\boldsymbol\alpha \mathbf T$ for the rotation matrix T returned by rotate_loadings(fit, level, method)$T (verified to reconstruct the rotated mean score to machine precision). There is no separate exported helper that performs this composition; it must be assembled from rotate_loadings()'s T and the estimate column here.

For native TMB fits, std.error is populated from a positive-definite sdreport() when available. Axis-effect SEs come from the fixed-parameter block for alpha_lv_B; trait-effect SEs come from TMB's delta-method ADREPORT(B_lv_unit) output. lower and upper are Wald intervals using conf.level. Broad repeated-sampling coverage has not been established, so finite bounds remain experimental rather than nominally calibrated. For Gaussian, Poisson, NB2, Gamma, Beta, and binomial logit/probit/cloglog engine = "julia" bridge fits, ci_method = "none" exposes point estimates only (std.error, lower, and upper are NA). When the bridge supplies a retained delta-method uncertainty data, extract_lv_effects() surfaces finite std.error, lower, and upper. Those Julia bridge values are not coverage-calibrated intervals.

Value

A data frame. For type = "axis_effect", columns are level, axis, predictor, estimate, std.error, lower, upper, rotation_status and uncertainty_status. For type = "trait_effect", columns are level, trait, predictor, estimate, std.error, lower, upper, and uncertainty_status. The uncertainty_status field states why standard errors are unavailable (for example, skipped sdreport(), a non-positive-definite Hessian, or non-finite standard errors) or labels finite Wald bounds as not yet coverage-calibrated.

See Also

extract_ordination()

Description

Sum of $\boldsymbol\Sigma_\text{tier}$ matrices across the user- selected tiers. The canonical PGLLVM use case is the species-level three-piece decomposition:

$\boldsymbol\Omega \;=\; \boldsymbol\Sigma_\text{phy} \;+\; \boldsymbol\Sigma_\text{non,shared} \;+\; \boldsymbol\Psi_\text{non}.$

Here tiers = c("phy", "B") returns $\boldsymbol\Sigma_\text{phy} + \boldsymbol\Sigma_\text{B}$ where $\boldsymbol\Sigma_B = \boldsymbol\Lambda_B \boldsymbol\Lambda_B^{\!\top} + \boldsymbol\Psi_B$ aggregates the non-phylogenetic shared and unique components.

Usage

extract_Omega(fit, tiers = NULL, link_residual = c("auto", "none"))

Arguments

Details

For two-level behavioural-syndrome fits, tiers = c("B", "W") returns the phenotypic trait covariance $\boldsymbol\Sigma_P = \boldsymbol\Sigma_B + \boldsymbol\Sigma_W$ .

Template for an existing fitted model or data set; replace the named inputs with your own objects before running this code.

# PGLLVM three-piece total
om <- extract_Omega(fit, tiers = c("phy", "B"))
round(om$R_Omega, 2)

# Phenotypic covariance from a two-level behavioural-syndromes fit
om2 <- extract_Omega(fit, tiers = c("B", "W"))

Value

A list with Omega (T × T summed covariance), R_Omega (correlation), tiers_used (which tiers were actually summed), note (notes from each underlying extract_Sigma() call), and (when link_residual = "auto") residual_split — the per-trait $\sigma^2_d / \sigma^2_e / \sigma^2_\text{total}$ data frame from extract_residual_split().

See Also

extract_Sigma(); extract_phylo_signal(); extract_proportions(); extract_residual_split().

Description

Extract ordination scores and loadings from a fitted multivariate model

Usage

extract_ordination(
  fit,
  level = "unit",
  component = c("total", "innovation", "mean")
)

Arguments

Value

A list with scores (units or within-unit observations in rows, latent axes in columns) and loadings (traits in rows, axes in columns).

Examples

sim <- simulate_site_trait(
    n_sites = 20, n_species = 6, n_traits = 4,
    mean_species_per_site = 4, seed = 1
  )
  fit <- gllvmTMB(
    value ~ 0 + trait +
            latent(0 + trait | site, d = 2),
    data  = sim$data,
    trait = "trait",
    unit  = "site"
  )
  ord <- extract_ordination(fit, level = "unit")
  head(ord$scores)
  ord$loadings

Description

For a phylogenetic species-level GLLVM, decomposes each trait's between-species latent variance into three additive components that sum to one:

$H_t^2 \;=\; \frac{[\boldsymbol\Sigma_\text{phy}]_{tt}}{V_{\eta,t}}, \qquad C^2_{\text{non},t} \;=\; \frac{[\boldsymbol\Sigma_\text{non,shared}]_{tt}}{V_{\eta,t}}, \qquad \psi_t \;=\; \frac{[\boldsymbol\Psi_\text{non}]_{tt}}{V_{\eta,t}},$

where $V_{\eta,t} = [\boldsymbol\Sigma_\text{phy}]_{tt} + [\boldsymbol\Sigma_\text{non,shared}]_{tt} + [\boldsymbol\Psi_\text{non}]_{tt}$ is the total between-species latent variance for trait $t$ .

Usage

extract_phylo_signal(
  fit,
  ci = FALSE,
  conf_level = 0.95,
  method = c("profile", "wald", "bootstrap"),
  nsim = 500L,
  seed = NULL,
  link_residual = c("none", "auto")
)

Arguments

Details

Interpretation:

$H_t^2$

phylogenetic signal — proportion of between- species latent variance attributable to phylogenetically structured variation ("evolutionary conservatism"). When the model uses the folded phylo_latent(..., unique = TRUE) decomposition, $\boldsymbol\Sigma_\text{phy}$ is the sum $\boldsymbol\Lambda_\text{phy} \boldsymbol\Lambda_\text{phy}^{\!\top} + \boldsymbol\Psi_\text{phy}$ and $H_t^2$ reflects the total phylogenetic variance.

$C^2_{\text{non},t}$

non-phylogenetic communality — proportion of variance attributable to shared non-phylogenetic axes ("coordinated tip-level lability" across traits).

$\psi_t$

uniqueness — proportion of variance not captured by any shared axis ("relative modularity").

Requires phylo_latent() (optionally with unique = TRUE for a phylogenetic Psi) plus a species-level latent() term, which carries its diagonal Psi companion by default. If the species-tier Psi is absent (a latent(..., unique = FALSE) subset), $\psi_t = 0$ for all traits and a cli::cli_inform() advisory fires.

Template for an existing fitted model or data set; replace the named inputs with your own objects before running this code.

fit <- gllvmTMB(
  value ~ 0 + trait + phylo_latent(species, d = 2) +
                      latent(0 + trait | species, d = 2) +
                      indep(0 + trait | species),
  data       = df,
  trait      = "trait",
  unit       = "species",
  cluster    = "species",
  phylo_tree = tree
)
extract_phylo_signal(fit)

Value

A data frame with columns trait, H2, C2_non, Psi, V_eta (the denominator), one row per trait. The three proportions sum to 1.0 by construction. With link_residual = "auto" an extra link_residual proportion column is inserted before V_eta, the four proportions sum to 1.0, and V_eta includes $\sigma^2_d$. When ci = TRUE, three additional columns are added: H2_lower, H2_upper, H2_method.

Interval calibration

The point estimates are the supported claim. When ci = TRUE, the interval methods are provided for exploration: their empirical coverage is not certified for this estimand, so treat the intervals as exploratory rather than coverage-calibrated. See NEWS.md for the current coverage status.

Which variances enter the denominator

V_eta is the species-level latent variance: only components whose grouping is the cluster column contribute. Concretely $V_\eta = \sigma^2_{phy} + \sigma^2_{non}$, where $\sigma^2_{non}$ collects the species-grouped non-phylogenetic variance.

• When unit == cluster (the usual unit = "species" setup) the ordinary latent() term and its Psi companion are species-level, so both enter.

• When unit != cluster (a crossed ⁠site x species⁠ design) Lambda_B and its Psi are site-level and do not enter; the species-level non-phylogenetic variance is the cluster-tier q_it term (⁠indep(0 + trait | <cluster>)⁠, reported as sd_q^2).

The denominator is defined at the species level. Other defensible definitions exist — for instance a total-variance denominator that also absorbs site-level variation, which would give a smaller H2 for the same fit. Before 2026-07-08 this function used the unit tier unconditionally, which silently reported H2 = 1 for every trait in a crossed design because the q_it variance was never read. Compare extract_proportions() if you want every component reported separately rather than folded into a species-level ratio.

Liability-scale heritability for categorical families

With the default link_residual = "none" the denominator contains only the estimated species-level components, so a phylogenetic-only fit reports H2 = 1 for every trait. For threshold/liability families the conventional phylogenetic heritability (Mizuno et al. 2025 J. Evol. Biol., eq 4/18/19; de Villemereuil & Nakagawa 2014) instead includes the distribution-specific latent residual $\sigma^2_d$ in the denominator. Request it with link_residual = "auto", which adds each trait's fixed link-scale residual (see extract_Sigma()'s "Family-aware link residuals" table) as a fourth component:

• ordinal_probit() (and binomial probit): $\sigma^2_d = 1$ exactly, so a phylogenetic-only fit reports the paper's $H^2 = V_a / (V_a + 1)$ per trait.

• multinomial(): each of the K-1 baseline-category contrasts is a logit, so $\sigma^2_d = \pi^2/3$ per contrast and a phylogenetic-only fit reports $H^2_{(k)} = V_{a(k)} / (V_{a(k)} + \pi^2/3)$ for contrast $k$. Contrast heritabilities are baseline-referenced: the softmax link residual is the full matrix $(\pi^2/6)(\mathbf{I} + \mathbf{J})$ — the shared baseline category couples the contrasts with $\pi^2/6$ off-diagonals — so each $H^2_{(k)}$ is the heritability of the "category k vs baseline" liability under that baseline choice, not of "the trait". Do not average or otherwise collapse the per-contrast values into a single scalar heritability; changing the baseline changes the contrasts.