Choose a response family

Choose a family from the values the response can take and the process that generated those values. The family sets the likelihood and link; covariance terms such as latent() or indep() describe dependence among responses. Changing the family can therefore require different response coding, diagnostics, and interpretation—not merely a different argument value.

Start with the simplest scientifically defensible family, inspect fit and predictive diagnostics, and change it when the observed mean–variance or zero pattern contradicts its assumptions.

Start from the response

Continuous responses

Observed values and sampling process	Family and link	Main assumption or warning
Any real value; roughly symmetric conditional errors	gaussian(); identity	Residual variation is Gaussian on the response scale.
Any real value; occasional genuine heavy-tailed observations	student(); identity	df = NULL estimates degrees of freedom; student(df = 3) fixes them. The fitted scale is not the response SD, and the variance is infinite when df <= 2.
Strictly positive; multiplicative errors are plausible	lognormal(); log	The linear predictor describes the mean of log(y), so exp(eta) is the conditional median, not the arithmetic mean.
Strictly positive; variance grows approximately with the square of the mean	Gamma(link = "log"); log	Uses a mean–shape parameterization with constant conditional coefficient of variation within a trait.
Non-negative continuous values with exact zeros from one compound process	tweedie(); log	Assumes Var(Y) = phi * mu^p, with 1 < p < 2. It is a one-likelihood alternative to an explicit occurrence/positive decomposition.
Non-negative continuous values with a scientifically meaningful occurrence process and positive part	delta_lognormal() or delta_gamma(); logit + log	The current standard hurdle model uses the same linear predictor for occurrence and positive mean. See the specialist boundary below before using covariance terms.

Gamma() and lognormal() can both fit right-skewed positive data. Choose from the sampling mechanism and residual pattern rather than the histogram alone: lognormal errors are additive after taking logs, whereas Gamma variation is described through a conditional coefficient of variation.

For exact zeros, ask whether one compound process is plausible or whether zero versus positive observations represent distinct biological events. That scientific distinction is more important than which likelihood gives the smallest AIC in one dataset.

Binary, trial, proportion, and ordinal responses

Observed values and sampling process	Family and link	Main assumption or warning
A 0/1 event or successes out of known trials	binomial(); logit, probit, or cloglog	Use cbind(successes, failures) for trials. A flat response plus weights = n_trials is also accepted.
Successes out of known trials with extra-binomial variation	betabinomial(); logit	Uses the same trial encoding as binomial data and estimates a precision parameter for each trait.
A continuous proportion strictly inside (0, 1)	Beta(); logit	Genuine zeros or ones are not Beta observations; do not silently nudge them inward without a scientific measurement model.
Ordered categories coded consecutively as 1, ..., K, with K >= 3	ordinal_probit(); probit	Category order is part of the model. For two categories, use binomial(link = "probit"). Recover estimated thresholds with extract_cutpoints().

ordinal_probit() uses an underlying liability

$$y^* = \eta + \varepsilon, \qquad \varepsilon \sim N(0, 1),$$

and records category $k$ when $\tau_{k-1} < y^* \leq \tau_k$. gllvmTMB fixes $\tau_1 = 0$ as its location convention and estimates the remaining $K-2$ free boundaries. extract_cutpoints(fit) returns those boundaries from cutpoint_2 onward. A cutpoint is a boundary between adjacent categories on the latent probit scale, not a regression effect or an observed category-unit change. A predictor shifts the same latent liability across all boundaries, which is the proportional-threshold assumption.

Prepare ordinal responses deliberately. Use meaningful ordered categories coded as consecutive integers 1, ..., K; a numeric code does not make a nominal category ordinal. The current engine infers $K$ from the largest observed code for each trait, so the highest intended category must occur in the fitted data and gaps must be resolved before fitting. Two-category input is accepted as the probit special case, but new binary analyses should use binomial(link = "probit") directly.

This lookup stops at family choice and cutpoint interpretation. Ordinal variance, correlation, or heritability requires a replicated or externally structured design that identifies the relevant variance tier; a converged fit with one random effect per observation is not evidence that such a component has been recovered.

