Package {gllvmTMB}

Type:	Package
Title:	Generalised Linear Latent Variable Models
Version:	0.7.1
Description:	Fits multivariate generalised linear latent variable models using either long-format data or a wide data frame. Models support several response distributions and trait covariance structures, and can estimate correlations, shared latent axes, and trait-specific variances. Binary presence-absence data can be used for joint species distribution analyses. Optional animal, phylogenetic, spatial, and kernel covariance components are available. Estimation uses 'Template Model Builder'. Spatial effects can use sparse Gaussian Markov random field constructions. This 0.7.1 release remains experimental: evidence for estimates is model- and regime-specific. Broad interval coverage is not certified; a 0.94 coverage floor is documented only for one two-sided Gaussian total-variance profile setting.
License:	GPL-3
Language:	en-GB
Copyright:	inst/COPYRIGHTS
URL:	https://itchyshin.github.io/gllvmTMB/, https://github.com/itchyshin/gllvmTMB
BugReports:	https://github.com/itchyshin/gllvmTMB/issues
Depends:	R (≥ 4.1.0)
Imports:	assertthat, cli, fmesher, generics, graphics, lifecycle, Matrix, methods, rlang, stats, tidyselect, TMB (≥ 1.8.0), utils
LinkingTo:	BH, RcppEigen, TMB
Suggests:	ape, DHARMa, detectseparation (≥ 0.4.0), dplyr, future, future.apply, ggforce, ggplot2, galamm, glmmTMB, gllvm, JuliaCall, jsonlite, knitr, lme4, MASS, MCMCglmm, mgcv, mirt, mvtnorm, nadiv, rmarkdown, sf, testthat (≥ 3.0.0), tibble, tidyr, tweedie, vegan, vdiffr, withr
VignetteBuilder:	knitr
ByteCompile:	true
Config/testthat/edition:	3
Config/testthat/parallel:	true
Encoding:	UTF-8
SystemRequirements:	GNU make, C++17
Config/roxygen2/version:	8.0.0
NeedsCompilation:	yes
Packaged:	2026-10-16 10:49:07 UTC; z3437171
Author:	Shinichi Nakagawa [aut, cre, cph]
Maintainer:	Shinichi Nakagawa <itchyshin@gmail.com>
Repository:	CRAN
Date/Publication:	2026-10-28 20:00:08 UTC

gllvmTMB: Generalised Linear Latent Variable Models with TMB

Description

Fits multivariate generalised linear latent variable models using either long-format data or a wide data frame. Models support several response distributions and trait covariance structures, and can estimate correlations, shared latent axes, and trait-specific variances. Binary presence-absence data can be used for joint species distribution analyses. Optional animal, phylogenetic, spatial, and kernel covariance components are available. Estimation uses 'Template Model Builder'. Spatial effects can use sparse Gaussian Markov random field constructions. This 0.7.1 release remains experimental: evidence for estimates is model- and regime-specific. Broad interval coverage is not certified; a 0.94 coverage floor is documented only for one two-sided Gaussian total-variance profile setting.

Current limitations and boundaries

The online current limitations and boundaries page follows the development version and may describe changes made after 0.7.1. For the scope and restrictions of 0.7.1, use this release's NEWS and the help and vignettes installed with the package.

Author(s)

Maintainer: Shinichi Nakagawa itchyshin@gmail.com (ORCID) [copyright holder]

Authors:

• Shinichi Nakagawa itchyshin@gmail.com (ORCID) [copyright holder]

See Also

Useful links:

• https://itchyshin.github.io/gllvmTMB/

• https://github.com/itchyshin/gllvmTMB

• Report bugs at https://github.com/itchyshin/gllvmTMB/issues

Structural AGHQ eligibility gate (internal)

Description

Measures the connected-component structure and an upper bound on the treewidth of each random-effect block in a fitted TMB object's Hessian sparsity graph, and routes each to "quadrature" or "laplace".

Usage

.aghq_gate(obj, data, threshold = 4L)

Arguments

Value

A data.frame with one row per distinct (parameter-name signature, size, treewidth, route) combination found among the random-effect Hessian's connected components: block, size, n_components, treewidth, route ("quadrature" | "laplace"), reason.

Multinomial K-1 contrast pseudo-trait degeneracy screen

Description

multinomial() (family_id 16) is fitted as K-1 baseline-category contrast pseudo-traits per categorical response (expand_multinomial_response(), named "<base>:<category>"). Three failure mechanisms are specific to that expansion and invisible to the generic weak-axis/near-zero-psi rows, which screen every trait's loadings pooled together rather than a single categorical response's own K-1 contrasts against each other: The variance-collapse screen (contrast_variance_collapse) flags one contrast whose loading energy (rowSums(Lambda^2)) collapses to ~0, absolutely or relative to its sibling contrasts. The contrast-rail screen (contrast_rail) flags two contrasts of the same response that load almost perfectly on the same axis (⁠|rho| ~ 1⁠ in the implied contrast-level covariance), a rank collapse that only shows up when the tier carries d >= 2 – at d = 1 every healthy fit reaches rho = +-1 by row proportionality, so d = 1 is exempt by construction. The spatial-range screen (spatial_range_collapse) flags a fitted spatial practical range that (sqrt(8) / kappa) collapses relative to the coordinate domain the multinomial response was fit over. Lambda matrices are read directly from object$report rather than through .gllvmTMB_latent_specs(), which deliberately skips Lambda_phy under use$phylo_dep – irrelevant here, since the implied covariance Lambda %*% t(Lambda) is rotation-invariant, and a phylo_dep rail is precisely a failure this row must be able to see.

Usage

.gllvmTMB_multinomial_degeneracy_row(
  object,
  multinomial_collapse_floor = 1e-10,
  multinomial_collapse_rel_thresh = Inf,
  multinomial_rail_thresh = 0.99,
  multinomial_range_collapse_thresh = 0.02
)

Arguments

Value

A one-row data frame in the check_gllvmTMB() row shape, or NULL when the fit has no multinomial (fid 16) contrast pseudo-traits.

Ordinal-probit loading degeneracy screen

Description

Pure ordinal_probit() fits have no auto-Psi block, so the near-zero-Psi rows in check_gllvmTMB() do not apply. This diagnostic row reports two loading summaries, but both thresholds default to Inf and the screen is disabled by default. A PASS with those defaults does not establish that an ordinal fit is free of degeneracy.

Usage

.gllvmTMB_ordinal_degeneracy_row(
  object,
  ordinal_loading_runaway_thresh = Inf,
  ordinal_loading_absolute_thresh = Inf
)

Arguments

Details

A mechanism probe of 24 degenerate ordinal fits found category-level separation, not link saturation: gll_log_pnorm_diff's cutpoint-underflow condition (both bracketing cutpoints more than 8.2924 from eta on the same side) never fired on any observed row of any degenerate fit (flat-row share exactly 0/24 fits), while dichotomising every degenerate fit's response to binary at the middle cutpoint and refitting as binomial(link = "probit") fired the package's EXISTING binomial_prevalence_loading detector on 24/24 refits. A flat-fit/saturation arm therefore has no empirical basis in this probe. This row is instead modeled directly on .gllvmTMB_binomial_prevalence_loading_row()'s loading arms, because the probe found the pathology concentrated in a single trait's loading column (worked example: one trait's loading 44.2 against a true ⁠max|Lambda| = 4.79⁠ while sibling traits stayed near truth) – consistent with per-trait quasi-separation in these fits.

Two arms, both loading-only. Unlike the binomial row there is no prevalence/saturation conjunct: an ordinal trait has no single Bernoulli "prevalence" to test against. The tested ordinal thresholds did not meet the predeclared sensitivity and false-positive targets, so neither arm is enabled by default. Finite thresholds are user-selected exploratory cutoffs, not calibrated defaults.

