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.

First workflow: Poisson or NB2 for insect counts?

Suppose pitfall traps record aphids, thrips, and weevils along an environmental gradient. Poisson is a sensible starting point because the responses are counts. It also makes a strong assumption: after accounting for the gradient, the conditional variance equals the conditional mean.

The example below deliberately simulates extra variation that grows with the square of the mean. We first prepare the same counts that an ecological survey would provide, without using the known simulation settings during fitting.

library(gllvmTMB)

set.seed(20260817)
n_unit <- 200L
trait_names <- c("aphids", "thrips", "weevils")
unit <- factor(seq_len(n_unit))
environment <- rnorm(n_unit, sd = 0.45)

eta <- cbind(
  aphids = 1.0 + 0.8 * environment,
  thrips = 1.4 - 0.4 * environment,
  weevils = 0.8 + 0.5 * environment
)
phi_true <- c(aphids = 1.5, thrips = 2, weevils = 1.2)
Y <- matrix(
  rnbinom(
    length(eta),
    mu = exp(eta),
    size = rep(phi_true, each = n_unit)
  ),
  nrow = n_unit,
  dimnames = list(NULL, trait_names)
)

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

Fit Poisson and NB2 with the same trait-specific environmental effects. This keeps the family comparison focused on the conditional count distribution, before adding latent or other covariance terms.

fit_poisson <- gllvmTMB(
  count ~ 0 + trait + trait:environment,
  data = counts_long,
  trait = "trait",
  unit = "unit",
  family = poisson(),
  silent = TRUE
)

fit_nb_long <- gllvmTMB(
  count ~ 0 + trait + trait:environment,
  data = counts_long,
  trait = "trait",
  unit = "unit",
  family = nbinom2(),
  silent = TRUE
)

Randomized-quantile residuals should have mean near 0 and standard deviation near 1 when the conditional distribution is adequate. Compare that diagnostic with the fitted log-likelihood:

resid_poisson <- residuals(
  fit_poisson,
  type = "randomized_quantile",
  seed = 1
)
resid_nb2 <- residuals(
  fit_nb_long,
  type = "randomized_quantile",
  seed = 1
)

data.frame(
  family = c("Poisson", "NB2"),
  logLik = round(c(
    as.numeric(logLik(fit_poisson)),
    as.numeric(logLik(fit_nb_long))
  ), 1),
  residual_mean = round(c(
    mean(resid_poisson$residual),
    mean(resid_nb2$residual)
  ), 2),
  residual_sd = round(c(
    sd(resid_poisson$residual),
    sd(resid_nb2$residual)
  ), 2)
)
#>    family  logLik residual_mean residual_sd
#> 1 Poisson -1459.8         -0.11        1.48
#> 2     NB2 -1291.0          0.00        0.99

The Poisson residual standard deviation is about 1.48, so the counts are much more variable than its conditional distribution allows. NB2 brings that value to about 0.99, keeps the residual mean near 0, and has the better likelihood. For these data, choose NB2. In real data, also plot residuals against fitted means and revisit the predictors or grouping structure if a pattern remains. Do not choose NB2 only because it has the larger likelihood: the ecological sampling process and the residual pattern should support the change.

The same NB2 model in long and wide data

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

counts_wide <- data.frame(
  unit = unit,
  environment = environment,
  Y,
  check.names = FALSE
)

fit_nb_wide <- gllvmTMB(
  traits(aphids, thrips, weevils) ~ 1 + environment,
  data = counts_wide,
  unit = "unit",
  family = nbinom2(),
  silent = TRUE
)

data.frame(
  input = c("long", "wide"),
  logLik = c(
    as.numeric(logLik(fit_nb_long)),
    as.numeric(logLik(fit_nb_wide))
  )
)
#>   input    logLik
#> 1  long -1291.038
#> 2  wide -1291.038

Equivalent long and wide designs give the same likelihood. Choose the shape that best matches the rest of the analysis; it does not change the family or its interpretation.

Interpret the NB2 dispersion

When counts vary more than a Poisson mean allows even after the mean structure and covariance terms, nbinom2() adds one dispersion parameter per trait. Its conditional variance is

$\operatorname{Var}(Y) = \mu + \frac{\mu^2}{\phi},$

so $\phi$ controls the quadratic extra-Poisson variance: smaller $\phi$ means stronger overdispersion, and as $\phi \to \infty$ the family collapses back to Poisson. This $\phi$ is the negative binomial size parameter on the natural scale—the value passed to rnbinom(mu = , size = ) in the first workflow—not a variance, a standard deviation, or a log-scale quantity.

