Choosing response families

This article helps users choose response families for distributional regression data: continuous outcomes, counts, proportions, percentages, positive measurements, shape-sensitive outcomes, and ordinal scores. Examples emphasize ecological, evolutionary, and environmental applications.

The implementation roadmap lives in docs/design/06-distribution-roadmap.md.

Meta-analysis is not listed as a separate family. It uses family = gaussian() with meta_V(V = V), while meta_known_V(V = V) remains a compatibility alias.

At a glance

Start from the measurement process, not from the name of the distribution. The table below summarizes the implemented first-pass families and the main quantity each sigma formula models.

Response type	First family to try	Distributional parameters	What sigma means	Main current limit
Continuous, symmetric residuals	gaussian()	mu, sigma	residual standard deviation	most mature path; non-Gaussian extensions are separate
Continuous with heavy tails	student()	mu, sigma, nu	Student-t scale, not exactly residual SD	fixed effects only
Positive multiplicative response	lognormal()	mu, sigma	SD on the log-response scale	fixed effects only
Positive response with mean-CV interpretation	Gamma(link = "log")	mu, sigma	coefficient of variation	fixed effects only
Continuous proportion in (0, 1)	beta()	mu, sigma	public scale mapped to beta precision	no 0/1 boundary mass yet
Successes out of known trials	beta_binomial()	mu, sigma	extra-binomial variation	fixed effects only
Count baseline or rate with exposure	poisson(link = "log")	mu, optional zi	no modelled sigma	ordinary non-zero-inflated mu random intercepts and independent slopes implemented; zi random effects planned
Overdispersed count or rate	nbinom2()	mu, sigma, optional zi	extra-Poisson dispersion scale	ordinary mu random intercepts and independent slopes for non-zero-inflated NB2; sigma and zi random effects planned
Positive count, zeros absent by design	truncated_nbinom2()	mu, sigma	NB2 dispersion for the untruncated component	fixed effects only
Count with a separate zero process	truncated_nbinom2() plus hu ~	mu, sigma, hu	NB2 dispersion for nonzero counts	fixed effects only; hu random effects planned
Ordered categories	cumulative_logit()	mu, cutpoints	fixed latent logistic scale	location-only; ordinal random effects and scale/discrimination planned
Two continuous Gaussian responses	c(gaussian(), gaussian())	mu1, mu2, sigma1, sigma2, rho12	residual SDs for each response	matching labelled mu1/mu2, sigma1/sigma2, and one same-response mu/sigma random-intercept pair implemented; richer covariance planned

The source-level family map in docs/design/02-family-registry.md records the links, shape or coscale slots, random-effect allowance, and test evidence behind each row.

For non-Gaussian random effects, the current fitted path is deliberately narrow. Ordinary Poisson and NB2 mu models can fit unlabelled random intercepts and independent numeric random slopes, such as (1 | site) and (0 + effort | site), when there is no zero-inflation formula. Other non-Gaussian random-effect requests are still boundary checks: non-Gaussian sigma random effects, Student-t nu random effects, zero-inflation or hurdle random effects, ordinal mixed models, structured non-Gaussian effects, and correlated or labelled count random-slope blocks need their own likelihood and recovery evidence before they become runnable syntax. Fixed-effect Wald intervals are available where the fixed coefficient covariance is available; profile targets remain limited to direct optimized parameters listed by profile_targets().

The current design keeps constructor names close to base R where possible. For example, exposure in count models is written with the standard formula term offset(log(exposure)), and Gamma models use stats::Gamma(link = "log") because gamma() already names the base R gamma function.

Why sigma is not always phi

drmTMB uses sigma as the public variability-facing scale slot. This does not mean every likelihood is written internally with a standard deviation. Instead, the package keeps the modelling grammar stable and records the family-specific conversion. The rule is simple: larger fitted sigma should mean larger modelled variability.

That rule is different from some comparator packages. For example, glmmTMB’s sigma() documentation reports a beta precision phi, where larger phi decreases the variance, and an NB2 size parameter often written theta or k, where larger values also lower the variance. In drmTMB, these precision-like parameters are internal or comparator quantities:

