Unordered categories with multinomial()

Some responses are a single choice among several categories that have no natural order: a bird’s diet guild (frugivore, insectivore, granivore), a plant’s pollination syndrome, a patch’s dominant land cover. multinomial() models such a response as a set of baseline-category log-odds and recovers how a predictor moves the odds of each category against a common reference.

This article shows when to reach for it, how to read its coefficients and predicted probabilities, and — just as important — where it stops in the current release.

When to use `multinomial()`

Match the family to how the categories relate, not just to the fact that the response is a factor.

Response	Family	Why
Unordered categories, `K >= 3`	`multinomial()`	No category is “between” the others; each gets its own contrast against a reference.
Ordered categories, `K >= 3`	`ordinal_probit()`	The order is information. A single latent axis crosses ordered thresholds; see Choose a response family.
Two categories	`binomial(link = "logit")`	One log-odds is enough. `multinomial()` on a 2-level response errors and points you here.

Two cautions before you start:

multinomial() is a response family. It is not the same as categorical(), which is a constructor for imputing a missing categorical predictor. The two are unrelated.
“Ordered-looking” numeric codes are not ordinal by themselves, and genuinely ordered categories are not multinomial. Diet guilds are unordered; abundance classes (low < medium < high) are ordered and belong to ordinal_probit().

The model: baseline-category log-odds

Pick a reference category (the first factor level by default). For a response with categories $1, \dots, K$ and reference $1$ ,

$\eta_{1} \equiv 0, \qquad \eta_{k} = \beta_{0k} + x^{\top}\beta_{k}, \quad k = 2, \dots, K,$

$\Pr(y = k) = \frac{\exp(\eta_{k})}{\sum_{j=1}^{K} \exp(\eta_{j})} \quad\text{(softmax)} .$

The reference is pinned at $\eta_1 = 0$ so the model is identified, and each non-reference category carries its own intercept and slopes. A coefficient is a contrast against the reference: $\beta_{0k}$ is the log-odds of category $k$ versus the reference at $x = 0$ , and $\beta_{k}$ is the change in that log-odds per unit of $x$ . A single categorical trait therefore occupies $K-1$ linear predictors.

Which category is the reference is only a labelling choice: the fitted probabilities and the maximised likelihood are identical whichever you pick, and only the coefficient labels change. You can set it with multinomial(baseline = ...) (shown below).

A worked example: diet guild

We simulate the diet guild of 250 bird species — frugivore, insectivore, or granivore — as a function of one scaled morphological axis, bill (larger values = deeper bills). The data-generating odds make deeper-billed species more likely to be granivores and finer-billed species more likely to be insectivores, relative to frugivores.

library(gllvmTMB)

set.seed(20260716)
n_species <- 250L
bill <- as.numeric(scale(rnorm(n_species)))

# True baseline-category log-odds, reference = frugivore.
eta <- cbind(
  frugivore   = 0,
  insectivore = 0.3 - 0.9 * bill,
  granivore   = -0.2 + 1.1 * bill
)
prob <- exp(eta - apply(eta, 1L, max))
prob <- prob / rowSums(prob)

guild <- c("frugivore", "insectivore", "granivore")[
  apply(prob, 1L, function(p) sample.int(3L, 1L, prob = p))
]

diet <- data.frame(
  species = factor(seq_len(n_species)),
  trait   = factor("diet"),
  guild   = factor(guild, levels = c("frugivore", "insectivore", "granivore")),
  bill    = bill
)

table(diet$guild)
#> 
#>   frugivore insectivore   granivore 
#>          55         116          79

This worked example has one categorical response column, so it uses the long form directly; there is no useful wide traits(...) twin for this isolated fixed-effect fit. The separate partial cross-family route uses a shared long-format mixed-family model and is demonstrated in Cross-family trait correlations, including nominal responses.

fit <- gllvmTMB(
  guild ~ 0 + trait + (0 + trait):bill,
  data   = diet,
  trait  = "trait",
  unit   = "species",
  family = multinomial(),
  silent = TRUE
)