• O1 (runaway_loading): relative_loading (a trait's largest loading divided by the typical loading among the OTHER ordinal traits, via ⁠.gllvmTMB_max_loading_by_trait(object, reference_traits = <ordinal trait ids>)⁠ – the family-scoped denominator, so a gaussian/binomial partner trait can neither mask nor manufacture an ordinal runaway) at or above ordinal_loading_runaway_thresh.

• O2 (extreme_magnitude): max_loading_unit (unit tiers ONLY – never the SPDE tier; see .gllvmTMB_max_loading_by_trait()'s own comment on the sqrt(4*pi)*kappa normalisation hazard, measured at 6.5e6 against 66 on the unit tier of the same fit) at or above ordinal_loading_absolute_thresh. This is scale-free for ordinal_probit by the same argument that justifies binomial's absolute arm: the probit-liability residual variance is EXACTLY 1 under the Wright/Falconer/Hadfield threshold convention (R/extract-sigma.R, sigma_d^2 = 1 fixed, no free scale parameter), so a loading IS the trait's latent SD in liability units.

Both thresholds default to Inf (disabled). No tested threshold met the predeclared calibration targets. Setting a finite threshold makes the row report that user-selected cutoff; it does not make the cutoff calibrated.

A third statistic is computed for the diagnostic table but is NOT wired into flag or status: cutpoint_span (a trait's fitted cutpoint span, max(tau) - min(tau) over its free cutpoints ⁠tau_2 .. tau_{K-1}⁠, via .gllvmTMB_ordinal_cutpoint_span_by_trait()) and the derived loading_over_span (max_loading_unit / cutpoint_span). The calibration analysis found cutpoint span associated with the degeneracy label, so this ratio is not treated as an independent scale correction or screening arm. K = 2 traits (no free cutpoint under the Hadfield tau_1 = 0 convention) return NA for cutpoint_span rather than dividing by zero, and loading_over_span is NA wherever cutpoint_span is NA or non-positive.

Value

A one-row data frame in the check_gllvmTMB() row shape, or NULL when the fit has no ordinal_probit() (family_id 14) trait.

Parametric-bootstrap CIs for per-(trait, component) variance proportions

Description

Simulate nsim response vectors via simulate(fit, nsim = nsim), refit the same formula on each, and recompute extract_proportions(). Percentile-quantile bounds across replicates that converged. Reuses the formula-reconstruction and auxiliary-args plumbing from bootstrap_Sigma() so phylogenetic / spatial / lambda-constraint fits refit cleanly.

Usage

.proportions_bootstrap_ci(
  fit,
  components = NULL,
  trait_idx = NULL,
  level = 0.95,
  nsim = 500L,
  seed = NULL
)

Arguments

Value

A data frame with columns trait, component, proportion, lower, upper, method. method = "bootstrap" for percentile rows, "(unavailable)" for link-residual rows.

Wald (delta-method) CIs for per-(trait, component) variance proportions

Description

Numerical-Jacobian Wald CI on the logit-transformed proportion ⁠p_c,t = Sigma_c[tt] / Sigma_total[tt]⁠. The Jacobian is built by perturbing each fixed parameter and re-reading the relevant ⁠Lambda_*⁠ / ⁠sd_*⁠ slots from fit$tmb_obj$report(); the asymptotic covariance is ⁠J Cov(theta_fix) J^T⁠ from fit$sd_report$cov.fixed. Bounds are computed on z = qlogis(p) (symmetric Wald) and back-transformed via plogis(), so they always lie in [0, 1] and respect the bounded-support asymmetry.

Usage

.proportions_wald_ci(fit, components = NULL, trait_idx = NULL, level = 0.95)

Arguments

Value

A data frame with columns trait, component, proportion, lower, upper, method. method = "wald" for delta-method rows, "(unavailable)" for link-residual rows (bounds collapse to point for fixed-scale families; bounds are NA for mean-dependent families).

Response-family constructors used by gllvmTMB

Description

Build family descriptors for likelihoods that are not supplied directly by stats::family(). The descriptor records the family name, link functions, and any fixed shape control consumed by gllvmTMB().

Usage

Beta(link = "logit")

lognormal(link = "log")

gengamma(link = "log")

gamma_mix(link = "log", p_extreme = NULL)

lognormal_mix(link = "log", p_extreme = NULL)

nbinom2_mix(link = "log", p_extreme = NULL)

nbinom2(link = "log")

nbinom1(link = "log")

truncated_poisson(link = "log")

truncated_nbinom2(link = "log")

truncated_nbinom1(link = "log")

student(link = "identity", df = NULL)

tweedie(link = "log", p = NULL)

censored_poisson(link = "log")

delta_gamma(link1, link2 = "log", type = c("standard", "poisson-link"))

delta_gamma_mix(link1 = "logit", link2 = "log", p_extreme = NULL)

delta_gengamma(link1, link2 = "log", type = c("standard", "poisson-link"))

delta_lognormal(link1, link2 = "log", type = c("standard", "poisson-link"))

delta_lognormal_mix(
  link1,
  link2 = "log",
  type = c("standard", "poisson-link"),
  p_extreme = NULL
)

delta_truncated_nbinom2(link1 = "logit", link2 = "log")

delta_truncated_nbinom1(link1 = "logit", link2 = "log")

delta_poisson_link_gamma(link1 = "log", link2 = "log")

delta_poisson_link_lognormal(link1 = "log", link2 = "log")

betabinomial(link = "logit")

delta_beta(link1 = "logit", link2 = "logit")

Arguments

Details

These functions construct descriptors; successful construction does not by itself mean that a route is in the dependable release core. Ordinary Gaussian, Poisson, NB2, and binomial Laplace routes have the strongest first-release evidence. Other families and delta, mixture, truncated, ordinal, or mixed-family routes have narrower evidence. For version 0.7.1, consult its NEWS and installed help for the route-specific evidence limits. The online limitations page is supplementary and may differ as development continues.

Mixed-family input

Supply a list of descriptors when response traits use different families. The data column family selects an element for each row. Set attr(families, "family_var") to select a differently named column.

Ordering contract. Name the list to match the selector column's values, e.g. list(gaussian = gaussian(), student = student()) – matching by name is unambiguous and always accepted. An unnamed list is matched positionally to the selector column's levels (levels() for a factor, alphabetical sort(unique(...)) otherwise); this is validated, not merely assumed, by checking that reading against each family's own name (e.g. a "student" level naturally means student()) agrees with where the list put it. When the two readings agree – which is true whenever the list is already written in level order – the fit proceeds silently. When they disagree, that disagreement is precisely the silent-swap failure mode (issue #1120: list(student(), gaussian()) against family = rep(c("student","gaussian"), each = n) used to pair the student rows with gaussian() and vice versa, with no warning and a converged fit), and the list is refused with an error showing both readings. A selector column whose values carry no name evidence at all (arbitrary labels such as "count"/"binary") still uses list order, but the resolved pairing is reported once so it stays auditable.

Ordinarily the family may vary between traits but not within one: all rows of a given trait must share one family and link, and a fit that mixes them stops with an error. Poisson-log with binomial-logit, for instance, is not a coherent common scale merely because it converges.

Integrated multi-source models (experimental)

To integrate more than two sources – a portal stream, digitised literature records, checklists, a structured survey – declare every source and its observation law with isdm_sources(). The two-source contract below is the two-source case of the same rule and keeps working unchanged.

Integrated two-source models (experimental)

There is one admitted exception, for combining opportunistic presence-only records with a structured detection/non-detection survey of the same species. Here the family genuinely does vary within a trait: the portal rows are Poisson-log counts and the survey rows are Bernoulli with a complementary log-log link, which is what makes both arms consistent with one shared underlying intensity. Poisson and Bernoulli carry no dispersion parameter, so nothing per-trait becomes ambiguous.

To reach it, the fit must match this contract exactly:

• family = list(gbif = poisson(), survey_pa = binomial("cloglog")), with attr(family, "family_var") <- "isdm_family";

• a source column holding only "gbif" and "survey", with both present;

• an isdm_family column equal to "gbif" on portal rows and "survey_pa" on survey rows.

Within that contract an offset() is admitted on the cloglog arm, where it is a known change-of-support (sampled area or visit effort) rather than an ordinary log-rate exposure. Anything short of the full contract keeps the ordinary one-family-per-trait refusal.

Everything such a fit reports is relative intensity. Presence-only data cannot identify absolute abundance, occupancy, or detectability, so those quantities are not estimated and must not be read off the output. This interface is experimental and may change.

Delta defaults and compatibility wrappers

When link1 is omitted, standard delta constructors use "logit" and Poisson-link constructors use "log". The older delta_poisson_link_gamma() and delta_poisson_link_lognormal() names remain as deprecated wrappers.

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

## Mixed-family example: per-row family chosen by a `family` column
sim <- simulate_site_trait(n_sites = 40, n_species = 10, n_traits = 2,
                           mean_species_per_site = 5)
sim$data$family <- factor(
  ifelse(sim$data$trait == levels(sim$data$trait)[1], "gaussian", "binomial"),
  levels = c("gaussian", "binomial")
)
fam <- list(gaussian(), binomial())
fit <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | site, d = 1),
  data   = sim$data,
  family = fam,
  trait  = "trait",
  unit   = "site"
)

## Integrated two-source model: portal counts + survey detections, one
## shared ecological predictor. `dat` is long, one row per
## (cell, species, source); `log_support` is log sampled area or effort.
dat$isdm_family <- factor(
  ifelse(dat$source == "gbif", "gbif", "survey_pa"),
  levels = c("gbif", "survey_pa")
)
isdm_fam <- list(gbif = poisson(), survey_pa = binomial(link = "cloglog"))
attr(isdm_fam, "family_var") <- "isdm_family"
fit_isdm <- gllvmTMB(
  value ~ 0 + trait + trait:env + trait:isdm_gbif + trait:bias +
    offset(log_support) + latent(0 + trait | cell_id, d = 1),
  data   = dat,
  family = isdm_fam,
  trait  = "trait",
  unit   = "cell_id"
)

gengamma() uses the Prentice parameterisation. Its fitted gengamma_Q value controls shape; the limiting case at zero is lognormal.

Several constructors remain exported for compatibility although the current multivariate fitter rejects them. In particular, this applies to generalized-gamma, mixture, censored Poisson, truncated NB1, and most delta variants. A rejected constructor is not silently mapped to another family.

Names ending in ⁠_mix()⁠ describe two mean components that share a scale parameter. p_extreme is the probability assigned to the component with the larger mean; leave it NULL to estimate that probability.

NB2 has variance \mu + \mu^2 / \phi; its extra variance therefore grows quadratically with the mean.

NB1 has variance \mu(1 + \phi), so its extra variance grows linearly with the mean.

student() estimates degrees of freedom unless df is supplied.

Value

A family object. One-part descriptors contain family, link, linkfun, and linkinv, plus constructor-specific metadata where applicable. Delta descriptors contain their two component families and delta = TRUE.

References

Prentice, R.L. 1974. A log gamma model and its maximum likelihood estimation. Biometrika 61(3): 539–544. doi:10.1093/biomet/61.3.539

Thorson, J.T., Stewart, I.J., and Punt, A.E. 2011. Accounting for fish shoals in single- and multi-species survey data using mixture distribution models. Can. J. Fish. Aquat. Sci. 68(9): 1681–1693. doi:10.1139/f2011-086.

Hilbe, J. M. 2011. Negative binomial regression. Cambridge University Press.

Thorson, J.T. 2018. Three problems with the conventional delta-model for biomass sampling data, and a computationally efficient alternative. Canadian Journal of Fisheries and Aquatic Sciences, 75(9), 1369-1382. doi:10.1139/cjfas-2017-0266

Examples

Beta(link = "logit")
lognormal(link = "log")
gamma_mix(link = "log")
lognormal_mix(link = "log")
nbinom2_mix(link = "log")
nbinom2(link = "log")
nbinom1(link = "log")
truncated_poisson(link = "log")
truncated_nbinom2(link = "log")
truncated_nbinom1(link = "log")
student(link = "identity") # estimate df
student(link = "identity", df = 3) # fix df at 3
tweedie(link = "log")
censored_poisson(link = "log")
delta_gamma()
delta_gamma_mix()
delta_gengamma()
delta_lognormal()
delta_lognormal_mix()
delta_truncated_nbinom2()
delta_truncated_nbinom1()
betabinomial(link = "logit")
delta_beta()

Variance partition by source

Description

Decomposes the marginal trait variance into contributions from each active component of the model: between-unit shared (latent_B), between-unit unique (unique_B), within-unit shared (latent_W), within-unit unique (unique_W), phylogenetic (phylo_scalar / phylo_latent), non-phylogenetic species, spatial (spatial), and Gaussian/lognormal observation residual where present. Non-Gaussian link-implicit residual shares are handled by extract_proportions().

Usage

VP(fit)

Arguments

Details

Mirrors gllvm::VP() / gllvm::plotVP().

Value

An ⁠n_traits × n_components⁠ matrix of variance shares (rows sum to 1). Columns are only those active in fit$use.

See Also

extract_proportions() for the canonical per-trait variance-share decomposition with explicit B / W / phy / link-residual columns.

Add UTM coordinates to a data frame

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"))

Full unstructured animal-model trait covariance: animal_dep(0 + trait | id)

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"
)

Independent per-trait animal-model random intercepts: animal_indep(0 + trait | id)

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"
)

Reduced-rank animal-model latent factors: animal_latent(id, d = K)

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"
)

Single-shared-variance animal-model random effect: animal_scalar(id)

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"
)

Animal-model random slope on a continuous covariate

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"
)

Per-trait independent animal-model random intercepts: animal_unique(id)

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"
)

Build a block-diagonal sampling-variance matrix V

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)

Bootstrap covariance, correlation, communality, and ICC summaries

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

Baseline-category softmax unordered categorical missing-predictor family

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())

Check whether link_residual = "auto" is coherent for this fit

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().

Check convergence, Hessian, gradients, and interval readiness

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".

Compare fitted Sigma-table rows with a known truth matrix

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")

Canonical likelihood-based cross-check for the paired phylogenetic decomposition

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.

Cheap diagonal cross-check for the paired phylogenetic decomposition (large T)

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().

Compare two loading matrices after Procrustes alignment

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.

Confidence intervals for a fitted gllvmTMB model

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)

Inspect profile confidence-interval shape for a fitted model

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.

Build a confirmatory lambda_constraint matrix from group membership

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

Cumulative-logit ordered categorical missing-predictor family

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())

Full unstructured trait covariance: dep(0 + trait | g)

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().

Deviance of a fitted gllvmTMB model

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.

Trait-specific unique variance term: unique(0 + trait | group)

Description

Arguments

Details

The unique() formula keyword is soft-deprecated as compatibility syntax in gllvmTMB 0.2.0. For standalone marginal diagonal tiers, use indep() instead; for the scalar standalone marginal case, use indep(..., common = TRUE). Ordinary latent() now carries the diagonal \boldsymbol\Psi companion by default; paired latent() + unique() remains accepted compatibility syntax. Removal is a later API-change decision while the parser and exports remain live. The legacy paired unique(..., common = TRUE) parsimony knob is still accepted as compatibility syntax; new ordinary intercept-only code should use latent(..., common = TRUE).

A formula keyword for adding trait-specific independent random effects to a gllvmTMB() fit, complementing the reduced-rank shared-variance term latent(0 + trait | group, d = K). Using the two together implements the standard psychometric / behavioural-syndrome decomposition

\boldsymbol\Sigma_g = \boldsymbol\Lambda \boldsymbol\Lambda^\top + \boldsymbol\Psi,

where \boldsymbol\Lambda \boldsymbol\Lambda^\top is the low-rank shared component (from latent()) and \boldsymbol\Psi is a diagonal matrix of trait-specific unique variances (now included by ordinary latent() by default, or by legacy explicit unique()). This is the decomposition every published GLLVM treatment uses (Bartholomew et al. 2011; McGillycuddy et al. 2025).

Why ordinary latent() now carries Psi for Gaussian / lognormal / Gamma fits

In earlier development builds, latent(0 + trait | site, d = K) without unique(0 + trait | site) fit only \boldsymbol\Lambda \boldsymbol\Lambda^\top. That constrained each trait's between-site variance to come entirely from the K shared latent factors and inflated correlations when trait-specific residual variance was present. Ordinary latent() now includes the diagonal \boldsymbol\Psi companion by default, so the decomposition is the default instead of a two-keyword spelling.

Set latent(..., unique = FALSE) only when you deliberately want the no-residual subset where communality reaches 1 by construction.

When you do not need unique()

• Binary responses with a probit / logit / cloglog link. The link function fixes an implicit residual variance on the latent scale (1 for probit, \pi^2/3 for logit, \pi^2/6 for cloglog), which acts as the implicit unique component. Adding an explicit unique() term on top is typically not identified.

• Source-specific latent terms. Use phylo_latent(..., unique = TRUE), animal_latent(..., unique = TRUE), spatial_latent(..., unique = TRUE), or kernel_latent(..., unique = TRUE) when the folded term itself should carry its source-specific diagonal \boldsymbol\Psi companion. The explicit compatibility spellings phylo_latent(..., unique = FALSE) + phylo_unique(), animal_latent(..., unique = FALSE) + animal_unique(), spatial_latent(..., unique = FALSE) + spatial_unique(), and kernel_latent(..., unique = FALSE) + kernel_unique() remain accepted. See the Pitfalls article, section 5.

• Confirmatory factor models. Sometimes domain knowledge tells you the latent-only model is correct (no trait-specific residuals); confirm with a likelihood-ratio test.

Two-level (between + within) models