Family	drmTMB public scale	Internal or comparator parameter	Variability direction
gaussian()	residual SD sigma	residual SD	larger sigma means larger residual variance
Gamma(link = "log")	coefficient of variation sigma	shape 1 / sigma^2	larger sigma means larger response CV
beta()	public scale sigma	precision phi = 1 / sigma^2	larger sigma means lower precision and larger variance
beta_binomial()	extra-binomial scale sigma	precision phi = 1 / sigma^2	larger sigma means more among-trial probability variation
nbinom2()	overdispersion scale sigma	size theta = 1 / sigma^2	larger sigma means more extra-Poisson variation
student()	scale sigma; shape nu	degrees of freedom nu	larger sigma widens the core scale; larger nu lightens tails

So a sigma coefficient answers a variability question in the same direction across implemented mean-scale families. If a paper, package, or diagnostic uses phi, theta, k, shape, precision, or variance, convert explicitly before comparing numerical values.

Reporting variation

The formula name sigma is stable across families, but the reported variation is family-specific. Start with predict(fit, dpar = "mu") and sigma(fit), then transform those fitted quantities for the scale your reader needs.

Family	Variation-facing summary
gaussian()	residual SD is sigma; residual variance is sigma^2
student()	model scale is sigma; residual variance is sigma^2 * nu / (nu - 2) when nu > 2
lognormal()	log-scale variance is sigma^2; response variance is (exp(sigma^2) - 1) * exp(2 * mu + sigma^2)
Gamma(link = "log")	residual SD is mu * sigma; variance is (mu * sigma)^2
beta()	variance is mu * (1 - mu) * sigma^2 / (1 + sigma^2)
beta_binomial()	proportion variance is mu * (1 - mu) * (1 + trials * sigma^2) / (trials * (1 + sigma^2))
poisson(link = "log")	variance is mu; there is no modelled sigma
zero-inflated Poisson	variance is (1 - zi) * mu * (1 + zi * mu)
nbinom2()	variance is mu + sigma^2 * mu^2
zero-inflated nbinom2()	variance is (1 - zi) * (mu + sigma^2 * mu^2) + zi * (1 - zi) * mu^2
truncated_nbinom2() and hurdle models	use the zero-truncated or hurdle mean and variance described in the count-family section
bivariate Gaussian	marginal variances are sigma1^2 and sigma2^2; residual covariance is rho12 * sigma1 * sigma2

For Gaussian location-scale examples, sigma^2 is a residual variance. Do not carry that shortcut to Gamma, beta, count, zero-inflated, hurdle, or bivariate models without the family-specific transformation.

Implemented univariate families

Gaussian location-scale models use mu and sigma:

$$\begin{aligned} y_i \mid \mu_i, \sigma_i &\sim \operatorname{Normal}(\mu_i, \sigma_i^2),\\ \mu_i &= \beta_0 + \beta_1 \text{temperature}_i,\\ \log(\sigma_i) &= \gamma_0 + \gamma_1 \text{treatment}_i. \end{aligned}$$

drmTMB(
  bf(y ~ temperature, sigma ~ treatment),
  family = gaussian(),
  data = dat
)

Student-t location-scale-shape models add nu, the degrees-of-freedom or tail-shape parameter:

$$\begin{aligned} y_i \mid \mu_i, \sigma_i, \nu_i &\sim \operatorname{Student\text{-}t}(\mu_i, \sigma_i, \nu_i),\\ \mu_i &= \beta_0 + \beta_1 \text{temperature}_i,\\ \log(\sigma_i) &= \gamma_0 + \gamma_1 \text{treatment}_i,\\ \nu_i &= 2 + \exp(\delta_0). \end{aligned}$$

drmTMB(
  bf(y ~ temperature, sigma ~ treatment, nu ~ 1),
  family = student(),
  data = dat
)

Use this first Student-t path when a continuous response is mostly well described by a location-scale model but has heavier-tailed residuals than a Gaussian model. Here sigma is the Student-t scale parameter; when nu > 2, the residual standard deviation is sigma * sqrt(nu / (nu - 2)). Random effects, meta-analytic known covariance, phylogenetic terms, and bivariate Student-t models are planned later. The tutorial “Robust continuous responses” shows the equation, syntax, diagnostics, and a Gaussian comparison for this family.