For multi-trial binomial and beta-binomial data, weights means trial counts. For other families, weights are likelihood multipliers. Prefer the explicit cbind(successes, failures) form when possible because it makes the denominator visible in the formula.

Counts

Sampling process	Family and link	Conditional variance
Counts with no remaining dispersion after the mean structure and random effects	poisson(); log	Var(Y) = mu
Overdispersion that grows approximately linearly with the mean	nbinom1(); log	Var(Y) = mu * (1 + phi)
Overdispersion that grows approximately quadratically with the mean	nbinom2(); log	Var(Y) = mu + mu^2 / phi
Positive counts because zero cannot be observed by design	truncated_poisson() or truncated_nbinom2(); log	Conditions on Y >= 1; every response must be a positive integer.

Do not choose a zero-truncated family merely because zeros were deleted during cleaning. It is appropriate only when the observation process could not produce or record a zero. Likewise, an observation-level random effect is not an automatic repair for Poisson overdispersion: first check whether the missing structure is better represented by NB1/NB2, a predictor, a grouping term, or a different sampling process.

One family in long and wide data

A single family applies to every response in the fit. The long and traits(...) wide calls below fit the same Poisson model; only the data shape and formula shorthand differ.

library(gllvmTMB)

set.seed(20260711)
n_unit <- 40L
trait_names <- c("beetles", "moths", "spiders")
unit <- factor(seq_len(n_unit))
z <- rnorm(n_unit, sd = 0.45)

eta <- cbind(
  beetles = 0.2 + 0.8 * z,
  moths = 0.6 - 0.4 * z,
  spiders = 0.1 + 0.5 * z
)
Y <- matrix(
  rpois(length(eta), lambda = exp(eta)),
  nrow = n_unit,
  dimnames = list(NULL, trait_names)
)

counts_wide <- data.frame(unit = unit, Y, check.names = FALSE)
counts_long <- data.frame(
  unit = rep(unit, each = length(trait_names)),
  trait = factor(rep(trait_names, times = n_unit), levels = trait_names),
  count = as.vector(t(Y))
)

fit_long <- gllvmTMB(
  count ~ 0 + trait + latent(0 + trait | unit, d = 1),
  data = counts_long,
  trait = "trait",
  unit = "unit",
  family = poisson(),
  silent = TRUE
)
#> Warning: ! Ordinary `latent()` now includes a per-trait Psi by default (Sigma = Lambda
#>   Lambda^T + Psi).
#> ℹ This changed in gllvmTMB 0.2.0; earlier `latent()` was loadings-only (Lambda
#>   Lambda^T).
#> → Pass `latent(..., unique = FALSE)` for the old rotation-invariant
#>   loadings-only fit.

fit_wide <- gllvmTMB(
  traits(beetles, moths, spiders) ~ 1 + latent(1 | unit, d = 1),
  data = counts_wide,
  unit = "unit",
  family = poisson(),
  silent = TRUE
)

data.frame(
  input = c("long", "wide"),
  logLik = c(as.numeric(logLik(fit_long)), as.numeric(logLik(fit_wide)))
)
#>   input    logLik
#> 1  long -184.6423
#> 2  wide -184.6423

Equivalent long and wide designs should give the same likelihood. This check does not show that Poisson is the right family; that requires diagnostics and a scientific mean–variance argument.

Different families in one fit

A mixed-family model needs long data because the family selector belongs to each stacked response row. Use a factor with explicit levels and a named family list; names are matched to selector levels, so accidental list ordering cannot silently swap likelihoods.

mixed_long$family <- factor(
  mixed_long$family,
  levels = c("continuous", "presence", "count")
)

fam <- list(
  continuous = gaussian(),
  presence = binomial(),
  count = poisson()
)
attr(fam, "family_var") <- "family"

fit_mixed <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | unit, d = 1),
  data = mixed_long,
  trait = "trait",
  unit = "unit",
  family = fam
)

An unnamed list is still accepted, but its order must match the selector factor levels exactly. A named list is safer. Keep one family within each trait; in particular, an ordinal trait must own all of its rows because its thresholds are estimated per trait.

Check that rule explicitly before fitting:

families_per_trait <- tapply(
  as.character(mixed_long$family),
  mixed_long$trait,
  function(x) length(unique(x))
)
stopifnot(all(families_per_trait == 1L))

Mixed-family fitting does not make every response directly comparable on its raw observed scale. Define the scale and scientific target before interpreting cross-trait covariance.

Covariance scale and uncertainty

The fitted covariance tier and the observed responses are different objects. For extract_Sigma():

• link_residual = "none" returns the covariance components fitted in the selected tier;
• link_residual = "auto" adds a family-specific diagonal convention before converting to correlations.

The automatic addition is exact for fixed-scale links such as probit and an approximation for several count, proportion, and two-part families. Neither setting generally gives the correlation between raw observed responses. Choose between them by stating the estimand, not by treating "auto" as universally better.

For a fitted mixed-family model, keep the two point-estimate targets separate:

# Covariance fitted at the selected random-effect tier. Entries retain the
# units of their trait-specific linear predictors.
structural <- extract_Sigma(
  fit_mixed,
  level = "unit",
  part = "total",
  link_residual = "none"
)

# Pairwise family/link-scale correlations after adding each trait's
# link-residual convention to its diagonal.
adjusted <- extract_Sigma(
  fit_mixed,
  level = "unit",
  part = "total",
  link_residual = "auto"
)

structural$Sigma
adjusted$R
adjusted$note

Raw covariance magnitudes are not directly comparable when one predictor is in identity-link units, another in logit units, and another in log units. The correlations in adjusted$R are dimensionless, but they remain conditional on the family, link, fitted mean trajectory, and residual convention for each trait. For example, the Poisson addition uses a trait-level approximation based on its fitted mean. These are pairwise model/link-scale summaries—not a single universal latent scale and not observed-response correlations.

For mixed-family fits, point estimates can describe dependence on the selected latent or link scale, but broad interval coverage has not been established for this target. Keep the public workflow point-estimate-only; do not use mixed-family correlation intervals as calibrated uncertainty.

Specialist boundary for hurdle families

delta_lognormal() and delta_gamma() currently use one shared predictor:

$$\Pr(Y > 0) = \operatorname{logit}^{-1}(\eta), \qquad E(Y \mid Y > 0) = \exp(\eta).$$

Use these constructors only when that shared-predictor assumption is scientifically defensible. The current reader-facing route is the standard fixed-effect hurdle fit. The package does not yet supply a response-scale covariance or correlation interpretation for hurdle fits, so do not describe fitted-tier correlations as correlations of total biomass, abundance, or another observed two-part response.

Diagnose the fitted choice

Family choice is provisional until the fitted model behaves sensibly.

check_gllvmTMB(fit_long)

# Exact randomized-quantile residuals are currently available for
# Gaussian, Poisson, NB1, and NB2 rows.
residuals(fit_long, type = "randomized_quantile", seed = 1)

# Rootograms currently support Poisson and NB2 count rows.
predictive_check(fit_long, type = "rootogram", ndraws = 100, seed = 1)

Systematic residual or rootogram structure can indicate a wrong family, link, mean model, offset, or random-effect structure. It does not identify which one automatically. Compare alternatives that represent credible sampling processes, then recheck convergence and prediction rather than choosing solely by AIC.

Common input failures have direct remedies:

• Gamma or lognormal data contain zero: revisit the measurement process; use a family that allows zero rather than adding an arbitrary constant.
• Counts are negative or non-integer: correct the response definition; a count likelihood is not appropriate for transformed or continuous values.
• Beta data include 0 or 1: model the boundary process explicitly or choose another response representation.
• Ordinal categories are not consecutive integers starting at 1: recode an ordered factor with as.integer() and verify the level order.
• A link is rejected: use the accepted link shown above; the constructor may allow a base-R link that the multivariate engine does not.
• Mixed-family names do not match selector levels: name every family-list element exactly as the factor level.
• A constructor is absent from this page: it is not a supported multivariate response family; choose an available likelihood rather than relying on an exported experimental constructor.

For the full diagnostic workflow, continue to Fit diagnostics. For covariance syntax, see the Formula keyword grid.