The fitted dispersions live in fit$report$phi_nbinom2, one value per trait in trait factor-level order, already on the natural $\phi$ scale:

data.frame(
  trait = trait_names,
  phi_true = as.numeric(phi_true),
  phi_hat = round(as.numeric(fit_nb_long$report$phi_nbinom2), 2)
)
#>     trait phi_true phi_hat
#> 1  aphids      1.5    2.01
#> 2  thrips      2.0    2.00
#> 3 weevils      1.2    1.45

Read these against the variance formula, not as effect sizes. For example, $\hat\phi \approx 2$ gives extra variance $\mu^2 / 2$ , which already exceeds the Poisson variance when $\mu = 5$ . When the data contain little quadratic overdispersion, $\phi$ can be weakly identified and run away to a huge value (1e6 or more). A runaway estimate means “this trait looks Poisson after the modelled predictors and random effects”, not a precise dispersion estimate; refit with poisson() and compare rather than reporting the number.

This example demonstrates point estimation and fitted diagnostics without a covariance tier. Interval coverage for negative binomial fits has not been certified, and NB2 latent models at small sample sizes have documented dispersion failures. See Current limitations and boundaries before adding a latent term or reporting uncertainty from a negative binomial fit.

Gamma’s `phi_gamma` is a shape, not a dispersion

Reach for Gamma(link = "log") when a positive response has a roughly constant coefficient of variation as its mean changes. The trap: fit$report$phi_gamma is the shape of the mean-shape parameterization, not a variance or an SD, and a larger value means less relative variability. Convert with scale = mu / phi_gamma and CV(y) = 1 / sqrt(phi_gamma).

# Shared plumbing for the compact fits below: two traits per unit in long
# format, the `value ~ 0 + trait + latent(0 + trait | unit)` shape used
# throughout this article.
two_trait_long <- function(unit, trait_names, y_matrix, response = "value") {
  out <- data.frame(
    unit = rep(unit, each = 2),
    trait = factor(rep(trait_names, times = length(unit)), levels = trait_names)
  )
  out[[response]] <- as.vector(t(y_matrix))
  out
}

set.seed(20260818)
n_gam <- 150L
gam_traits <- c("leaf_area", "root_mass")
z_gam <- rnorm(n_gam, sd = 0.4)
eta_gam <- cbind(leaf_area = 1.5 + 0.6 * z_gam, root_mass = 1.0 - 0.3 * z_gam)
shape_true_gam <- c(leaf_area = 6, root_mass = 3)
Y_gam <- matrix(
  rgamma(length(eta_gam), shape = rep(shape_true_gam, each = n_gam),
         scale = exp(eta_gam) / rep(shape_true_gam, each = n_gam)),
  nrow = n_gam, dimnames = list(NULL, gam_traits)
)
gam_long <- two_trait_long(factor(seq_len(n_gam)), gam_traits, Y_gam)

fit_gam <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | unit, d = 1),
  data = gam_long, trait = "trait", unit = "unit",
  family = Gamma(link = "log"), 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.
phi_gam_hat <- fit_gam$report$phi_gamma
data.frame(trait = gam_traits, shape_true = shape_true_gam,
           phi_hat = round(phi_gam_hat, 2),
           cv_hat = round(1 / sqrt(phi_gam_hat), 2))
#>               trait shape_true phi_hat cv_hat
#> leaf_area leaf_area          6    7.79   0.36
#> root_mass root_mass          3    3.23   0.56

leaf_area’s larger phi_hat (about 7.8, shape simulated at 6) corresponds to the smaller CV (about 0.36); root_mass’s smaller phi_hat (about 3.2) has the larger CV (about 0.56). Report the CV, not phi_hat itself, whenever “more dispersed” is the claim you want to make.

Student-t’s scale is not the response SD

Reach for student() when a real-valued response has occasional genuine heavy-tailed observations that a Gaussian residual would under-predict. The trap: fit$report$sigma_student is the scale of the underlying $t$ -distribution, not the response SD, and the two diverge as df shrinks: $\mathrm{SD}(y) = \sigma \sqrt{df / (df - 2)}$ , undefined for $df \leq 2$ . Leave df = NULL to estimate degrees of freedom, or fix them with student(df = 3).