Lognormal location-scale models are implemented for positive continuous responses such as biomass, body mass, concentration, time, and area:

$$\begin{aligned} \log(y_i) \mid \mu_i, \sigma_i &\sim \operatorname{Normal}(\mu_i, \sigma_i^2),\\ \mu_i &= \beta_0 + \beta_1 \text{habitat}_i,\\ \log(\sigma_i) &= \gamma_0 + \gamma_1 \text{treatment}_i,\\ E[y_i] &= \exp(\mu_i + \sigma_i^2 / 2). \end{aligned}$$

drmTMB(
  bf(biomass ~ habitat, sigma ~ treatment),
  family = lognormal(),
  data = dat
)

For this family, predict(fit, dpar = "mu") returns the mean of log(y), sigma(fit) returns the standard deviation of log(y), and fitted(fit) returns the arithmetic response mean. Use this path when multiplicative variation is scientifically natural and the response values are positive and finite. Random effects, known sampling covariance, phylogenetic terms, and bivariate lognormal models are planned later.

Gamma mean-CV models are implemented for positive continuous responses where relative variability is the scale target:

$$\begin{aligned} y_i \mid \mu_i, \sigma_i &\sim \operatorname{Gamma}(\text{shape}_i, \text{scale}_i),\\ \log(\mu_i) &= \beta_0 + \beta_1 \text{habitat}_i,\\ \log(\sigma_i) &= \gamma_0 + \gamma_1 \text{treatment}_i,\\ \text{shape}_i &= 1 / \sigma_i^2,\\ \text{scale}_i &= \mu_i \sigma_i^2,\\ E[y_i] &= \mu_i,\\ \operatorname{Var}(y_i) &= \mu_i^2\sigma_i^2. \end{aligned}$$

drmTMB(
  bf(biomass ~ habitat, sigma ~ treatment),
  family = Gamma(link = "log"),
  data = dat
)

For this family, predict(fit, dpar = "mu") and fitted(fit) return the response mean. sigma(fit) returns the coefficient of variation, not the residual standard deviation; the residual standard deviation is mu * sigma. Use this path for positive responses such as biomass, body mass, metabolic rate, or concentration when predictors may change relative variability. The first implementation requires stats::Gamma(link = "log"); stats::Gamma() with its default inverse link is rejected, and drmTMB does not export gamma() because base::gamma() already exists.

Beta mean-scale models are implemented for continuous responses strictly inside (0, 1), such as continuously measured cover proportions or rates that are not generated as successes out of known trials:

$$\begin{aligned} y_i \mid \mu_i, \sigma_i &\sim \operatorname{Beta}(\alpha_i, \beta_i),\\ \operatorname{logit}(\mu_i) &= \beta_0 + \beta_1 \text{habitat}_i,\\ \log(\sigma_i) &= \gamma_0 + \gamma_1 \text{treatment}_i,\\ \phi_i &= 1 / \sigma_i^2,\\ \alpha_i &= \mu_i\phi_i,\\ \beta_i &= (1 - \mu_i)\phi_i,\\ E[y_i] &= \mu_i,\\ \operatorname{Var}(y_i) &= \frac{\mu_i(1 - \mu_i)\sigma_i^2}{1 + \sigma_i^2}. \end{aligned}$$

drmTMB(
  bf(cover ~ habitat, sigma ~ treatment),
  family = beta(),
  data = dat
)

For this family, predict(fit, dpar = "mu") and fitted(fit) return the mean proportion. sigma(fit) returns the public scale parameter, not beta precision; internally phi = 1 / sigma^2, so larger sigma means more variation around the mean. Responses equal to 0 or 1 need a later zero/one-inflated beta or ordered-beta route. Percentages derived from counts should wait for beta_binomial() so the denominator remains part of the data. The tutorial Proportions and success rates works through the strict beta route and keeps boundary values separate.

Poisson mean models are implemented for count responses:

$$\begin{aligned} y_i \mid \mu_i &\sim \operatorname{Poisson}(\mu_i),\\ \log(\mu_i) &= \log(\text{trap\_nights}_i) + \beta_0 + \beta_1 \text{habitat}_i,\\ E[y_i] &= \operatorname{Var}(y_i) = \mu_i. \end{aligned}$$