For repeated-measures / behavioural-syndrome data, ordinary latent() now emits the diagonal \boldsymbol\Psi companion by default, so the recommended pattern is two latent() terms:

value ~ 0 + trait +
        latent(0 + trait | individual, d = d_B) +
        latent(0 + trait | obs_id,     d = d_W)

giving \boldsymbol\Sigma_B = \boldsymbol\Lambda_B \boldsymbol\Lambda_B^\top + \boldsymbol\Psi_B and \boldsymbol\Sigma_W = \boldsymbol\Lambda_W \boldsymbol\Lambda_W^\top + \boldsymbol\Psi_W.

For ordinary Gaussian random-regression fits, the same default applies to augmented latent(1 + x | unit, d = K) terms. The model estimates \boldsymbol\Lambda_{\text{aug}}\boldsymbol\Lambda_{\text{aug}}^\top + \boldsymbol\Psi_{\text{aug}} over the ⁠(intercept, slope) x trait⁠ coefficient vector, and extract_Sigma(level = "unit_slope", part = "unique") reports the fitted diagonal \boldsymbol\Psi_{\text{aug}}. Explicit augmented unique(1 + x | unit) remains compatibility syntax and is currently Gaussian-only.

Phylogenetic + non-phylogenetic species-level models

For a species-level fit with phylogeny, the natural three-component decomposition is

\boldsymbol\Omega = \boldsymbol\Sigma_\mathrm{phy} + \boldsymbol\Sigma_\mathrm{non}.

Here \boldsymbol\Sigma_\mathrm{non} is the non-phylogenetic species-level covariance:

\boldsymbol\Sigma_\mathrm{non} = \boldsymbol\Lambda_\mathrm{non}\boldsymbol\Lambda_\mathrm{non}^\top + \boldsymbol\Psi_\mathrm{non}.

In gllvmTMB syntax:

value ~ 0 + trait +
        phylo_latent(species, d = K_phy) +             # Sigma_phy
        latent(0 + trait | species, d = K_non)         # Sigma_non

extract_Sigma() with level = "phy" returns \Sigma_\mathrm{phy}; level = "unit" returns \boldsymbol\Sigma_\mathrm{non} (the non-phylogenetic species-level covariance). Their sum is \boldsymbol\Omega.

Per-row indep() / legacy unique() and sigma_eps: auto-suppression

For Gaussian and lognormal fits, the engine also estimates a single observation-scale residual sigma_eps (the sigma_eps of the response). Ordinary Gamma carries its own per-trait shape phi_gamma and does not use sigma_eps. In new code, write observation-level diagonal residual terms as indep(0 + trait | g). The legacy unique(0 + trait | g) spelling is still accepted as compatibility syntax. If that grouping g has one row per (trait, g) cell (i.e. the diagonal random effects are at the per-row / per-observation level), the diagonal-variance parameters and sigma_eps are jointly unidentifiable – only the sum \mathrm{sd}_g[t]^2 + \sigma_\varepsilon^2 is identified.

In that case the engine auto-suppresses sigma_eps (fixed at \approx 10^{-3} of sd(y)) so the diagonal random effects fully absorb the row-level residual variance, and emits a one-shot message announcing the suppression. This matches the user's intent when they write a per-row indep() term: they want the diagonal variance to represent the row-level residual, not to compete with sigma_eps for it.

If you have multiple rows per (trait, g) cell (e.g. indep(0 + trait | site) with several species per site), sigma_eps is the within-cell residual and the diagonal random effects are the between-cell per-trait variance – both are separately identified and both are estimated.

Family-aware interpretation

The same standalone diagonal tier has a slightly different meaning across families, depending on whether the response carries an observation-layer residual:

For Gaussian / lognormal / Gamma fits, standalone indep() (or legacy standalone unique()) estimates the trait-specific residual variance on the (log-)response scale. For binomial fits with probit / logit / cloglog links, the link function fixes a distribution-specific implicit residual; an explicit diagonal term is identifiable only when there are repeated rows per cell. For Poisson and other log-link families, indep(0 + trait | unit_obs) doubles as an observation-level random effect (OLRE / additive overdispersion).

Across families the unifying rule is: the standalone diagonal tier is the effective per-trait variance parameter at tier g, with the family determining what counts as observation-layer versus latent-scale residual. See the link_residual argument of extract_Sigma() for how the family-specific implicit residual is added to the diagonal of the reported Sigma.

Note on the function name

R's base base::diag() returns the diagonal of a matrix. The formula keyword unique(0 + trait | group) inside a gllvmTMB() formula is recognised by the formula parser at fit time; it is not a call to base::diag(). This help page documents the formula keyword.

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

# Behavioural-syndrome / two-level pattern:
fit <- gllvmTMB(
  value ~ 0 + trait +
          latent(0 + trait | individual, d = 2) +
          latent(0 + trait | obs_id,    d = 1),
  data     = df,
  trait    = "trait",
  unit     = "individual",
  unit_obs = "obs_id"
)
extract_Sigma(fit, level = "unit", part = "shared")$Sigma # Lambda Lambda^T
extract_Sigma(fit, level = "unit", part = "unique")$s     # diag(Psi)
extract_Sigma(fit, level = "unit", part = "total")$Sigma  # both, summed

Value

A formula marker; never evaluated as a call.

References

• Bartholomew, Knott & Moustaki (2011) Latent Variable Models and Factor Analysis: A Unified Approach. Wiley. ISBN 978-0-470-97192-5.

• McGillycuddy, Popovic, Bolker & Warton (2025) Parsimoniously Fitting Large Multivariate Random Effects in glmmTMB. J. Stat. Softw. 112(1). doi:10.18637/jss.v112.i01

• Westneat, Wright & Dingemanse (2015) The biology hidden inside residual within-individual phenotypic variation. Biological Reviews 90: 729–743. doi:10.1111/brv.12131

See Also

extract_Sigma() — pull \boldsymbol\Sigma, \boldsymbol\Lambda \boldsymbol\Lambda^\top, or \boldsymbol\Psi from the fit; phylo_scalar(), phylo_latent(), spatial_indep(), re_int().

Diagnose fixed-kernel separability before fitting

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)

Extract report-ready tables from diagnostic objects

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")

Extract a cross-lineage Gamma block

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.

Site / individual-level ICC per trait

Description

\mathrm{ICC}_t \;=\; \frac{(\boldsymbol\Sigma_B)_{tt}}{(\boldsymbol\Sigma_B)_{tt} + (\boldsymbol\Sigma_W)_{tt}}.

Usage

extract_ICC_site(fit, link_residual = c("auto", "none"))

Arguments

Details

Calls extract_Sigma() internally for both levels with part = "total", so the diagonal of each \boldsymbol\Sigma includes the \boldsymbol\Psi component when it is present. Ordinary latent() includes \boldsymbol\Psi by default; if either level uses latent(..., unique = FALSE), the corresponding advisory message fires and the ICC is computed against the no-Psi diagonal.

For binomial fits the implicit link residual is included in the within-unit variance by default (matching the marginal latent-scale ICC convention); set link_residual = "none" to suppress.

Value

Numeric vector indexed by trait, or NULL if either Sigma_B or Sigma_W is unavailable.

See Also

extract_proportions() for the canonical per-trait variance decomposition (B / W / phy / link-residual shares); extract_Sigma() for the unified covariance API; extract_communality().

Total trait covariance Omega summed across requested tiers

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().

Extract the implied trait covariance / correlation at one tier

Description

Implements the decomposition

\boldsymbol\Sigma_\text{tier} \;=\; \underbrace{\boldsymbol\Lambda_\text{tier}\boldsymbol\Lambda_\text{tier}^\top}_{\text{shared (latent)}} \;+\; \underbrace{\boldsymbol\Psi_\text{tier}}_{\text{unique}},

where ordinary latent() now carries both \boldsymbol\Lambda and the diagonal \boldsymbol\Psi companion by default. This is the same decomposition the behavioural-syndromes / phenotypic-integration literature uses (Bartholomew et al. 2011).

Usage

extract_Sigma(
  fit,
  level = c("unit", "unit_slope", "unit_obs", "phy", "phy_slope", "column_slope",
    "spatial", "spde_slope", "cluster", "cluster2", "B", "B_slope", "W", "spde"),
  part = c("total", "shared", "unique", "psi"),
  link_residual = c("auto", "none"),
  .skip_warn = FALSE
)

Arguments

Details

Legacy covariance wrappers

extract_Sigma_B() and extract_Sigma_W() are soft-deprecated compatibility wrappers. They still return their historical Sigma_B / R_B and Sigma_W / R_W field names, respectively, but new code should call this function with level = "unit" or level = "unit_obs".

When a fit has no Psi component

If the formula deliberately uses latent(0 + trait | unit, d = K, unique = FALSE), the engine fits only the \boldsymbol\Lambda \boldsymbol\Lambda^\top component. Calling extract_Sigma(fit, level, part = "total") then returns just the shared component. This is useful for an intentionally loadings-only target. It should not be described as a full latent-liability covariance unless the relevant family/link residual is also included. Whether the diagonal is incomplete depends on the estimand and response family.

For Gaussian / lognormal / Gamma fits this function emits an advisory note when a reduced-rank tier has no Psi component. Use the ordinary latent(..., unique = TRUE) default for \boldsymbol\Lambda\boldsymbol\Lambda^\top + \boldsymbol\Psi.

For non-Gaussian families (binomial, Poisson, Gamma) the latent-scale residual variance has a closed-form approximation that should be added to the diagonal of \boldsymbol\Sigma – see the link_residual argument below.

The part argument

"total" (default)

\boldsymbol\Sigma_\text{tier} = \boldsymbol\Lambda \boldsymbol\Lambda^\top + \boldsymbol\Psi – the matrix users almost always want for reporting correlations.

"shared"

\boldsymbol\Lambda \boldsymbol\Lambda^\top only – the reduced-rank, rotation-invariant component. Diagonals are \sum_\ell \Lambda_{t\ell}^2; these are not the trait variances, they are the shared part of the trait variances.

"unique"

\boldsymbol\Psi_\text{tier} only – the trait-specific unique variances, returned as a length-T named numeric vector (the diagonal of \boldsymbol\Psi).

Caveat: "shared" vs "unique" partition is only weakly identified

The total \boldsymbol\Sigma_\text{tier} = \boldsymbol\Lambda \boldsymbol\Lambda^\top + \boldsymbol\Psi is rotation-invariant, so part = "total" avoids arbitrary loading rotation. Rotation invariance does not guarantee that the data identify every covariance entry; fit health, sample design, and uncertainty still matter. The split between \boldsymbol\Lambda \boldsymbol\Lambda^\top and \boldsymbol\Psi is only weakly identified – different optimiser starts can flow trait t's variance more into the shared ("shared") or unique ("unique") component, with the same total likelihood. In a synthetic Poisson + OLRE recovery run with target \sigma^2_S = (0.5, 0.4, 0.3, 0.6), gllvmTMB returned (0, 0.67, 0, 0.57) – total trait variances correct, but the partition arbitrary.

This is the standard rotation-and-shift indeterminacy of factor models (gllvm and Hmsc have it too); communality c_t^2 computed from a single fit therefore inherits the same indeterminacy. For interpretation, lean on part = "total" for correlations and use extract_communality() only when paired with bootstrap_Sigma() or the variance-decomposition consistency checks in gllvmTMB_diagnose().

Family-aware link residuals

For non-Gaussian responses each row carries an implicit observation-level residual on the latent (link) scale. Adding it to the per-trait diagonal of \boldsymbol\Sigma defines one family-specific latent-scale convention. It can make the denominator explicit, but it does not make different response families directly commensurate without additional scientific assumptions.

Per-family formulas (Nakagawa & Schielzeth 2010; Nakagawa, Johnson & Schielzeth 2017):

multinomial() is the exception to the scalar-per-trait table: its distribution-specific residual is not a scalar but the (K-1)\times(K-1) matrix (\pi^2/6)(\mathbf{I} + \mathbf{J}) over the category-contrast pseudo-traits – \pi^2/3 on the diagonal (each baseline contrast is a logit, exactly as for binomial(link = "logit")) and \pi^2/6 off-diagonal (the shared baseline category couples the contrasts). This is the residual covariance implied by the softmax's random-utility (Gumbel) representation (McFadden 1974), and it reduces to the binomial \pi^2/3 at K = 2. (It is distinct from MCMCglmm's arbitrary fixed identification residual (1/K)(\mathbf{I}+\mathbf{J}) and that package's MCMC-specific c^2 overdispersion correction, neither of which a direct softmax fit needs.) link_residual = "auto" (the default) applies this full-matrix residual: the (\pi^2/6)(\mathbf{I}+\mathbf{J}) block is added to the multinomial trait's (K-1)-contrast block of \boldsymbol\Sigma (the diagonal \pi^2/3 entry-by-entry with the other families, the \pi^2/6 off-diagonal as a block); "none" returns the latent-scale V unchanged. The addition is tier-agnostic – it applies to whichever tier \boldsymbol\Sigma is requested (level = "phy", a spatial or kernel tier, or an ordinary latent term) – and it makes the categorical block commensurable with any single-scale (e.g. Gaussian) trait sharing the fit. Relatedly, a diagonal among-category V is not independence: the independence null for the contrasts is itself (\mathbf{I}+\mathbf{J})-structured (equal diagonal, equal off-diagonal).