set.seed(20260818)
n_stu <- 150L
stu_traits <- c("bill_len", "wing_len")
z_stu <- rnorm(n_stu, sd = 0.4)
eta_stu <- cbind(bill_len = 10 + 1.5 * z_stu, wing_len = 15 - 1.0 * z_stu)
sigma_true_stu <- c(bill_len = 0.8, wing_len = 0.5)
df_fix <- 5
Y_stu <- matrix(
  eta_stu + rt(length(eta_stu), df = df_fix) * rep(sigma_true_stu, each = n_stu),
  nrow = n_stu, dimnames = list(NULL, stu_traits)
)
stu_long <- two_trait_long(factor(seq_len(n_stu)), stu_traits, Y_stu)

fit_stu <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | unit, d = 1),
  data = stu_long, trait = "trait", unit = "unit",
  family = student(df = df_fix), silent = TRUE
)
sigma_stu_hat <- fit_stu$report$sigma_student
data.frame(trait = stu_traits, sigma_true = sigma_true_stu,
           sigma_hat = round(sigma_stu_hat, 2),
           sd_hat = round(sigma_stu_hat * sqrt(df_fix / (df_fix - 2)), 2))
#>             trait sigma_true sigma_hat sd_hat
#> bill_len bill_len        0.8      0.80   1.03
#> wing_len wing_len        0.5      0.52   0.67

At df = 5, bill_len’s recovered scale (about 0.80) converts to an SD of about 1.03—29% larger than the scale itself, and the gap widens as df approaches 2 from above.

Beta’s `phi_beta` is a precision, and boundary values are not Beta data

Reach for Beta() when a response is a continuous proportion strictly inside (0, 1). fit$report$phi_beta is a precision, not a variance: larger values mean less spread, via shape1 = mu * phi_beta, shape2 = (1 - mu) * phi_beta with mu = plogis(eta). The other trap is at the boundary: genuine 0s or 1s are not Beta observations, and nudging them inward to make the likelihood accept them silently changes the measurement model rather than fixing a data-entry problem.

set.seed(20260818)
n_bet <- 150L
bet_traits <- c("cover_a", "cover_b")
z_bet <- rnorm(n_bet, sd = 0.4)
eta_bet <- cbind(cover_a = 0.2 + 0.8 * z_bet, cover_b = -0.3 - 0.5 * z_bet)
mu_bet <- plogis(eta_bet)
phi_true_bet <- c(cover_a = 12, cover_b = 20)
Y_bet <- matrix(
  rbeta(length(eta_bet), shape1 = mu_bet * rep(phi_true_bet, each = n_bet),
        shape2 = (1 - mu_bet) * rep(phi_true_bet, each = n_bet)),
  nrow = n_bet, dimnames = list(NULL, bet_traits)
)
bet_long <- two_trait_long(factor(seq_len(n_bet)), bet_traits, Y_bet)

fit_bet <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | unit, d = 1),
  data = bet_long, trait = "trait", unit = "unit",
  family = Beta(), silent = TRUE
)
data.frame(trait = bet_traits, phi_true = phi_true_bet,
           phi_hat = round(fit_bet$report$phi_beta, 2))
#>           trait phi_true phi_hat
#> cover_a cover_a       12   10.94
#> cover_b cover_b       20   22.24

Both precisions recover within simulation noise of their targets (about 10.9 against 12, and 22.2 against 20). If your data legitimately contain 0 or 1, model the boundary process explicitly rather than reporting a Beta fit against nudged values.

Multi-trial binomial: `cbind()` vs `weights`, and a `simulate()` gotcha

Reach for binomial() with a trial count whenever the response is successes out of known trials rather than a single 0/1 event. Two encodings fit identically: cbind(successes, failures) ~ ..., or a flat response of raw success counts with weights = n_trials—weights is the trial count, not a proportion, and passing successes / n_trials silently fits the wrong likelihood.

set.seed(20260818)
n_trl <- 150L
trl_traits <- c("infected", "damaged")
n_trials <- 10L
z_trl <- rnorm(n_trl, sd = 0.4)
eta_trl <- cbind(infected = 0.1 + 0.7 * z_trl, damaged = -0.4 - 0.5 * z_trl)
succ_trl <- matrix(
  rbinom(length(eta_trl), size = n_trials, prob = plogis(eta_trl)),
  nrow = n_trl, dimnames = list(NULL, trl_traits)
)
trl_long <- two_trait_long(factor(seq_len(n_trl)), trl_traits, succ_trl,
                            "successes")
trl_long$failures <- n_trials - trl_long$successes