c(convergence = fit$opt$convergence, pdHess = isTRUE(fit$sd_report$pdHess))
#> convergence      pdHess 
#>           0           1

Reading the contrasts

The fit holds $K - 1 = 2$ contrasts (insectivore and granivore, each versus the frugivore reference), each with an intercept and a bill slope. We pull them into a readable table.

sdf <- summary(fit$sd_report, "fixed")
est <- sdf[grepl("b_fix", rownames(sdf)), "Estimate"]
names(est) <- fit$X_fix_names

pick <- function(cat, slope = FALSE) {
  key <- if (slope) paste0(cat, ":bill") else paste0(cat, "$")
  unname(est[grepl(key, names(est))])
}

data.frame(
  contrast   = c("insectivore vs frugivore", "granivore vs frugivore"),
  intercept  = round(c(pick("insectivore"), pick("granivore")), 2),
  bill_slope = round(c(pick("insectivore", TRUE), pick("granivore", TRUE)), 2),
  truth_slope = c(-0.9, 1.1)
)
#>                   contrast intercept bill_slope truth_slope
#> 1 insectivore vs frugivore      0.54      -0.87        -0.9
#> 2   granivore vs frugivore     -0.06       1.05         1.1

Read the bill slopes as log-odds changes against the frugivore reference. The granivore slope is positive: deeper-billed species are more likely to be granivores than frugivores. The insectivore slope is negative: deeper-billed species are less likely to be insectivores than frugivores. Both recover the signs and rough magnitudes of the truth. An intercept is the log-odds of that guild versus frugivory at the mean bill (bill = 0).

Coefficients are contrasts, not category-level effects: a positive granivore slope does not by itself mean the probability of granivory rises everywhere, because every category shares the denominator. For probabilities, predict.

Predicted category probabilities

predict(type = "response") returns a long data frame with one row per observation and category; the est column is the fitted probability, and the probabilities sum to one within each observation.

p_hat <- predict(fit, type = "response")
head(p_hat)
#>   species trait    category       est
#> 1       1  diet   frugivore 0.2493847
#> 2       1  diet insectivore 0.5919154
#> 3       1  diet   granivore 0.1586999
#> 4       2  diet   frugivore 0.2713860
#> 5       2  diet insectivore 0.3127362
#> 6       2  diet   granivore 0.4158779

# probabilities sum to 1 within each observation
head(tapply(p_hat$est, p_hat$species, sum), 3)
#> 1 2 3 
#> 1 1 1

To see the trend, evaluate the fitted softmax across a grid of bill values directly from the estimated contrasts (the reference stays at $\eta = 0$ ):

bill_grid <- seq(-2, 2, by = 1)
eta_hat <- cbind(
  frugivore   = 0,
  insectivore = pick("insectivore") + pick("insectivore", TRUE) * bill_grid,
  granivore   = pick("granivore")   + pick("granivore",   TRUE) * bill_grid
)
p_grid <- exp(eta_hat) / rowSums(exp(eta_hat))
cbind(bill = bill_grid, round(p_grid, 2))
#>      bill frugivore insectivore granivore
#> [1,]   -2      0.09        0.90      0.01
#> [2,]   -1      0.18        0.76      0.06
#> [3,]    0      0.27        0.47      0.26
#> [4,]    1      0.23        0.16      0.61
#> [5,]    2      0.11        0.03      0.86

Fitted probability of granivory rises with bill depth while insectivory falls — the diet composition the data were built from.

Prediction is limited to the fitted data in this release: predict() with a newdata argument on a multinomial fit is not supported and fails loud rather than returning a silently wrong answer. Build the softmax grid from the coefficients (as above) when you need probabilities at new predictor values.

Changing the reference category

Any category can be the reference. Refitting with a different baseline relabels the contrasts but leaves the fit unchanged — same log-likelihood, same fitted probabilities.

fit_gran <- gllvmTMB(
  guild ~ 0 + trait + (0 + trait):bill,
  data   = diet,
  trait  = "trait",
  unit   = "species",
  family = multinomial(baseline = "granivore"),
  silent = TRUE
)