For mixed-family fits the residual is computed per trait when that trait has one unambiguous family/link, then added to the diagonal of \boldsymbol\Sigma entry-by-entry. The default link_residual = "auto" applies this; "none" returns the fitted model covariance without link-residual additions (\boldsymbol\Lambda \boldsymbol\Lambda^\top + \boldsymbol\Psi where a latent decomposition is present). For continuous-only Gaussian or lognormal fits "auto" is a no-op (their per-trait \sigma^2_d is zero).

Citations: Nakagawa & Schielzeth (2010); Nakagawa, Johnson & Schielzeth (2017); McFadden (1974) for the multinomial-logit random-utility residual — see References below.

Future: 3+ latent tiers

The engine currently supports two latent() tiers ("unit" and "unit_obs"); legacy level = "B" and "W" aliases are still accepted with a soft-deprecation message. level = "phy" extracts the phylogenetic implied \boldsymbol\Sigma_\text{phy} from phylo_latent(). If a future release adds 3+ latent tiers, level = "<colname>" will dispatch to the corresponding tier without API change. For now, custom strings error with a clear unsupported-level message.

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 | unit, d = 2),
  data  = df,
  trait = "trait",
  unit  = "unit"
)
extract_Sigma(fit, level = "unit", part = "total")$Sigma   # full T x T cov
extract_Sigma(fit, level = "unit", part = "shared")$Sigma  # rr-only
extract_Sigma(fit, level = "unit", part = "unique")$s      # diag(s_unit)

Value

For part = "total" or "shared": a list with components Sigma (T x T matrix), R (T x T correlation matrix; only for "total"), level, part, and note (character vector of advisory messages, e.g. about a no-Psi unique = FALSE fit).

For part = "unique": a list with s (length-T named numeric vector of unique variances), level, part, note.

For a phylo_dep(1 + x1 + ... + xs | species) fit with one or more slopes, call with level = "phy": the result is the single full unstructured ⁠(1+s)T x (1+s)T⁠ covariance over the trait-stacked (intercept, slope_1, ..., slope_s) random-effect columns – a list with Sigma and R carrying INTERLEAVED dimnames (per trait: ⁠intercept.<t>⁠, then ⁠slope.<t>⁠ for s == 1 or ⁠slope.<x_j>.<t>⁠ for s >= 2), level = "phy_dep", part = "dep", and a note. The part and link_residual arguments do not apply to this single unstructured block and are ignored. (The unit / unit_obs tiers return NULL for a dep-only fit, as it carries no between/within-unit covariance term.)

Current phylo_indep(1 + x | species) uses the same interleaved ⁠2T⁠ report channel with all cross-trait blocks fixed to zero. Its return has level = "phy_indep_slope" and part = "indep"; each trait retains its own intercept-slope block. This is distinct from the legacy shared ⁠2 x 2⁠ augmented phylo_unique channel.

For augmented spatial slopes, call with level = "spatial". spatial_unique(1 + x | coords) retains the shared ⁠2 x 2⁠ field-scale channel and returns level = "spde_base_slope". Current spatial_indep(1 + x | coords) returns the interleaved, cross-trait block-diagonal ⁠2T x 2T⁠ field covariance with level = "spde_indep_slope" and part = "indep"; spatial_dep returns the full unstructured block with level = "spde_dep" and part = "dep". These SPDE matrices are on the fitted field-covariance scale; the returned note records the marginal-scale conversion.

For an ordinary individual-level random-regression fit with latent(1 + x | unit, d = K), indep(1 + x | unit), or their long-form equivalents, call with level = "unit_slope": the result is the augmented ⁠2T x 2T⁠ covariance over trait-specific intercept and slope coefficients, with row names ⁠intercept.<trait>⁠ and ⁠slope.<x>.<trait>⁠. As for level = "unit", part = "shared" returns ⁠Lambda_aug Lambda_aug^T⁠, part = "unique" returns the augmented diagonal, and part = "total" returns their sum.

References

Nakagawa, S. & Schielzeth, H. (2010). Repeatability for Gaussian and non-Gaussian data: a practical guide for biologists. Biological Reviews 85, 935-956. doi:10.1111/j.1469-185X.2010.00141.x

Nakagawa, S., Johnson, P. C. D., & Schielzeth, H. (2017). The coefficient of determination R^2 and intra-class correlation coefficient from generalized linear mixed-effects models revisited and expanded. Journal of the Royal Society Interface 14(134), 20170213. doi:10.1098/rsif.2017.0213

McFadden, D. (1974). Conditional logit analysis of qualitative choice behavior. In P. Zarembka (Ed.), Frontiers in Econometrics (pp. 105-142). Academic Press, New York. (Random-utility origin of the multinomial-logit (\pi^2/6)(\mathbf{I}+\mathbf{J}) contrast residual.)

See Also

extract_communality() for the per-trait shared / unique variance share at one tier; extract_proportions() for the canonical per-trait variance decomposition across tiers; extract_Omega() for the multi-tier sum.

Between-site covariance matrix Sigma_B (soft-deprecated compatibility wrapper)

Description

Returns the implied between-unit trait covariance \Sigma_B = \Lambda_B \Lambda_B^\top + \boldsymbol\Psi_B and its correlation matrix. This is a thin wrapper around the unified extract_Sigma() function: equivalent to extract_Sigma(fit, level = "unit", part = "total") but with the historical Sigma_B / R_B field names.

Usage

extract_Sigma_B(fit)

Arguments

Details

Prefer extract_Sigma() for new code; that interface lets you ask for the rr-only "shared" component or the diag-only "unique" component, exposes the binomial-link implicit-residual option, and will be the entry point for 3+ rr tiers when the engine adds them.

Value

A list with Sigma_B (T x T covariance), R_B (correlation), or NULL if no rr/diag term is present at the between-unit tier.

Lifecycle

Soft-deprecated in 0.7.0. Existing calls continue to return the historical Sigma_B / R_B names, but new code should use extract_Sigma(fit, level = "unit").

See Also

extract_Sigma() — the canonical unified API; pass level = "unit" for the same matrix.

Within-site covariance matrix Sigma_W (soft-deprecated compatibility wrapper)

Description

Returns \Sigma_W = \Lambda_W \Lambda_W^\top + \boldsymbol\Psi_W and the correlation. Thin wrapper around extract_Sigma() with level = "unit_obs". Prefer the unified interface for new code.

Usage

extract_Sigma_W(fit)

Arguments

Value

A list with Sigma_W and R_W, or NULL.

Lifecycle

Soft-deprecated in 0.7.0. Existing calls continue to return the historical Sigma_W / R_W names, but new code should use extract_Sigma(fit, level = "unit_obs").

See Also

extract_Sigma() — the canonical unified API; pass level = "unit_obs" for the same matrix.

Extract a report-ready table of Sigma or correlation entries

Description

extract_Sigma() returns matrices because that is the most convenient interactive representation. extract_Sigma_table() returns the same covariance target as one row per matrix entry so articles, tables, and plot helpers can work without hand-indexing matrices. It can also turn a bootstrap_Sigma() result into the same row schema with bootstrap interval columns filled in.

Usage

extract_Sigma_table(
  fit,
  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

The helper reshapes the same covariance and correlation matrices returned by extract_Sigma(). A bootstrap_Sigma() object can supply interval columns. Julia bridge fits currently expose only point estimates for the ordinary unit tier. The helper does not create profile or Wald intervals.

The table is a point-estimate view over extract_Sigma(). It does not compute confidence intervals from a fitted model directly. Use extract_correlations() when you need pairwise correlation intervals, or bootstrap_Sigma() followed by extract_Sigma_table() when you need bootstrap intervals for Sigma or correlation matrix entries.

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 | unit, d = 2),
  data  = df,
  trait = "trait",
  unit  = "unit"
)
extract_Sigma_table(fit, level = "unit")
extract_Sigma_table(fit, level = "unit", measure = "correlation")
extract_Sigma_table(fit, level = c("unit", "unit_obs"), entries = "all")
boot <- bootstrap_Sigma(fit, n_boot = 50, level = "unit",
                        what = "Sigma", progress = FALSE)
extract_Sigma_table(boot, level = "unit", entries = "upper")

Value

A data frame with one row per requested entry and stable columns: estimand, trait_i, trait_j, integer indices i and j, level, component, matrix, estimate, lower, upper, interval_method, interval_status, scale, diagonal, and triangle. Interval columns are NA with interval_method = "none" because this helper is point-estimate only.

See Also

extract_Sigma() for the underlying matrix extractor; extract_correlations() for pairwise correlation intervals; bootstrap_Sigma() for bootstrap uncertainty.

Extract cross-lineage coevolutionary modules

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.

Communality of each trait

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)

Extract cross-trait correlations

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)

Cross-family correlations between a nominal (multinomial) trait and partners

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)

Extract ordinal-probit cutpoints from a fitted gllvmTMB model

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.

Extract the trait loading matrix

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.

Predictor effects on latent-score means

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()

Extract ordination scores and loadings from a fitted multivariate model

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

Phylogenetic-signal proportions per trait

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.

When a categorical family (fid 14/16) is present and link_residual = "none", an advisory fires pointing at link_residual = "auto"; the default is unchanged for backward compatibility. ci = TRUE is not yet implemented for the liability-scale denominator and refuses with a typed error.

See Also

extract_Sigma(); extract_Omega(); extract_proportions(); extract_repeatability(); extract_communality(); extract_correlations().

Per-trait proportion-of-variance decomposition across all model components

Description

For each trait, returns the proportion of total latent variance attributable to each component present in the model: shared (rr) and unique (diag) at each tier (unit, unit_obs, phy), diagonal cluster / cluster2 components when present, plus optionally the binomial link's implicit residual.

Usage

extract_proportions(
  fit,
  link_residual = c("auto", "none"),
  format = c("long", "wide")
)

Arguments

Details

This is the most general proportion-decomposition function in the package. extract_phylo_signal() is the focused PGLLVM H^2 / C^2_\text{non} / \psi^2 convenience wrapper for the species-level case; extract_communality() gives the rr-only "shared" proportion at one tier; extract_ICC_site() gives the between-vs-within proportion. extract_proportions() returns all of these in one tidy frame.

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

extract_proportions(fit)                           # long format
extract_proportions(fit, format = "wide")          # one row per trait

Value

Long format: data frame with columns trait, component (e.g. "shared_unit", "unique_unit", "unique_cluster", "shared_phy", "link_residual"), variance (the absolute variance), proportion (the share of the total). Wide format: data frame with one column per component plus a total_variance column; the per-trait proportions sum to 1.

OLRE interpretation: for fits with a genuine observation-level diagonal term, usually written as per-row indep() in new code (see extract_residual_split()), the historical unique_W component in this output corresponds to \sigma^2_e (the estimated OLRE variance) and the link_residual component corresponds to \sigma^2_d (the distribution-specific latent residual).

See Also

extract_phylo_signal() — the PGLLVM-specific shortcut; extract_communality(); extract_ICC_site(); extract_residual_split() — explicit \sigma^2_d / \sigma^2_e decomposition for OLRE fits.

Per-trait repeatability with confidence intervals

Description

Returns the per-trait repeatability

R_t = v_{B,t} / (v_{B,t} + v_{W,t}),

where

v_{B,t} = [\Lambda_B\Lambda_B^\top]_{tt} + \sigma^2_{B,t}

and

v_{W,t} = [\Lambda_W\Lambda_W^\top]_{tt} + \sigma^2_{W,t} + \sigma^2_{d,t}.

The first two terms include the shared latent and diagonal companion variance at the unit and observation tiers; \sigma^2_{d,t} is the family-specific link residual used for non-Gaussian traits. The function is intended for a fit returned by gllvmTMB() with the relevant ordinary unit and observation-level components. Also known as the intraclass correlation coefficient (ICC) at the unit level.

Usage

extract_repeatability(
  fit,
  level = 0.95,
  method = c("wald", "profile", "bootstrap"),
  nsim = 500L,
  seed = NULL
)

Arguments