fit_trl_cbind <- gllvmTMB(
  cbind(successes, failures) ~ 0 + trait + latent(0 + trait | unit, d = 1),
  data = trl_long, trait = "trait", unit = "unit",
  family = binomial(), silent = TRUE
)
fit_trl_weights <- gllvmTMB(
  successes ~ 0 + trait + latent(0 + trait | unit, d = 1),
  data = trl_long, trait = "trait", unit = "unit",
  weights = rep(n_trials, nrow(trl_long)), family = binomial(), silent = TRUE
)
c(cbind_logLik = as.numeric(logLik(fit_trl_cbind)),
  weights_logLik = as.numeric(logLik(fit_trl_weights)))
#>   cbind_logLik weights_logLik 
#>      -583.0369      -583.0369

The two log-likelihoods match exactly (about -583.04). Separately: simulate() treated every multi-trial binomial response as Bernoulli regardless of trial count until it was fixed on 2026-08-17; any saved simulation output for a multi-trial binomial fit produced before that fix is suspect and should be regenerated.

Ordinal cutpoints depend on which categories appear

Reach for ordinal_probit() when categories are ordered and coded as consecutive integers 1, ..., K. gllvmTMB fixes $\tau_1 = 0$ as its location convention, and extract_cutpoints() returns only $\tau_2, \ldots, \tau_{K-1}$ . The hazard: K is inferred per trait from the largest observed category, not declared, so a category that happens not to occur in the fitted data silently shrinks that trait’s model.

set.seed(20260818)
n_ord <- 150L
ord_traits <- c("condition", "vigor")
z_ord <- rnorm(n_ord, sd = 0.4)
taus_true <- c(0, 1.0, 2.2) # tau_1 = 0, tau_2 = 1.0, tau_3 = 2.2 -> K = 4
eta_ord <- cbind(condition = 0.5 + 0.9 * z_ord, vigor = 0.2 - 0.6 * z_ord)
draw_ordinal <- function(eta, taus) {
  ystar <- eta + rnorm(length(eta))
  as.integer(cut(ystar, breaks = c(-Inf, taus, Inf), labels = FALSE))
}
Y_ord <- matrix(draw_ordinal(as.vector(eta_ord), taus_true), nrow = n_ord,
                 dimnames = list(NULL, ord_traits))
Y_ord[, "vigor"] <- pmin(Y_ord[, "vigor"], 3L) # category 4 never occurs here
ord_long <- two_trait_long(factor(seq_len(n_ord)), ord_traits, Y_ord)

fit_ord <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | unit, d = 1),
  data = ord_long, trait = "trait", unit = "unit",
  family = ordinal_probit(), silent = TRUE
)
extract_cutpoints(fit_ord, quiet = TRUE)[, c("trait", "cutpoint_label",
                                              "tau_estimate")]
#>                        trait cutpoint_label tau_estimate
#> ordinal_cutpoints  condition     cutpoint_2    1.0447795
#> ordinal_cutpoints1 condition     cutpoint_3    2.2778460
#> ordinal_cutpoints2     vigor     cutpoint_2    0.9586544

condition keeps all four simulated categories and returns both cutpoint_2 and cutpoint_3 (about 1.04 and 2.28, against 1.0 and 2.2). vigor was simulated the same way but its top category never appears after recoding; the fitted model has only three categories and cutpoint_3 does not exist for it. K is a property of the observed data, not a declared setting.

Truncated counts condition on Y >= 1

Reach for truncated_poisson() or truncated_nbinom2() only when the sampling process itself cannot produce or record a zero—not because zeros were removed during cleaning. The likelihood conditions on $Y \geq 1$ , renormalizing the count CDF by $1 - F(0)$ ; every response must be a positive integer. truncated_nbinom2()’s dispersion lives in fit$report$phi_truncnb2, a vector kept separate from phi_nbinom2— fitting the untruncated family on the same data and reading its phi_nbinom2 would be a different, wrong, parameter.

set.seed(20260818)
rtnbinom2 <- function(n, mu, size) {
  p0 <- dnbinom(0, mu = mu, size = size)
  qnbinom(runif(n, p0, 1), mu = mu, size = size)
}
n_trn <- 300L
trn_traits <- c("visits", "clutch")
z_trn <- rnorm(n_trn, sd = 0.4)
eta_trn <- cbind(visits = 1.2 + 0.6 * z_trn, clutch = 1.4 - 0.4 * z_trn)
phi_true_trn <- c(visits = 3, clutch = 2)
Y_trn <- matrix(
  rtnbinom2(length(eta_trn), mu = exp(eta_trn),
            size = rep(phi_true_trn, each = n_trn)),
  nrow = n_trn, dimnames = list(NULL, trn_traits)
)
trn_long <- two_trait_long(factor(seq_len(n_trn)), trn_traits, Y_trn)