c(default_logLik  = as.numeric(logLik(fit)),
  granivore_logLik = as.numeric(logLik(fit_gran)))
#>   default_logLik granivore_logLik 
#>        -216.4474        -216.4474

The two log-likelihoods match; only whether coefficients read “versus frugivore” or “versus granivore” changes. Choose the reference that makes your contrasts easiest to interpret.

Current scope boundary

The fixed-effect route demonstrated above is validated. Two later extensions are deliberately partial, not universal categorical covariance support:

The phylogenetic route admits one phylo_latent() term and reports the standalone $(K-1) \times (K-1)$ among-category phylogenetic covariance $V$ described below. Its recovery is data-hungry.
The cross-family route permits one multinomial trait to share an ordinary latent() block with other response families. The nominal trait remains an explicit $K-1$ baseline-contrast block; extract_cross_correlations() supplies the reference-invariant partner summary. See Cross-family trait correlations, including nominal responses.
A fixed-effect-only fit still has no latent covariance, so extract_correlations() and extract_Sigma() refuse it rather than fabricate a matrix. On that cross-family route, point extraction is the default; Wald/bootstrap summaries are target-specific and uncalibrated, and nonlinear profile requests are withdrawn.
Multiple multinomial traits, augmented slopes, explicit multinomial unique()/indep(), and unlisted animal, spatial, kernel, or phylogenetic tiers remain blocked and fail loud.

Thus coefficient recovery is validated only for the fixed-effect route, while the phylogenetic and cross-family covariance routes retain their narrower partial labels. None of these rows establishes interval coverage.

The phylogenetic route: a standalone among-category table

The standalone $(K-1) \times (K-1)$ table is most natural in a phylogenetic setting, and gllvmTMB now reports it there (design 84), following the phylogenetic multinomial GLMM of Mizuno, Drobniak, Williams, Lagisz & Nakagawa (2025, Journal of Evolutionary Biology, doi:10.1093/jeb/voaf116). A phylo_latent() factor on the $K-1$ category contrasts estimates the among-category covariance $V$ — how the category liabilities coevolve across the tree — read out with extract_Sigma(fit, level = "phy", part = "shared", link_residual = "none"). That explicit call returns the fitted $V_{\mathrm{phy}}$ ; the default total extraction with link_residual = "auto" adds the fixed softmax residual $(\pi^2/6)(\mathbf I + \mathbf J)$ . A categorical trait already occupies $K-1$ latent dimensions, so a low-rank factor on those contrasts is exactly the right object; no new likelihood is needed beyond the softmax you have here.

The honest caveat is that this table is data-hungry. A single categorical draw carries little information about the $(K-1)$ -dimensional liability, so recovering $V$ needs either per-species replication or a large number of species; one observation per species is too weakly informative to pin the among-category correlation reliably (the point estimate is high-variance and can reach the $\pm 1$ boundary). Read it as recovery-oriented, and prefer replicated or large samples — the same discipline the coefficient recovery above is held to.

Two scales are in play, and it helps to name them. $V$ above is the covariance of the phylogenetic factor on the latent (log-odds) scale. To place a categorical trait on the same footing as ordinary single-scale traits you also need its observation-scale residual, and for the softmax that residual is a matrix, not a scalar: each category contrast carries $\pi^2/3$ — the same distribution-specific link variance a binomial logit uses (Nakagawa & Schielzeth 2010, Biological Reviews 85, 935–956; Nakagawa, Johnson & Schielzeth 2017, J. R. Soc. Interface 14, 20170213) — with a $\pi^2/6$ covariance between contrasts induced by their shared baseline category, i.e. $(\pi^2/6)(\mathbf{I} + \mathbf{J})$ . That is the covariance of the differenced Gumbel errors in the softmax’s random-utility representation (McFadden 1974, in Frontiers in Econometrics, ed. P. Zarembka, pp. 105–142); it reduces to the binomial’s $\pi^2/3$ when $K = 2$ , and it is the scale on which a categorical trait can eventually be made commensurable with single-scale traits.

Within those limits, multinomial() gives you an honest, identified baseline-category model for a single unordered categorical response.

When to use multinomial()