drmTMB(
  bf(count ~ habitat + offset(log(trap_nights))),
  family = poisson(link = "log"),
  data = dat
)

For this family, predict(fit, dpar = "mu") and fitted(fit) return the count mean. There is no fitted sigma distributional parameter; sigma(fit) returns a fixed unit dispersion vector only for base-R method compatibility. The offset(log(trap_nights)) term is the standard R exposure form: the model estimates a rate per trap night while the expected count remains proportional to sampling effort. Use this path as a baseline count model and as a comparator for later overdispersed count families. Biological count data with extra-Poisson variation will usually need nbinom2(); COM-Poisson remains a later route for underdispersion or dispersion patterns that NB2 does not describe well.

Zero-inflated Poisson models are implemented by adding a zi formula:

$$\begin{aligned} y_i \mid \mu_i, zi_i &\sim \operatorname{ZIP}(\mu_i, zi_i),\\ \log(\mu_i) &= \log(\text{trap\_nights}_i) + \beta_0 + \beta_1 \text{habitat}_i,\\ \operatorname{logit}(zi_i) &= \gamma_0 + \gamma_1 \text{survey\_method}_i,\\ E[y_i] &= (1 - zi_i)\mu_i. \end{aligned}$$

drmTMB(
  bf(count ~ habitat + offset(log(trap_nights)), zi ~ survey_method),
  family = poisson(link = "log"),
  data = dat
)

For this model, predict(fit, dpar = "mu") returns the conditional Poisson mean, predict(fit, dpar = "zi") returns the structural-zero probability, and fitted(fit) returns (1 - zi) * mu. There is no separate zi_poisson() constructor in the current public API.

Negative-binomial 2 mean-dispersion models are implemented for overdispersed count responses:

$$\begin{aligned} y_i \mid \mu_i, \sigma_i &\sim \operatorname{NB2}(\mu_i, \text{size}_i),\\ \log(\mu_i) &= \log(\text{trap\_nights}_i) + \beta_0 + \beta_1 \text{habitat}_i,\\ \log(\sigma_i) &= \gamma_0 + \gamma_1 \text{treatment}_i,\\ \text{size}_i &= 1 / \sigma_i^2,\\ E[y_i] &= \mu_i,\\ \operatorname{Var}(y_i) &= \mu_i + \sigma_i^2\mu_i^2. \end{aligned}$$

drmTMB(
  bf(count ~ habitat + offset(log(trap_nights)), sigma ~ treatment),
  family = nbinom2(),
  data = dat
)

For this family, predict(fit, dpar = "mu") and fitted(fit) return the count mean. sigma(fit) returns the overdispersion scale in the variance equation, not a residual standard deviation. Larger sigma means greater extra-Poisson variation. The implementation uses the equivalent stats::dnbinom(mu = mu, size = 1 / sigma^2) parameterization internally. The tutorial Count abundance and extra zeros works through this sigma-to-size conversion with a soil-invertebrate example.

Zero-inflated NB2 models use the same nbinom2() family and add a zi formula for structural zeros:

$$\begin{aligned} y_i \mid \mu_i, \sigma_i, zi_i &\sim \operatorname{ZINB2}(\mu_i, \sigma_i, zi_i),\\ \log(\mu_i) &= \log(\text{trap\_nights}_i) + \beta_0 + \beta_1 \text{habitat}_i,\\ \log(\sigma_i) &= \gamma_0 + \gamma_1 \text{treatment}_i,\\ \operatorname{logit}(zi_i) &= \delta_0 + \delta_1 \text{survey\_method}_i,\\ \text{size}_i &= 1 / \sigma_i^2,\\ E[y_i] &= (1 - zi_i)\mu_i,\\ \operatorname{Var}(y_i \mid \text{count component}) &= \mu_i + \sigma_i^2\mu_i^2. \end{aligned}$$

drmTMB(
  bf(count ~ habitat + offset(log(trap_nights)), sigma ~ treatment, zi ~ survey_method),
  family = nbinom2(),
  data = dat
)