fit_trn <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | unit, d = 1),
  data = trn_long, trait = "trait", unit = "unit",
  family = truncated_nbinom2(), silent = TRUE
)
data.frame(trait = trn_traits, phi_true = phi_true_trn,
           phi_hat = round(fit_trn$report$phi_truncnb2, 2))
#>         trait phi_true phi_hat
#> visits visits        3    2.28
#> clutch clutch        2    1.66

Both traits recover their dispersion in the right neighborhood (about 2.3 against 3, and 1.7 against 2), and every simulated value is at least 1—no zero ever enters the truncated likelihood, by construction of the process.

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.

set.seed(20260711)
n_mixed <- 60L
mixed_unit <- factor(seq_len(n_mixed))
z_mixed <- rnorm(n_mixed, sd = 0.45)

mixed_long <- rbind(
  data.frame(unit = mixed_unit, trait = "size", family = "continuous",
             value = rnorm(n_mixed, 0.3 + 0.7 * z_mixed)),
  data.frame(unit = mixed_unit, trait = "present", family = "presence",
             value = rbinom(n_mixed, 1, plogis(0.2 + 0.9 * z_mixed))),
  data.frame(unit = mixed_unit, trait = "abundance", family = "count",
             value = rpois(n_mixed, exp(0.4 + 0.5 * z_mixed)))
)
mixed_long$trait <- factor(mixed_long$trait)

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.

delta_gamma()’s dispersion trap sits next to Gamma()’s and points the opposite way: fit$report$phi_gamma_delta is the positive part’s coefficient of variation directly, not a shape—its near-identical neighbor phi_gamma is the inverse-square of the CV, so the two report names read alike but move oppositely for the same underlying spread. delta_lognormal()’s positive-part dispersion, sigma_lognormal_delta, is already a scale in the usual sense.

set.seed(20260818)
n_hur <- 150L
hur_traits <- c("biomass_a", "biomass_b")
z_hur <- rnorm(n_hur, sd = 0.4)
eta_hur <- cbind(biomass_a = 0.3 + 0.8 * z_hur, biomass_b = -0.2 - 0.5 * z_hur)
cv_true_hur <- c(biomass_a = 0.5, biomass_b = 0.7)
draw_delta_gamma <- function(p, mu, cv) {
  occ <- rbinom(length(p), 1, p)
  ifelse(occ == 1, rgamma(length(p), shape = 1 / cv^2, scale = mu * cv^2), 0)
}
Y_hur <- cbind(
  biomass_a = draw_delta_gamma(plogis(eta_hur[, 1]), exp(eta_hur[, 1]), cv_true_hur[1]),
  biomass_b = draw_delta_gamma(plogis(eta_hur[, 2]), exp(eta_hur[, 2]), cv_true_hur[2])
)
hur_long <- two_trait_long(factor(seq_len(n_hur)), hur_traits, Y_hur)

fit_hur <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | unit, d = 1),
  data = hur_long, trait = "trait", unit = "unit",
  family = delta_gamma(), silent = TRUE
)
data.frame(trait = hur_traits, cv_true = cv_true_hur,
           cv_hat = round(fit_hur$report$phi_gamma_delta, 2))
#>               trait cv_true cv_hat
#> biomass_a biomass_a     0.5   0.39
#> biomass_b biomass_b     0.7   0.61

The recovered CVs (about 0.39 against 0.5, and 0.61 against 0.7) are read off directly; there is no shape-to-CV conversion here, unlike Gamma() above.

Diagnose the fitted choice

Family choice is provisional until the fitted model behaves sensibly.

check_gllvmTMB(fit_nb_long)

# Exact randomized-quantile residuals are currently available for Gaussian,
# binomial, Poisson, lognormal, Gamma, NB1, NB2, Beta, betabinomial, student,
# truncated Poisson, truncated NB2, and ordinal_probit rows -- every family
# above except tweedie, delta_lognormal, delta_gamma, and multinomial.
residuals(fit_nb_long, type = "randomized_quantile", seed = 1)

# Rootograms currently support Poisson, NB1, and NB2 count rows.
predictive_check(fit_nb_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.