Details

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) +
          latent(0 + trait | site_species, d = 1),
  data     = df,
  trait    = "trait",
  unit     = "site",
  unit_obs = "site_species"
)
extract_repeatability(fit)
boot <- bootstrap_Sigma(fit, n_boot = 50, level = c("unit", "unit_obs"),
                        what = "ICC", progress = FALSE)
extract_repeatability(boot)

Value

A data frame with columns trait, R (point estimate), lower, upper, method.

Interval calibration

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

Method choice

• "wald" (default): Gaussian-approximation CI via the delta method on log(v_B) - log(v_W), transformed with plogis().

• "profile": withdrawn. It aborts rather than silently substituting a different interval method or estimand.

• "bootstrap": parametric bootstrap via bootstrap_Sigma().

References

Nakagawa, S. & Schielzeth, H. (2010) Repeatability for Gaussian and non-Gaussian data: a practical guide for biologists. Biological Reviews 85, 935-956. doi:10.1111/j.1469-185X.2010.00141.x

See Also

extract_communality, extract_correlations, extract_phylo_signal, confint.gllvmTMB_multi.

Implied trait covariance or correlation

Description

Canonical snake_case accessors in the ⁠extract_*()⁠ family for the implied trait covariance and correlation at level = "unit" or "unit_obs". getResidualCov() and getResidualCor() are accepted gllvm-compatibility spellings; each pair returns the same matrix.

Usage

extract_residual_cov(fit, level = "unit")

extract_residual_cor(fit, level = "unit")

Arguments

Value

An ⁠n_traits x n_traits⁠ matrix.

See Also

extract_Sigma() for the unified Sigma API; getResidualCov() for the gllvm-style spellings.

Separate OLRE residual variance from the distribution-specific latent residual

Description

For an additive overdispersion (OLRE) model

\eta_{it} = \mathbf{X}\boldsymbol\beta + \ldots + e_{it}, \quad e_{it} \sim N(0, \sigma^2_e),

the total latent-scale residual variance for trait t is

\sigma^2_{d,t} + \sigma^2_{e,t},

where \sigma^2_d is the distribution-specific (theoretical) component that depends only on the family/link, and \sigma^2_e is the estimated OLRE variance – the per-trait diagonal of the within-unit residual covariance \boldsymbol{\Psi}_W.

Usage

extract_residual_split(fit)

Arguments

Details

The function detects whether the fit includes a genuine observation-level random effect: normally a per-row ⁠indep(0 + trait | <obs-level>)⁠ term. When every (trait, obs) cell has one row per observation level per trait, sigma2_e is populated; otherwise it is zero.

Terminology note

Nakagawa & Schielzeth (2010) use \sigma^2_d for both components. Nakagawa, Johnson & Schielzeth (2017) §7 refine the terminology: \sigma^2_d (distribution-specific) applies only to binomial-type families whose link function introduces a fixed latent-scale variance; \sigma^2_\varepsilon (observation-level) applies to overdispersed Poisson / NB / Gamma and is estimated from the data. gllvmTMB keeps the colloquial sigma2_d column name for compatibility but documents the distinction here (NJS 2017 §7).

Per-family \sigma^2_d table

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

## Add a site_species column (one level per row) as the obs-level grouping.
df$site_species <- factor(seq_len(nrow(df)))
fit <- gllvmTMB(
  value ~ 0 + trait + indep(0 + trait | site_species),
  data     = df,
  trait    = "trait",
  unit     = "site",
  unit_obs = "site_species",
  family   = poisson()
)
extract_residual_split(fit)

Value

A data frame with one row per trait and columns:

trait

Factor of trait names.

sigma2_d

Theoretical / parameter-dependent distribution-specific latent residual (computed by the internal link_residual_per_trait() helper; see the per-family table above; zero for gaussian and lognormal).

sigma2_e

Estimated OLRE variance per trait – the per-trait diagonal of \boldsymbol{\Psi}_W when the fit has a genuine observation-level diagonal term, else 0.

sigma2_total

sigma2_d + sigma2_e.

References

Nakagawa, S. & Schielzeth, H. (2010) Repeatability for Gaussian and non-Gaussian data: a practical guide for biologists. Biological Reviews 85(4): 935-956. doi:10.1111/j.1469-185X.2010.00141.x

Nakagawa, S., Johnson, P. C. D. & Schielzeth, H. (2017) The coefficient of determination R^2 and intra-class correlation coefficient from generalized linear mixed-effects models revisited and expanded. Journal of the Royal Society Interface 14(134): 20170213. doi:10.1098/rsif.2017.0213

See Also

extract_Omega() (returns residual_split as a list component when link_residual = "auto"); extract_Sigma(); extract_proportions() (the unique_W component in its output is \sigma^2_e and the link_residual component is \sigma^2_d for OLRE-style fits).

Extract a tidy rotated loading table

Description

Returns one row per trait and latent axis from a fitted gllvmTMB() model. This is the report-ready companion to rotate_loadings() and to plot(type = "ordination"): the same rotation, axis-ordering, sign-anchor, and optional loading-standardization conventions are recorded in columns instead of left inside a matrix.

Usage

extract_rotated_loadings_table(
  fit,
  level = "unit",
  method = c("varimax", "promax", "none"),
  order_axes = TRUE,
  sign_anchor = c("auto", "none"),
  anchor_traits = NULL,
  loading_scale = c("raw", "standardized")
)

Arguments

Details

Scope: covered for fitted latent() components at level = "unit" or level = "unit_obs" with raw or standardized point-estimate loadings. This table does not attach uncertainty intervals to loadings. Bootstrap or simulation-based uncertainty is not yet implemented and remains possible future work for a later inference slice.

Value

A data frame with columns level, trait, axis, loading, abs_loading, axis_variance, axis_share, rotation, order_axes, sign_anchor, anchor_trait, and loading_scale. axis_variance and axis_share are computed from the rotated raw loading matrix before any optional standardization; these are the quantities used to order axes.

See Also

rotate_loadings(), extract_ordination(), extract_Sigma(), plot.gllvmTMB_multi().

Examples

sim <- simulate_site_trait(
  n_sites = 60, n_species = 12, n_traits = 4,
  Lambda_B = matrix(c(1.0, 0.7, -0.3, 0.5,
                      0.3, -0.5, 0.8, 0.2),
                    nrow = 4, ncol = 2),
  psi_B = c(0.3, 0.3, 0.3, 0.3),
  seed = 1
)
fit <- gllvmTMB(value ~ 0 + trait +
                        latent(0 + trait | site, d = 2),
                data  = sim$data,
                trait = "trait",
                unit  = "site")
extract_rotated_loadings_table(
  fit,
  level = "unit",
  method = "varimax",
  loading_scale = "standardized"
)

Fitted values from a fitted gllvmTMB model

Description

A thin wrapper over predict.gllvmTMB_multi() at the training rows (newdata = NULL). Returns the same long data frame predict() returns – one row per training observation, carrying its unit/species/trait identifiers alongside the fitted value – rather than a bare vector, so row identity against the original data is preserved. ... is forwarded to predict() unchanged, so its other arguments (se.fit, re_form) work here too.

Usage

## S3 method for class 'gllvmTMB_multi'
fitted(object, type = c("response", "link"), ...)

Arguments

Value

A data frame with the original row identifiers (unit, species, trait columns) plus an est column on the requested scale – identical in shape to predict(object, newdata = NULL, type = type, ...). For a multinomial() fit this is .predict_multinomial()'s output: K rows per observation (one per category, including the baseline), not one.

Flag loadings whose CI overlaps a "biologically negligible" region

Description

A simple decision aid in the spirit of Zientek & Thompson (2007, doi:10.3758/BF03193163): given per-entry CIs on Lambda, flag the entries whose conf_level confidence interval does not exclude a user-supplied "null region" — a band around zero considered biologically negligible on the interval's recorded loading_scale. Loadings flagged as unreliable = TRUE are the ones for which the data do not provide evidence that the species responds non-trivially to the axis.

Usage

flag_unreliable_loadings(
  fit,
  null_region = c(-0.1, 0.1),
  level = c("unit", "unit_obs"),
  method = "wald",
  conf_level = 0.95,
  loading_scale = NULL
)

Arguments

Details

Pinned entries (set explicitly by lambda_constraint) are reported with unreliable = NA because no inference is being made about them.

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

## Grounded in tests/testthat/test-loading-ci.R. flag_unreliable_loadings()
## accepts either a fit or a loading_ci() data frame directly.
fit <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | site, d = 2L),
  data = df, family = stats::binomial(link = "probit")
)
## From a fit (computes per-entry Lambda CIs internally):
flag_unreliable_loadings(fit, null_region = c(-0.1, 0.1))
## Or pass a precomputed loading_ci() data frame:
ci <- loading_ci(fit, level = "unit")
flag_unreliable_loadings(ci, null_region = c(-0.1, 0.1))

Value

The loading_ci() data frame with one extra logical column unreliable: TRUE if the CI overlaps null_region, FALSE if it lies entirely outside, NA for pinned entries.

See Also

loading_ci().

Extract latent-variable scores from a fitted multivariate model

Description

Returns latent scores from a fit returned by gllvmTMB(). This is a small compatibility wrapper around extract_ordination() for readers familiar with gllvm::getLV().

Usage

getLV(fit, level = "unit", rotate = c("none", "varimax", "promax"), se = FALSE)

Arguments

Details

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

getLV(fit, level = "unit")
getLV(fit, level = "unit", se = TRUE)

Value

When se = FALSE (default): a matrix with one row per unit (level = "unit") or one row per within-unit observation (level = "unit_obs"), and one column per latent factor. When se = TRUE: a list with scores (that same matrix) and se (a matrix of identical shape and dimnames holding the score uncertainty described below).

Score uncertainty

For an ordinary Laplace fit, se = TRUE reads the marginal standard error of every unit-level (or within-unit) latent-score random effect – z_B at level = "unit", z_W at level = "unit_obs" – from the fit's TMB sdreport() (sqrt(sd_report$diag.cov.random)), and reshapes it with the identical matrix(..., nrow = d, ncol = n) then transpose convention that extract_ordination() uses for the point estimates, so scores[i, k] and se[i, k] always refer to the same (unit, axis) cell. This value is mathematically equivalent to inverting the fit's full joint precision matrix (TMB::sdreport(getJointPrecision = TRUE)) and reading the diagonal of the same block; the two routes were verified to agree to machine precision in development checks. Requirements:

• The fit must carry a valid sdreport() (gllvmTMBcontrol(se = TRUE), the default) with a positive-definite Hessian; otherwise se = TRUE raises an error (no sdreport) or a warning with NA standard errors (non-positive-definite Hessian).

• rotate must be "none": rotating scores changes their covariance, which is not currently propagated, so se = TRUE together with rotate != "none" raises an error rather than silently pairing rotated scores with un-rotated standard errors.