For this model, predict(fit, dpar = "mu") and sigma(fit) describe the conditional NB2 count component. predict(fit, dpar = "zi") returns the structural-zero probability, and fitted(fit) returns the unconditional response mean (1 - zi) * mu. There is no separate zi_nbinom2() constructor in the current public API. Fit this route only when the data story has a plausible structural-zero process; otherwise start with ordinary NB2 and diagnostics.

Zero-truncated NB2 models are implemented for positive counts where zeros are absent by design, such as clutch size among breeding individuals or parasite load among infected hosts:

$$\begin{aligned} y_i \mid y_i > 0, \mu_i, \sigma_i &\sim \operatorname{NB2}_{+}(\mu_i, \sigma_i),\\ \log(\mu_i) &= \beta_0 + \beta_1 \text{habitat}_i,\\ \log(\sigma_i) &= \gamma_0 + \gamma_1 \text{treatment}_i,\\ \text{size}_i &= 1 / \sigma_i^2,\\ \Pr_{+}(y_i) &= \Pr_{\operatorname{NB2}}(y_i) / \left(1 - \Pr_{\operatorname{NB2}}(0)\right),\\ E[y_i \mid y_i > 0] &= \mu_i / \left(1 - \Pr_{\operatorname{NB2}}(0)\right),\\ q_i &= 1 - \Pr_{\operatorname{NB2}}(0),\\ \operatorname{Var}(y_i \mid y_i > 0) &= \frac{\mu_i + (1 + \sigma_i^2)\mu_i^2}{q_i} - \left(\frac{\mu_i}{q_i}\right)^2. \end{aligned}$$

drmTMB(
  bf(count ~ habitat, sigma ~ treatment),
  family = truncated_nbinom2(),
  data = dat
)

For this family, predict(fit, dpar = "mu") and sigma(fit) describe the untruncated NB2 count component. fitted(fit) returns the expected observed positive count, mu / (1 - Pr_NB2(0)). Without a hu formula, the implementation rejects zeros because the sampling model is conditional on positive counts.

Hurdle NB2 models are implemented by adding hu ~ predictors to the same family route. Use this when zeros are modelled by a separate process and all nonzero counts come from a zero-truncated count distribution:

$$\begin{aligned} \operatorname{logit}(hu_i) &= \delta_0 + \delta_1 \text{survey\_method}_i,\\ \Pr(y_i = 0) &= hu_i,\\ \Pr(y_i = k > 0) &= (1 - hu_i)\Pr_{+}(k \mid \mu_i, \sigma_i),\\ E[y_i] &= (1 - hu_i)\mu_i / \left(1 - \Pr_{\operatorname{NB2}}(0)\right),\\ \operatorname{Var}(y_i) &= (1 - hu_i)\operatorname{Var}_{+i} + hu_i(1 - hu_i)E_{+i}^2. \end{aligned}$$

drmTMB(
  bf(count ~ habitat, sigma ~ treatment, hu ~ survey_method),
  family = truncated_nbinom2(),
  data = dat
)

Here hu is the hurdle-zero probability. predict(fit, dpar = "mu") still returns the untruncated NB2 component mean, while fitted(fit) returns the unconditional response mean. The design deliberately uses hu as a formula component, parallel to the implemented zi component, instead of adding a separate hurdle_nbinom2() constructor.

Use zi when the count distribution can still generate sampling zeros and you want an extra structural-zero process. Use hu when zeros are modelled separately and all nonzero counts come from a zero-truncated count distribution.

Bounded and ordered families

Percentages derived from counts should keep their denominator. Use beta_binomial() when the response is a count of successes out of known trials rather than a continuous percentage. Write the response as cbind(successes, failures), where trials_i = successes_i + failures_i; do not use cbind(successes, trials).

$$\begin{aligned} y_i \mid n_i, p_i &\sim \operatorname{Binomial}(n_i, p_i),\\ p_i \mid \mu_i, \sigma_i &\sim \operatorname{Beta}(\alpha_i, \beta_i),\\ \operatorname{logit}(\mu_i) &= \beta_0 + \beta_1 \text{habitat}_i,\\ \log(\sigma_i) &= \gamma_0 + \gamma_1 \text{treatment}_i,\\ \phi_i &= 1 / \sigma_i^2,\\ \alpha_i &= \mu_i\phi_i,\\ \beta_i &= (1 - \mu_i)\phi_i,\\ E[y_i / n_i] &= \mu_i,\\ \operatorname{Var}(y_i / n_i) &= \frac{\mu_i(1 - \mu_i)(1 + n_i\sigma_i^2)} {n_i(1 + \sigma_i^2)}. \end{aligned}$$

drmTMB(
  bf(cbind(successes, failures) ~ habitat, sigma ~ treatment),
  family = beta_binomial(),
  data = dat
)

For beta-binomial models, mu_i is the fitted success probability and sigma describes extra-binomial variation, not residual standard deviation. fitted(fit) returns probabilities; multiply by successes_i + failures_i when the scientific summary is the expected number of successes. The tutorial Proportions and success rates gives a denominator-aware seed-germination example.

Ordinal models are implemented first as fixed-effect univariate models with ordered cutpoints and a fixed latent logistic scale:

$$\begin{aligned} \Pr(y_i \le k) &= \operatorname{logit}^{-1}(\theta_k - \mu_i),\\ \mu_i &= \beta_1 \text{habitat}_i,\\ \theta_1 &< \theta_2 < \cdots < \theta_{K-1}. \end{aligned}$$

Although the R formula below is bf(score ~ habitat), the cumulative-logit fit drops the location intercept internally because a free location intercept and free cutpoints are not identifiable together.

set.seed(1)
n <- 120
habitat <- rnorm(n)
eta <- 0.8 * habitat
p_low <- plogis(-0.9 - eta)
p_medium <- plogis(0.7 - eta) - p_low
score_id <- vapply(seq_len(n), function(i) {
  sample.int(
    3,
    size = 1,
    prob = c(p_low[i], p_medium[i], 1 - p_low[i] - p_medium[i])
  )
}, integer(1))
dat <- data.frame(
  score = ordered(c("low", "medium", "high")[score_id],
    levels = c("low", "medium", "high")
  ),
  habitat = habitat
)

drmTMB(
  bf(score ~ habitat),
  family = cumulative_logit(),
  data = dat
)

Use an ordinal model for ordered scores such as disease severity, breeding condition, or habitat quality classes where category order matters but distances between categories are not numeric measurements.

For an ecology/evolution example, nest success can be recorded as ordered fledging categories. The location equation models expected reproductive success on the latent ordinal scale. fitted() returns the expected ordered-category score, sum_k k * Pr(y_i = k), which is useful for plotting the direction of a predictor but should not be treated as a measured continuous outcome. With Pr(y_i <= k) = logit^{-1}(theta_k - mu_i), larger mu_i shifts probability toward higher ordered categories.

For fitted location coefficients, confint(fit) returns Wald intervals on the latent mu coefficient scale. Internal ordered-cutpoint profile targets are visible through profile_targets(), but transformed cutpoint or category-probability intervals are not a polished user-facing surface yet.

Ordinal scale or discrimination formulas remain planned. A future extension may add sigma ~ mismatch with Pr(y_i <= k) = logit^{-1}((theta_k - mu_i) / sigma_i) or expose a direct discrimination parameter. That decision needs a separate formula-grammar and interpretation check before implementation.

Bivariate ordinal and mixed-response correlation models remain research features because the latent correlation and residual rho12 interpretation need dedicated simulation tests.

The practical next extensions are zero-one-inflated beta or ordered beta for bounded continuous responses with exact 0 or 1 values, and an ordinal scale or discrimination formula after the direction of interpretation is documented. COM-Poisson and generalized Poisson remain valuable, but they should wait until the mean-dispersion contract and comparator checks are designed.

Bivariate Gaussian families

For implemented bivariate Gaussian models, the preferred public spelling is to combine two Gaussian response families:

family = c(gaussian(), gaussian())
family = list(gaussian(), gaussian())

The all-Gaussian composed case is implemented for both c() and list() spellings and currently routes to the same likelihood as biv_gaussian().

Mixed-response composed families remain future work. For example, family = c(gaussian(), poisson()) is a planned direction for ecological examples such as body mass with fecundity counts, survival with dispersal counts, or leaf area with seed number. It is not a supported fitting path yet. The package rejects those mixed-family requests before fitting for both c() and list() spellings until the joint likelihood, residual-association parameter, prediction, simulation, intervals, and examples are designed.