• Predictor-informed latent(..., lv = ~ x) fits at level = "unit" are not yet supported (the score mean's own uncertainty is not yet propagated) and raise an error.

• engine = "julia" bridge fits are not yet supported and raise an error.

For a gllvmTMB_va (variational) fit, se is instead the per-unit variational posterior SD read from the fit's own variational distribution at its optimum – not a Wald standard error, and not calibrated (the returned matrix carries a "uncertainty_basis" and "calibrated" attribute making this explicit in the object itself, not only here – though those attributes are silently dropped by se[i, ], head(se) and as.data.frame(se)).

Two independent gates restrict it. By tier, an explicitly requested "jj" (Jaakkola-Jordan) fit is refused pending its own measurement. By mechanism, any fit whose per-unit SD turns out to be constant across units is refused outright: the array would carry one row per unit while containing no per-unit information. That degeneracy is provable under "ac" (Albert-Chib) and was measured on a Gaussian fit (coefficient of variation 1.6e-15) which the tier gate alone did not catch.

See Also

extract_ordination() for scores and loadings together.

Extract the loading matrix from a fitted multivariate model

Description

Returns the trait loading matrix from a fit returned by gllvmTMB(). This is a small compatibility wrapper around extract_ordination() for readers familiar with gllvm::getLoadings(). The canonical snake_case spelling in this package is extract_loadings(), which forwards here.

Usage

getLoadings(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.

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

Value

An ⁠n_traits × d⁠ numeric matrix.

See Also

extract_ordination() for the row-and-column interface that returns scores and loadings together.

Standard errors of random-effect blocks without a scores accessor

Description

getLV() exposes standard errors for the ordinary latent-variable score blocks (z_B / z_W). Several other random-effect families fitted by gllvmTMB() have the identical machinery available – a marginal variance for every random-effect coordinate, already computed by the fit's TMB sdreport() – but no accessor exposed it. getREsd() closes that gap for the blocks whose point-estimate reshape convention is already established elsewhere in the package (see block below).

Usage

getREsd(
  fit,
  block = c("diag_unit", "diag_unit_obs", "diag_species", "phylo", "re_int", "equalto")
)

Arguments

Details

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

getREsd(fit, block = "re_int")
getREsd(fit, block = "diag_unit")

Value

As described per block above. Standard errors are always NA (with a warning) when the fit's Hessian is not positive-definite.

What this is, and is not

Every number returned here is a Wald / sdreport() marginal standard error, propagated through the delta method the same way getLV()'s se = TRUE already is – not a resampled, profiled, or simulation-based quantity. Coverage of intervals built from it has not been measured for these blocks; do not describe or advertise it as exact, and do not construct a CI from it and call the CI "certified" without separately measuring its coverage. diag.cov.random reflects fixed-effect-uncertainty -propagated marginal variance for each coordinate individually; it does NOT give cross-block or random-vs-fixed covariance (that needs sdreport(getJointPrecision = TRUE), which this fit's single production sdreport() call does not request).

Because of that fixed-effect-uncertainty propagation, a near-saturated or otherwise degenerate design (e.g. an indep() diagonal term at a grouping where every group has one observation, with its residual variance driven to the boundary) can return SEs one or more orders of magnitude larger than the block's own conditional (GMRF) variance would suggest – in that regime the random effect deterministically tracks the fixed-effect residual, so its delta-method SE inherits the fixed effect's own uncertainty almost entirely. This is correct Wald behaviour, not a bug, but it means the magnitude is not always "the block's own noise" in an intuitive sense; treat a surprisingly large SE as a cue to check whether the corresponding variance component sits at or near a boundary.

Blocks not covered here (s_B_slope, s_W_slope, r_c2 / diagonal cluster2, g_phy, g_phy_diag, the SPDE spatial fields ⁠omega_spde*⁠, the dense-kernel blocks ⁠g_kernel*⁠, every augmented random-slope block, and the missing-predictor latent blocks x_mis / u_mi_group / g_x) are, in principle, reachable the same way – the underlying sdreport() call already covers their marginal variances too – but either their point-estimate reshape convention is not yet established anywhere else in the package (so getting the row/column order right cannot be cross-checked against existing code), or their dimensions are not carried as top-level fields on the fit object. Adding them safely needs that groundwork first; this function raises a clear error naming the block rather than silently mis-reshaping.

See Also

getLV() for the ordinary latent-score standard errors this function generalises.

Extract implied trait covariance or correlation

Description

Returns the implied trait covariance at level = "unit" or level = "unit_obs". For a reduced-rank plus unique tier this is \Sigma_X = \Lambda_X \Lambda_X^\top + \Psi_X. getResidualCor() returns the corresponding correlation matrix. The canonical snake_case spellings in this package are extract_residual_cov() and extract_residual_cor(), which forward here.

Usage

getResidualCov(fit, level = "unit")

getResidualCor(fit, level = "unit")

Arguments

Value

An ⁠n_traits × n_traits⁠ matrix.

See Also

extract_Sigma() — the canonical unified API for between-/within-/phylogenetic Sigma at any tier.

Infer a UTM coordinate reference system

Description

Infer a UTM coordinate reference system

Usage

get_crs(dat, ll_names = c("longitude", "latitude"))

Arguments

Value

An EPSG code for the predominant UTM zone.

Fit multivariate response models from wide or long data

Description

Top-level entry point for models with several responses per site, individual, species, or study. Start from a wide data frame with traits() on the formula left-hand side, or from already-stacked long data with one row per ⁠(unit, trait)⁠ observation. Both routes estimate the same trait covariance, pairwise correlations, shared latent axes, and trait-specific variance. The formula syntax also supports fixed effects plus covariance-structure keywords organised by correlation source (none / animal / phylo / spatial / kernel) and three modes (independent / dependent / latent). The common = TRUE modifier on ⁠*_indep()⁠ gives the one-shared-variance special case:

Usage

gllvmTMB(
  formula,
  data,
  trait = "trait",
  unit = "site",
  unit_obs = NULL,
  cluster = NULL,
  cluster2 = NULL,
  family = gaussian(),
  weights = NULL,
  REML = FALSE,
  mesh = NULL,
  phylo_vcv = NULL,
  phylo_tree = NULL,
  known_V = NULL,
  lambda_constraint = NULL,
  Xcoef_fixed = NULL,
  control = gllvmTMBcontrol(),
  missing = miss_control(),
  impute = NULL,
  silent = TRUE,
  engine = c("tmb", "julia"),
  ci_method = c("none", "wald", "profile", "bootstrap"),
  ci_level = 0.95,
  ci_nboot = 200L,
  ci_seed = 0L,
  site = NULL,
  species = NULL,
  estimator = c("ml", "mspl")
)

Arguments

Details

The grid lists available syntax. Reliability depends on the family and uncertainty method. Gaussian indep() and dep() point fits are the narrow tested starting routes; structured-source and non-Gaussian routes need their named evidence boundary. General interval coverage is not claimed.

The three covariance modes (indep / dep / latent) encode covstruct intent across traits:

• latent: the decomposition mode \boldsymbol\Sigma = \boldsymbol\Lambda \boldsymbol\Lambda^\top + \boldsymbol\Psi: a low-rank shared part plus a diagonal \boldsymbol\Psi companion. Ordinary latent() carries \boldsymbol\Psi by default; the ⁠unique =⁠ argument controls it (latent(..., unique = FALSE) for the loadings-only subset, ⁠*_latent(..., unique = TRUE)⁠ to fold the \boldsymbol\Psi companion into a source-specific term).

• indep: the marginal-only mode: each trait gets its own variance, no cross-trait covariance.

• dep: the full unstructured mode: \boldsymbol\Sigma is free with T(T+1)/2 parameters via a Cholesky factor.

Plus the supporting phylo_slope() / animal_slope() (random slopes), meta_V() (known-V meta-analytic; meta_known_V() is a deprecated alias), and the engine-internal propto() / equalto() covstructs (used by the canonical keywords above; not typically called directly). See the formula keyword grid article for a one-paragraph tour of each cell.

gllvmTMB() parses the glmmTMB-style formula, converts wide traits() input to the same internal stacked-trait representation as explicit long data, and dispatches to the underlying TMB template. Covariance-structure terms (latent(), indep(), propto(), equalto(), spatial()) are processed by extending the formula parser and the TMB template.

Per the manuscript, when stacking traits one should set dispformula = ~ 0 so that no implicit residual variance competes with structured ⁠indep(0 + trait | …)⁠ terms. gllvmTMB() enforces this internally.

Multi-trial binomial. The TMB engine evaluates dbinom(y, n_trials, p) so binomial fits are not restricted to Bernoulli (size = 1). To pass a per-row trial count, use either:

• cbind(successes, failures) ~ ... on the LHS of the formula (canonical R / glmmTMB convention; recommended), or

• a flat response y with weights = n_trials (the alternative glmmTMB API).

Both interfaces produce identical fits. For Bernoulli data, omit weights and use a 0/1 response; the engine sets n_trials = 1 and the likelihood is identical to the previous Bernoulli-only behaviour.

Missing responses and predictors. NA response cells are treated as unobserved unit-trait cells. Under the default miss_control(response = "drop"), they are dropped before the TMB likelihood is built; with miss_control(response = "include"), they are kept in the model data, masked out of the likelihood, and available through predict_missing(). For cbind(successes, failures) binomial responses, a row is treated as missing when either response component is missing. Observation weights, when supplied, are subset to retained likelihood rows before validation. Ordinary missing predictors, grouping variables, or fixed-effect design values still error because the model cannot construct a design row; explicitly modelled missing predictors use mi(x), miss_control(predictor = "model"), and impute = list(...). The modelled-mi() route is native Laplace only; the variational route refuses mi() predictors. Dense missing-response masks have their own narrower VA fence and do not imply predictor-model support. These contracts are covered by the package's missing-data validation tests, both for responses and for predictors.

Delta (hurdle) families. delta_lognormal() and delta_gamma() use a single linear predictor for both components: presence is \Pr(y > 0) = \mathrm{invlogit}(\eta) and the positive-component mean is E[y \mid y > 0] = \exp(\eta). This matches sdmTMB's standard delta default (one fixed-effects matrix, one set of trait random effects). The shared-predictor scheme is the simplest hurdle formulation; a future release may decouple the two predictors so presence and abundance can have independent fixed and random effects. Only type = "standard" (the default) is currently wired in the multivariate engine; type = "poisson-link" is not implemented. Each delta family carries one per-trait dispersion of the positive component (no extra Bernoulli dispersion). The response must be non-negative.

Per-trait residual variance: when does it activate?

"Residual" in a mixed-effects model is scale-relative: what counts as residual variance shifts as you add levels to the model. In a Gaussian fit without a per-row diagonal term, the residual is row-level noise captured by a single shared sigma_eps. Once you add a per-row indep(0 + trait | obs) term, the row-level residual is now T per-trait random-effect variances and sigma_eps is auto-suppressed to avoid double-counting. If you also add a site-level diagonal term on top, the row-level term remains the residual and the site-level term is now an additional, higher-level random effect, not a residual at all. The dispatch table below records this explicitly for the configurations the engine supports today.

Gaussian and lognormal responses have one residual scale parameter, sigma_eps. Ordinary Gamma responses instead carry a per-trait shape phi_gamma (CV = 1 / sqrt(phi_gamma)). Per-trait residual variances only appear if you explicitly add a per-row indep(...) term. The dispatch is automatic:

No per-row indep, Gaussian/lognormal present

One shared sigma_eps across Gaussian/lognormal rows. Not per-trait.

No per-row indep, no Gaussian/lognormal rows

sigma_eps is mapped off; the family's intrinsic dispersion handles the residual (for ordinary Gamma, phi_gamma).

indep(0 + trait | g) where g has fewer levels than rows (e.g. g = "site")

sigma_eps is still estimated as the row-level residual; the diagonal term adds a per-trait random effect at level g on top.

indep(0 + trait | obs) at the per-row level (one level per row), Gaussian/lognormal rows fitted

sigma_eps is auto-suppressed (mapped off, fixed at a tiny stabiliser); the T per-trait diagonal random effects are the residual. A one-shot cli::cli_inform fires at fit time announcing the auto-suppression.

indep(0 + trait | obs) at the per-row level, non-Gaussian or mixed-family fit

Treated as observation-level random effects (OLRE). For Bernoulli traits the OLRE is statistically unidentifiable and is mapped off; for hurdle / delta families a warning is emitted (see "Per-family-aware OLRE selection" below).

Mnemonic: Gaussian/lognormal sigma_eps is the default; ordinary Gamma uses phi_gamma; per-row indep replaces the Gaussian/lognormal scalar residual with T per-trait residuals; non-per-row indep adds a higher-level random effect on top of the family residual; non-continuous families never carry sigma_eps regardless.

Per-family-aware OLRE selection. When ⁠indep(0 + trait | <unit_obs>)⁠ is at per-row resolution, i.e. one row per ⁠(trait, unit_obs)⁠ cell, the resulting per-trait random effects on the linear predictor are an observation-level random effect (OLRE). The engine now decides per trait what to do with the OLRE variance based on the trait's response family:

• single-trial Bernoulli (binomial(), all rows have n_trials == 1): theta_diag_W[t] and the corresponding s_W column are mapped off. The OLRE is statistically unidentifiable here (no within-cell replication) and the MLE is the trivial boundary \sigma_W \to 0; pinning the parameter at \sigma_W \approx 10^{-6} removes the spurious free parameter and emits a one-shot informational message. Multi-trial binomial (via cbind(succ, fail) or weights = n_trials) is identifiable and fit normally.

• delta / hurdle families (delta_lognormal(), delta_gamma()): OLRE is fit as before, but a one-shot warning is emitted because the OLRE acts on the shared linear predictor of the hurdle and mixes presence and positive-component noise.

• all other families: OLRE is fit as before.

In a mixed-family fit (per-row family list) the per-trait decision applies only to the affected trait. References: Nakagawa & Schielzeth (2010) Biol. Rev. 85: 935-956; Nakagawa, Johnson & Schielzeth (2017) J. R. Soc. Interface 14: 20170213.

Value

A gllvmTMB object. With no covariance-structure terms in the formula the result has class "gllvmTMB" (single-response engine); with latent() or other covstruct sugar it has class c("gllvmTMB_multi", "gllvmTMB") (multi-trait engine). S3 methods are registered on the subclasses, not on bare gllvmTMB: for a multi-trait fit (gllvmTMB_multi), tidy(), predict(), summary(), confint() and logLik() are available, and AIC() / BIC() dispatch through logLik() only when the fitted point is an unpenalised native-Laplace ML maximum. AGHQ uses a different integration objective and a loading ridge returns a penalised MAP point, so ordinary likelihood comparisons do not apply. An MSPL result has leading class "gllvmTMB_mspl"; likelihood comparison, standard errors, intervals, profiles, and hypothesis-test methods fail closed because this experimental surface is point estimation only. Trait-level extractors such as extract_ICC_site() and extract_communality() are multi-trait only.

vcov() and coef() are available too; see gllvmTMB_multi-vcov. coef() works on any fit; vcov() needs the sdreport() and raises the same typed errors confint() does when it is missing or non-finite.

History, 2026-08-04: this block used to say vcov() dispatched on gllvmTMB, which was never true of any release. It and coef() existed only for gllvmTMB_va, where they refuse. The wording was corrected first, then the two methods were added, which is why the promise now holds.

Dispersion parameters in the report slot

The fitted object's ⁠$report⁠ carries per-family dispersion quantities on the natural scale, each a vector of length n_traits indexed by trait (except sigma_eps, a single scalar). Their names follow the engine's internal parameterisation, not standard R distribution arguments; several are easy to misread. What each one IS, and how it maps to standard R arguments (mu = exp(eta) unless noted):

Objective provenance

A fitted object carries two different surfaces that both look like "the objective", and they answer different questions.

report$joint_nll_unpenalized and report$joint_nll_penalized are conditional joint quantities computed on the TMB tape (the negative log-likelihood of the data and random effects at their current values, before and after TMB-side MSPL penalty terms are added). They see only penalties applied inside the C++ template – the LA-MSPL penalties documented above – and never see anything computed at the R level.

fit$objective_components is the marginal surface: the value TMB's Laplace approximation integrates the random effects out of, decomposed into likelihood_nll (the unpenalised marginal negative log-likelihood), ridge_penalty (the R-level loading ridge from aghq_ridge, added outside the TMB objective – see .gllvmTMB_penalised_gradient()), optimization_nll (their sum, the objective the optimiser actually minimised), and optimizer_reported (the value nlminb() / optim() returned). This is the only place the R-level loading ridge appears; ⁠report$joint_nll_*⁠ cannot see it because it never reaches the TMB tape.

The two surfaces are not expected to match, on a ridged fit or otherwise, and a difference between them is not itself an error: they are conditional-joint and marginal quantities from different penalty regimes.

References

Sterzinger, P. and Kosmidis, I. (2023). Maximum softly-penalized likelihood for mixed effects logistic regression. Statistics and Computing, 33, 53. doi:10.1007/s11222-023-10217-3. This paper is the logistic mixed-effects antecedent for the fixed-effect Jeffreys component. The probit, complementary-log-log, reduced-rank GLLVM, and spatial penalties implemented here are gllvmTMB-specific extensions; that paper does not validate those extensions.

See Also

screen_gllvmTMB() for the opt-in fixed-design separation certificate to use before fitting; estimator = "mspl" is documented in this topic. traits() for wide data-frame formula input; gllvmTMB_wide() for wide matrix/data-frame input; simulate_site_trait() for generating recovery test data; extract_Sigma() for the unified post-fit covariance API; gllvmTMB_diagnose() for a one-stop convergence + identifiability health check; ordinal_probit() for the gllvmTMB-native ordinal threshold family.

Examples

set.seed(1)
sim <- simulate_site_trait(n_sites = 50, n_species = 8, n_traits = 3,
                           mean_species_per_site = 5)
# Long format: one response value per site-species-trait row.
fit_long <- gllvmTMB(
  value ~ 0 + trait + (0 + trait):env_1 + (0 + trait):env_2,
  data   = sim$data,
  family = gaussian(),
  trait  = "trait",
  unit   = "site_species"
)

# Wide format: one row per site-species, with the same three traits.
# `traits()` expands `1 + env_1 + env_2` to the explicit long terms above.
wide_values <- stats::reshape(
  sim$data[c("site_species", "trait", "value")],
  idvar = "site_species", timevar = "trait", direction = "wide"
)
wide <- merge(
  unique(sim$data[c("site_species", "env_1", "env_2")]), wide_values,
  by = "site_species", sort = FALSE
)
names(wide) <- sub("^value[.]", "", names(wide))
fit_wide <- gllvmTMB(
  traits(trait_1, trait_2, trait_3) ~ 1 + env_1 + env_2,
  data   = wide,
  family = gaussian(),
  unit   = "site_species"
)
summary(fit_long)
summary(fit_wide)

Run an advanced Laplace-consistency check

Description

gllvmTMB_check_consistency() is an advanced validation helper for the Laplace approximation, not a first-use diagnostic. Start with check_gllvmTMB() or gllvmTMB_diagnose(). Use this function when a fit looks numerically acceptable but you need a direct check of whether the approximate marginal score is centred.

Usage

gllvmTMB_check_consistency(fit, n_sim = 100L, seed = NULL, estimate = FALSE)

Arguments

Details

It simulates n_sim datasets from the fitted model under the joint parameter vector at the MLE, then tests whether the approximate marginal score function is centred at zero across the simulations. A non-centred score is a sign that the Laplace approximation is unreliable for this fit – typically because the random-effects posterior is far from Gaussian (saturated binomial / Beta, sparse counts, weakly identified random effects) or the data don't constrain the random effects well.

Scope: local Laplace-consistency diagnostics for a fitted model. This is a diagnostic signal rather than a substitute for bootstrap, profile, or external Bayesian calibration. Broader calibration evidence is not yet available.

This is a complementary signal to sanity_multi() (which checks Hessian definiteness, gradient magnitude, convergence flags) and to a target-explicit known-DGP study of rank-selection behaviour. Slower than the first-line checks – n_sim likelihood evaluations – but the only diagnostic that targets the Laplace approximation itself rather than the optimisation outcome or the parameter identification.

A large marginal-score p-value is an absence of evidence against centring in this finite simulation; it is not a certificate that the Laplace approximation is accurate. If the result is small (p_value <= 0.05) or large flagged-parameter list), consider:

• Refitting with a richer random-effects structure (so the conditional posterior is more Gaussian).

• Switching to a parametric bootstrap CI (which does not depend on the Laplace approximation being a good fit).

• Validating against tmbstan::tmbstan(fit$tmb_obj) (the audit's recommended Bayesian-comparison path).

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")
res <- gllvmTMB_check_consistency(fit, n_sim = 50L, seed = 1)
res
res$marginal_bias
res$flagged_parameters

Value

An object of class gllvmTMB_check_consistency with components:

⁠$marginal_p_value⁠

Joint chi-squared p-value across all marginal-score components. NA if TMB's information matrix could not be inverted (often the case on tiny fixtures).

⁠$marginal_bias⁠

Named numeric vector: per-parameter bias (mean of the simulated score, normalised by its SE). Large absolute values flag specific parameters where the Laplace approximation is biased.

⁠$joint_p_value⁠

Joint-score p-value. NA unless estimate = TRUE. Slower because each replicate refits.

⁠$flagged_parameters⁠

Character names of parameters whose marginal bias exceeds ⁠|bias| > 0.5⁠ (a heuristic; finite-n_sim noise can produce small spurious biases, but |0.5| is generally beyond noise for n_sim = 100).

⁠$n_sim⁠

The n_sim actually run.

⁠$threshold⁠

The marginal-bias threshold used.

⁠$diagnostics⁠

Character vector of one or more of: "no_centring_warning" (the finite simulation did not detect non-centring; this is not an accuracy certificate), "marginal_score_non_centred", "joint_score_non_centred", "information_matrix_singular" (TMB couldn't invert it; interpret with caution), "marginal_p_value_unavailable" (TMB returned NA without a captured warning; usually happens on tiny / weakly-identified fixtures).

⁠$raw⁠

The full TMB::checkConsistency() return value in case the user wants the per-replicate gradient matrix.

⁠$call⁠

match.call() of the invocation.

See Also

sanity_multi() (structural / convergence checks), check_auto_residual() (mixed-family safeguard), gllvmTMB_diagnose() (holistic fit summary), TMB::checkConsistency() (the underlying TMB call).

Diagnose a fitted model and suggest next actions

Description

This is the human-readable diagnostic to call right after fit <- gllvmTMB(...). It combines the quick numerical screen (sanity_multi()), the rotation identifiability advisory, and key biological summaries (correlation diagonals, ICCs, communalities) into a single report with explicit next-step hints for any WARN signal. Use check_gllvmTMB() when you need the same fit-health checks as a stable table for scripts or reports.

Usage

gllvmTMB_diagnose(
  object,
  gradient_thresh = 0.01,
  se_thresh = 100,
  big_corr_thresh = 0.5,
  verbose = TRUE
)

Arguments

Details

Scope: first-line convergence, Hessian, standard-error, restart, and rotation diagnostics. It reports risks and summaries but does not replace profile, bootstrap, or simulation calibration. Target-explicit validation will decide which interval warnings can be promoted to broader guarantees.

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")
gllvmTMB_diagnose(fit)

Value

Invisibly a list with components: sanity (the sanity_multi() flags), rotation (rotation-advisory list), Sigma_B, Sigma_W, ICC_site, communality_B, communality_W, and hints (character vector of suggested actions).

See Also

check_gllvmTMB(), sanity_multi(), suggest_lambda_constraint(), extract_Sigma(), extract_communality(), compare_dep_vs_two_psi() / compare_indep_vs_two_psi() for identifiability cross-checks on the paired phylogenetic fit.

Methods for Julia bridge fits

Description

Small S3 surface for fits returned by gllvmTMB(..., engine = "julia") or gllvm_julia_fit(). These methods expose the flat bridge payload that has already passed the R admission gates. They are point-estimate summaries: prediction and ordinary response/Pearson residuals are in-sample retained-payload reconstructions only. Simulation is conditional on retained fitted values for admitted one-family rows and complete balanced mixed-family rows. Unit-tier covariance honors native link_residual scale semantics, and raw ordination accessors are routed; richer extractor parity remains a separate bridge row. Confidence intervals are routed for admitted no-X Gaussian, Poisson, Bernoulli binomial, NB2, NB1, Beta, and Gamma rows; response-mask CIs are routed for the same non-Gaussian rows when X = NULL; complete-response fixed-effect-X CIs are routed for Gaussian, Poisson, Bernoulli binomial, NB2, Beta, and Gamma rows. They may be requested at fit time through gllvmTMB(ci_method = ...), or retrieved and recomputed through confint(). Predictor-informed latent(..., lv = ~ x) rows admit ci_method = "wald" at fit time; retained Wald payloads are surfaced as std.error, lower, and upper through extract_lv_effects() rather than confint().

Usage

## S3 method for class 'gllvmTMB_julia'
logLik(object, ...)

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

## S3 method for class 'gllvmTMB_julia'
coef(object, ...)

## S3 method for class 'gllvmTMB_julia'
predict(
  object,
  newdata = NULL,
  type = c("link", "response", "prob", "class"),
  ...
)

## S3 method for class 'gllvmTMB_julia'
fitted(object, type = c("response", "link", "prob", "class"), ...)

## S3 method for class 'gllvmTMB_julia'
residuals(object, type = c("response", "pearson"), ...)

## S3 method for class 'gllvmTMB_julia'
simulate(
  object,
  nsim = 1,
  seed = NULL,
  newdata = NULL,
  condition_on_RE = TRUE,
  ...
)

## S3 method for class 'gllvmTMB_julia'
confint(
  object,
  parm,
  level = 0.95,
  method = c("stored", "wald", "profile", "bootstrap"),
  ci_nboot = 200L,
  ci_seed = 0L,
  ...
)

## S3 method for class 'gllvmTMB_julia'
summary(object, ...)

## S3 method for class 'summary.gllvmTMB_julia'
print(x, digits = 3, ...)

Arguments

Value

logLik() returns a "logLik" object. coef() returns a named list of available point-estimate components. confint() returns a conventional confidence-interval matrix for stored or recomputed Julia CI payloads. predict() returns an in-sample data frame. Non-ordinal rows use trait, unit, and est columns; ordinal probability rows use trait, unit, category, and prob; ordinal class rows use trait, unit, and est. fitted() returns the in-sample fitted matrix with traits in rows and units in columns for non-ordinal rows, and an array of category probabilities or a modal-class matrix for ordinal rows. residuals() returns an in-sample residual matrix with the same shape as fitted() for non-ordinal rows and keeps masked response cells as NA. simulate() returns an ⁠n_obs x nsim⁠ matrix in the same trait-major cell order as predict() and keeps masked response cells as NA. summary() returns a list with header, coefficients, covariance, and status fields.

Methods on a fitted gllvmTMB model

Description

Standard model-object accessors for a multivariate fit returned by gllvmTMB(), whether the call started from wide traits(...) data or already-stacked long data. Internally the fit has class gllvmTMB_multi, which is what these S3 methods dispatch on, but you call print(fit) and summary(fit) as usual. Likelihood and inference accessors depend on the estimator, as described below.

Usage

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

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

## S3 method for class 'summary.gllvmTMB_multi'
print(x, digits = 3, ...)

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

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

Arguments

Details

• For ML fits, print() shows the active covstructs, the number of fixed effects, and the converged log-likelihood.

• For ML fits, summary() adds a fixed-effects table with SEs, the global and local trait correlation matrices, per-trait ICCs, and global / local communalities.

• For an unpenalised native-Laplace ML fit, logLik() is the converged maximum with df = length(opt$par) and nobs equal to the number of likelihood-contributing response cells. AGHQ has a distinct integration objective; a loading ridge is penalised MAP; and non-unit likelihood weights define an estimating objective. Ordinary likelihood comparison is unavailable on each of those restricted surfaces.

• For estimator = "mspl", print() and summary() identify the experimental softly penalised Laplace point estimator and show the unpenalised Laplace value at that point only as provenance. logLik(), AIC, BIC, likelihood-ratio tests, standard errors, intervals, and profiles fail closed because the point is not the ordinary likelihood maximum and repeated-sampling inference is not calibrated.

nobs() returns the number of likelihood-contributing observations – the observed-response cells. This equals fit$missing_data$counts$likelihood_rows and the nobs attribute of logLik(). Under the default miss_control(response = "drop") every fitted row is observed, so it equals length(fit$tmb_data$y); under response = "include" the masked rows are excluded. Original-row counts live in fit$missing_data, never in nobs().

Value

print() returns its input invisibly, after displaying the fit or summary. summary() returns a list of class summary.gllvmTMB_multi containing the model header, fixed-effect estimates, covariance summaries, and convergence diagnostics. logLik() returns a scalar of class logLik with df and nobs attributes for supported likelihood fits. nobs() returns the number of likelihood-contributing response cells as a numeric scalar. Estimator restrictions described below apply to these accessors.

Fixed-effect estimates and covariance for a multi-trait fit

Description

coef() returns the fixed-effect point estimates; vcov() returns their covariance matrix, taken from the fit's TMB sdreport().

Usage

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

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

Arguments

Details

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

coef(fit)
vcov(fit)

Value

coef(): a named numeric vector, one entry per fixed-effect term (object$X_fix_names).

vcov(): a square numeric matrix with those same names on both margins. Rows and columns for coefficients held fixed via Xcoef_fixed are NA – a fixed parameter was not estimated, so it has no sampling covariance.

Point estimates need no standard errors

coef() works on any fit, including one made with gllvmTMBcontrol(se = FALSE) – point estimates are this package's supported claim and do not depend on sdreport(). vcov() does depend on it, and raises the same typed errors confint() does when it is missing (gllvmTMB_confint_no_sdreport) or non-finite (gllvmTMB_confint_nonfinite_se), naming standard_errors() as the remedy in the first case. The conditions are shared deliberately: a caller that handles one should handle the other.

What this covariance is

The fixed-effect block of the single TMB sdreport() the fit already carries – a Wald covariance, with exactly the caveats any Wald quantity from this package carries. It is not a resampled or profiled quantity, and no interval built from it has certified coverage. vcov() refuses a fit with non-unit likelihood weights: that fit targets a weighted estimating objective and has no certified sandwich covariance.

See Also

standard_errors() to compute the sdreport() after fitting; confint() for intervals; summary() for a coefficient table.

Methods for variational (integration = "va") fits

Description

Methods for the object returned by gllvmTMB(..., control = gllvmTMBcontrol(integration = "va")).

Usage

## S3 method for class 'gllvmTMB_va'
print(x, digits = 3, ...)

## S3 method for class 'gllvmTMB_va'
summary(object, ...)

## S3 method for class 'summary.gllvmTMB_va'
print(x, digits = 3, ...)

## S3 method for class 'gllvmTMB_va'
nobs(object, ...)

## S3 method for class 'gllvmTMB_va'
logLik(object, ...)

## S3 method for class 'gllvmTMB_va'
confint(object, parm, level = 0.95, ...)

## S3 method for class 'gllvmTMB_va'
vcov(object, ...)

## S3 method for class 'gllvmTMB_va'
coef(object, ...)

## S3 method for class 'gllvmTMB_va'
residuals(object, ...)

## S3 method for class 'gllvmTMB_va'
fitted(object, ...)

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

## S3 method for class 'gllvmTMB_va'
df.residual(object, ...)

## S3 method for class 'gllvmTMB_va'
weights(object, ...)

## S3 method for class 'gllvmTMB_va'
predict(object, ...)

Arguments

Details

This is an opt-in research route. Its objective is an ELBO – a lower bound on the log-likelihood – not a log-likelihood, and its inverse Hessian is not automatically calibrated frequentist uncertainty. Accordingly logLik(), AIC(), and BIC() remain undefined. vcov() and confint() expose only fixed-effect VA-Wald uncertainty computed from the profiled Schur information matrix. The returned values are explicitly labelled uncalibrated; raw-loading Wald intervals are not available.

The ordination surface is the exception: extract_ordination(), getLV(), getLoadings(), and extract_loadings() all work for a variational fit and return latent scores and loadings as POINT ESTIMATES. getLV(se = TRUE) additionally returns a variational posterior SD – not a standard error – and only when the fit's eval_method resolves to "gh". See the integration argument of gllvmTMB() for the measured accuracy of these point estimates.

Value

print() and print.summary() return their argument invisibly. summary() returns an object of class "summary.gllvmTMB_va". nobs() returns an integer count of unit-by-response cells. vcov() returns the profiled-Schur beta covariance matrix and confint() its nominal Wald intervals; both carry calibrated = FALSE and uncertainty_basis = "VA-Wald profiled Schur information" attributes.

Fit a GLLVM from a wide unit × trait matrix

Description

Convenience wrapper that lets users supply a wide ⁠unit × trait⁠ response matrix Y (rows = units, columns = traits) and an optional unit-level predictor data frame X. Pivots to long format internally and dispatches to gllvmTMB().

Usage

gllvmTMB_wide(
  Y,
  X = NULL,
  d = 2,
  family = gaussian(),
  phylo_vcv = NULL,
  formula_extra = NULL,
  weights = NULL,
  REML = FALSE,
  ...
)

Arguments

Details

The "unit × trait" framing is generic: ⁠site × species⁠ (joint species distribution modelling), ⁠individual × trait⁠ (morphometrics / behavioural syndromes), ⁠species × trait⁠ (phylogenetic comparative), ⁠paper × outcome⁠ (meta-analysis), or any similar layout. The function does not assume the ecological special case; it is the unified matrix-in entry point for the stacked-trait GLLVM engine.

Value

A fit returned by gllvmTMB(). The column dimension of Y is exposed as the "trait" axis of the model (so extract_ordination() returns the trait loadings). For the site × species special case the columns are species, and the loadings are species loadings; the same machinery returns trait loadings for individual × trait morphometrics, study loadings for paper × outcome meta-analysis, etc.

Deprecation

gllvmTMB_wide() is soft-deprecated as of gllvmTMB 0.2.0. The recommended replacement is gllvmTMB() with the traits() LHS marker, which is more general (formula-native predictors, full covariance grammar, mixed-family fits) and matches the formula-first idiom used by lme4, glmmTMB, brms, drmTMB.

Migration:

## Old:
fit <- gllvmTMB_wide(Y, X = X_df, d = 2,
                     formula_extra = ~ env_temp + env_precip)

## New:
df_wide <- cbind(data.frame(unit = rownames(Y)),
                 as.data.frame(Y), X_df)
fit <- gllvmTMB(
  traits(<colnames of Y>) ~ 1 + env_temp + env_precip +
    latent(1 | unit, d = 2),
  data = df_wide,
  unit = "unit"
)

Per-cell weight matrices (the one path gllvmTMB_wide() uniquely supports) remain available via the long-format API by passing a long-format weights column aligned with ⁠(unit, trait)⁠ rows.

See Also

gllvmTMB() for the recommended formula-API entry point (long format or wide data frames via traits() LHS sugar); extract_Sigma() for post-fit covariance summaries; extract_ordination() for scores and loadings.

Control parameters for gllvmTMB()

Description

Control parameters for gllvmTMB()

Usage

gllvmTMBcontrol(
  d_B = NULL,
  d_W = NULL,
  spde_mode = c("per_trait", "shared"),
  n_init = 1L,
  optimizer = c("nlminb", "optim"),
  optArgs = list(),
  init_jitter = 0.3,
  init_strategy = c("default", "single_trait_warmup"),
  start_method = list(method = NULL, jitter.sd = 0),
  start_from = NULL,
  se = TRUE,
  verbose = FALSE,
  integration = c("laplace", "va"),
  va_H = 7L,
  va_eval_method = c("auto", "jj", "gh"),
  aghq = FALSE,
  aghq_iter_cap = 1L,
  aghq_n_adapt = 400L,
  aghq_ridge = 2,
  aghq_continuation = TRUE,
  aghq_shift_tol = 1e-04,
  aghq_grad_tol = 1e-04,
  aghq_grad_tol_rel = 1e-06,
  aghq_multistart = TRUE,
  aghq_f_tol = 1e-09,
  aghq_escalate_patience = 3L,
  aghq_rho_min = 1/64,
  warn_runaway = TRUE,
  allow_nongaussian_reml = FALSE,
  loading_ridge = NULL,
  ...
)

Arguments

Details

Recommended workflow for two-level latent models. Reduced-rank likelihoods are notoriously multimodal in two-level fits like value ~ 0 + trait + latent(0+trait|ID, d_B) + latent(0+trait|obs_ID, d_W). The recommended workflow is:

1. Run several restarts (n_init = 5 to 10) and keep the one with the lowest -logLik.

2. If that fails to converge cleanly, switch the optimiser to optim with BFGS: gllvmTMBcontrol(optimizer = "optim", optArgs = list(method = "BFGS")).

3. start_method = list(method = "res", jitter.sd = 0.2) is soft-deprecated as of 0.6.0 and warns once per session. It still fits: it decomposes the fixed-effect residual matrix into starting values for \Lambda and the latent scores. Do not reach for it.

   It was retired on measurement. Across 89 simulated fits (Gaussian, Poisson and negative-binomial, d = 1 to 3, three and five traits), it was never materially better than the default start (its three best margins were 0.07, 0.29 and 0.66 log-likelihood units) and was materially worse eight times, once by 14.6 units. At d >= 2 it made no difference at all. Every failure reported convergence == 0 and a positive-definite Hessian on both sides, so the worse fit was silent, and restarts did not detect it: in one such fit all five n_init = 5 restarts returned the same worse optimum.

   If you use it anyway, compare fit$opt$objective against a default-start fit and keep the lower value.

4. For Gaussian two-level models, prefer start_method = list(method = "indep") or manually fit a simpler model and pass it through start_from = simpler_fit. This is a GLMM warm start rather than a fixed-effect-only GLM residual start.

Value

A gllvmTMBcontrol list.

Examples

# gllvmTMBcontrol() is a pure-R constructor: it builds and validates a
# control list without fitting anything.
gllvmTMBcontrol()                                  # historical defaults

# Multi-start to guard against multimodal reduced-rank likelihoods, with
# jitter on the starting parameter vector across restarts.
gllvmTMBcontrol(n_init = 3, init_jitter = 0.2)

# Switch the optimiser to optim + BFGS for finicky two-level rr fits.
gllvmTMBcontrol(optimizer = "optim", optArgs = list(method = "BFGS"))

# Experimental, opt-in scale-aware ridge for the evaluated Bernoulli AGHQ
# scope. The route uses a fixed 9-node multi-start rule; this is not an
# interval-calibration certificate.
gllvmTMBcontrol(aghq_ridge = "auto")
# The same fixed loading MAP penalty, named independently of integration.
gllvmTMBcontrol(loading_ridge = 2)

# Experimental scalar variational route: public auto uses GH with H = 7.
gllvmTMBcontrol(integration = "va")
# H = 61 remains a diagnostic; JJ is explicit binomial-logit only.
gllvmTMBcontrol(integration = "va", va_H = 61, va_eval_method = "gh")

Fit a GLLVM with the Julia engine (GLLVM.jl bridge_fit).

Description

Fit a GLLVM with the Julia engine (GLLVM.jl bridge_fit).

Usage

gllvm_julia_fit(
  y,
  family = "gaussian",
  num.lv = 2L,
  N = NULL,
  X = NULL,
  X_lv = NULL,
  coef_fixed = NULL,
  mask = NULL,
  units_are_rows = FALSE,
  ci_method = c("none", "wald", "profile", "bootstrap"),
  ci_level = 0.95,
  ci_nboot = 200L,
  ci_seed = 0L,
  ...
)

Arguments