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. If your open question is about which scale, weights, or known variance you are modelling rather than which family, read Which scale are you modelling? first.

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 deprecated 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	ordinary mu random intercepts/slopes plus an exact recovery-grade q=1 mu ~ spatial(1 + x | ...) gate; the intercept-only spatial mu and nu ~ phylo() gates are diagnostic-only; other sigma/nu random effects and bivariate models are planned
Continuous with asymmetric residuals	skew_normal()	mu, sigma, nu	response standard deviation; nu is residual slant	ordinary mu random intercepts and independent slopes are fitted; the exact slope cell is inference-ready with caveats for true SD 0.50 and M>=16; sigma/nu random effects and structured routes are planned
Positive multiplicative response	lognormal()	mu, sigma	SD on the log-response scale	mu random intercepts/slopes and a sigma random intercept are fitted; only the exact ledger domain for the sigma intercept is inference-ready with caveats
Positive response with mean-CV interpretation	Gamma(link = "log")	mu, sigma	coefficient of variation	mu random intercepts/slopes and a sigma random intercept are fitted; its exact ordinary log-sigma intercept domain is inference-ready with caveats at true SD 0.40, n_each=12, and M >= 32 (M=16 borderline); sigma slopes and bivariate models are planned
Non-negative semicontinuous response with exact zeros	tweedie()	mu, sigma, nu	public scale with phi = sigma^2	ordinary mu random intercepts and independent slopes are fitted; the exact slope cell is inference-ready with caveats for true SD 0.50 and M>=16; nu ~ 1 remains intercept-only and other distributional random effects are planned
Continuous proportion in (0, 1)	beta()	mu, sigma	public scale mapped to beta precision	ordinary mu random intercepts/slopes and exact q=1 mu/sigma ~ animal() recovery-grade gates are fitted; exact 0/1 boundary values use zero_one_beta()
Continuous proportion in [0, 1] with structural exact 0 or 1 values	zero_one_beta()	mu, sigma, zoi, coi	public interior beta scale mapped to precision	ordinary mu random intercepts and independent slopes are fitted; the exact slope cell is inference-ready with caveats for true SD 0.50 and M>=16 under its generator-qualified evidence; sigma, zoi, and coi random effects are planned
Event indicator or successes out of known trials	stats::binomial(link = "logit")	mu	no modelled sigma; ordinary binomial sampling variation	ordinary mu random intercepts and independent slopes are fitted; use cbind(successes, failures) for trial totals
Successes out of known trials	beta_binomial()	mu, sigma	extra-binomial variation	mu random intercepts and independent slopes implemented (recovery grade); sigma random effects and bivariate models planned
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 plus q=1 phylo()/spatial()/animal()/relmat() mu intercept-plus-one-slope routes are fitted; one diagnostic-only q=1 zi ~ spatial(1 | id, coords = coords) intercept is also fitted; pure, multiple, or labelled structured slopes and other zi random effects are planned
Overdispersed count or rate	nbinom2()	mu, sigma, optional zi	extra-Poisson dispersion scale	ordinary non-zero-inflated mu random intercepts and independent slopes, grouped sigma random intercepts, q=1 phylo()/spatial()/animal()/relmat() mu intercept-plus-one-slope routes, and separate q=1 structured sigma intercept-plus-one-slope routes (recovery grade) are fitted; plain sigma random slopes, richer or labelled structured routes, and zi random effects are planned
Positive count, zeros absent by design	truncated_nbinom2()	mu, sigma	NB2 dispersion for the untruncated component	mu random intercepts and independent slopes implemented (recovery grade; not with hu); sigma random effects and bivariate models planned
Count with a separate zero process	truncated_nbinom2() plus hu ~	mu, sigma, hu	NB2 dispersion for nonzero counts	fixed-effect hurdle models plus one diagnostic-only q=1 hu ~ relmat(1 | id, K/Q = ...) intercept route are fitted; other hu random effects and structured routes are planned
Ordered categories	cumulative_logit()	mu, cutpoints	fixed latent logistic scale	ordinary mu random intercepts/slopes plus one diagnostic-only q=1 mu ~ phylo(1 | id, tree = tree) intercept are fitted; other structured ordinal effects and scale/discrimination are 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 same-response mu/sigma intercept or slope-only covariance pairs 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.

The supported non-Gaussian evidence goal is intentionally narrower than “everything non-Gaussian.” It closes the fixed-effect family evidence and the first count mixed-model slices listed here, while leaving blocked neighbours as future likelihood and simulation work.

For non-Gaussian random effects, the current fitted path is deliberately narrow. The eligible ordinary routes accept unlabelled mu random intercepts and independent numeric slopes at recovery grade, such as (1 | site) and (0 + effort | site), inside their family-specific boundaries. These routes include Student-t, skew-normal, lognormal, Gamma, Tweedie, beta, zero-one beta, beta-binomial, binomial, Poisson, NB2, truncated NB2, and cumulative-logit models. Zero-inflated and hurdle neighbours remain blocked except for the exact structured gates named below; in particular, an active hu formula does not admit ordinary count-side mu random effects. Ordinary Poisson and ordinary NB2 also have q=1 structured mu intercept-plus-one-slope routes using one of phylo(), spatial(), animal(), or relmat() on the log-mean scale. Ordinary NB2 separately has the same four q=1 structured sigma intercept-plus-one-slope routes at recovery grade. Ordinary NB2, lognormal, and Gamma also have first grouped scale routes, sigma ~ z + (1 | id), on their family-specific log-sigma scales. The lognormal and Gamma scale routes must be fitted separately from their mu random-effect routes. Tweedie still requires intercept-only nu ~ 1, and skew-normal nu remains fixed-effect residual slant. Other non-Gaussian random-effect requests are still boundary checks: Tweedie distributional random effects or predictor-dependent power, skew-normal sigma/nu random effects or structured effects, NB2/lognormal/Gamma plain sigma random slopes, combined lognormal/Gamma mu+sigma random effects, and NB2 sigma structured routes beyond the one-slope gate, Student-t nu random effects beyond the exact diagnostic-only q=1 nu ~ phylo(1 | id, tree = tree) intercept, zero-inflation random effects beyond the exact diagnostic-only q=1 Poisson zi ~ spatial() intercept, hurdle random effects beyond the exact diagnostic-only q=1 hu ~ relmat(1 | id, K/Q = ...) intercept route, ordinal mixed models beyond ordinary mu intercepts/slopes and the exact diagnostic-only q=1 mu ~ phylo(1 | id, tree = tree) intercept, other non-count spatial()/animal()/relmat() structured effects beyond the implemented recovery-grade routes (Gamma relmat() on mu, beta animal() on mu/sigma, Student-t spatial() on mu), pure, multiple, or labelled structured count slopes, simultaneous structured count effects beyond the exact crossed NB2 mu ~ spatial(1 | site, coords = coords) + relmat(1 | id, Q = Q) recovery-only gate, and correlated or labelled ordinary 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
tweedie()	public scale sigma	dispersion phi = sigma^2	larger sigma means larger semicontinuous-response variance
beta()	public scale sigma	precision phi = 1 / sigma^2	larger sigma means lower precision and larger variance
zero_one_beta()	interior beta scale sigma	precision phi = 1 / sigma^2 for the interior component	larger sigma means more variation among non-boundary proportions
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
skew_normal()	response SD sigma; slant nu	native skew-normal scale omega after moment transform	larger sigma widens residuals; positive nu gives right skew

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
skew_normal()	response residual SD is sigma; nu controls residual asymmetry
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
tweedie()	mean is mu; variance is sigma^2 * mu^nu, with 1 < nu < 2
beta()	variance is mu * (1 - mu) * sigma^2 / (1 + sigma^2)
zero_one_beta()	mean is (1 - zoi) * mu + zoi * coi; variance combines the interior beta second moment with exact boundary mass
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)). The mu predictor accepts random intercepts and independent numeric slopes at recovery grade (trust the point estimate, not the interval). Scale and shape random effects (sigma, nu), 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.

Skew-normal location-scale-shape models add nu for residual asymmetry:

$$\begin{aligned} y_i \mid \mu_i, \sigma_i, \nu_i &\sim \operatorname{SkewNormalMoment}(\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 &= \delta_0 + \delta_1 \text{habitat}_i. \end{aligned}$$

drmTMB(
  bf(y ~ temperature, sigma ~ treatment, nu ~ habitat),
  family = skew_normal(),
  data = dat
)

Use this first skew-normal path when the residual distribution is asymmetric after modelling the mean and residual scale. Here mu is the arithmetic response mean, sigma is the response standard deviation, and nu is the residual slant: positive values indicate right-skewed residuals, negative values indicate left-skewed residuals, and nu = 0 reduces to the Gaussian location-scale likelihood. Random effects, sd(group), known sampling covariance, phylogenetic or spatial terms, bivariate responses, residual rho12, and latent skew(id) syntax are planned later.

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. A mu random intercept and independent slopes run at recovery grade. An ordinary log-sigma random intercept is also fitted; only the exact tested domain in capability cell mc-0382 is inference-ready with caveats, with mildly anti-conservative rather than nominal coverage. Fit the mu and sigma random-effect routes separately; the package rejects their combination. Scale slopes, 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.

Gamma also fits one ordinary log-sigma random intercept. It is inference_ready_with_caveats only for the certified ML-Laplace profile-interval domain (true SD 0.40, n_each=12, and M >= 32; M=16 is borderline). Fit it separately from the ordinary mu random-effect route; combining mu and sigma random effects is rejected. Gamma sigma slopes, labels, REML, and any supported claim remain outside that evidence.

Tweedie mean-scale-power models are implemented for non-negative semicontinuous responses with exact zeros and positive continuous values:

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

drmTMB(
  bf(biomass ~ habitat, sigma ~ treatment, nu ~ 1),
  family = tweedie(),
  data = dat
)

For this family, fitted(fit) returns the unconditional response mean mu, including exact-zero mass. sigma(fit) returns public sigma, the square root of Tweedie dispersion phi. Use this path for biomass, cover, CPUE-like indices, or other non-negative field summaries where zeros and positive continuous values are generated by one measurement process. Random effects, predictor-dependent nu, structured effects, bivariate or mixed-response Tweedie models, zero-inflation aliases, and hurdle aliases remain planned.

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 use zero_one_beta() when those endpoints are structural outcomes. Percentages derived from counts should keep their denominator. Use stats::binomial(link = "logit") for ordinary Bernoulli/binomial event probabilities with no extra-binomial variation, or beta_binomial() when success probabilities vary beyond binomial sampling. The tutorial Proportions and success rates works through the bounded-response routes.

Zero-one beta mean-scale-boundary models are implemented for continuous proportions on [0, 1] when exact 0 or 1 values are structural outcomes rather than binomial count outcomes:

$$\begin{aligned} \Pr(y_i = 0) &= zoi_i(1 - coi_i),\\ \Pr(y_i = 1) &= zoi_i coi_i,\\ \Pr(0 < y_i < 1) &= 1 - zoi_i,\\ \operatorname{logit}(\mu_i) &= \beta_0 + \beta_1 \text{habitat}_i,\\ \log(\sigma_i) &= \gamma_0 + \gamma_1 \text{treatment}_i,\\ \operatorname{logit}(zoi_i) &= \delta_0 + \delta_1 \text{drought}_i,\\ \operatorname{logit}(coi_i) &= \kappa_0 + \kappa_1 \text{canopy}_i,\\ E[y_i] &= (1 - zoi_i)\mu_i + zoi_i coi_i. \end{aligned}$$

drmTMB(
  bf(cover ~ habitat, sigma ~ treatment, zoi ~ drought, coi ~ canopy),
  family = zero_one_beta(),
  data = dat
)

For this family, predict(fit, dpar = "mu") returns the interior beta mean. predict(fit, dpar = "zoi") returns the probability of an exact 0 or exact 1 response, and predict(fit, dpar = "coi") returns the probability that a boundary response is exactly 1. fitted(fit) returns the unconditional mean including boundary mass. Ordinary unlabelled mu random intercepts and independent numeric slopes are now fitted at recovery grade. Random effects in sigma, zoi, or coi, structured effects, known covariance, denominator syntax, and bivariate bounded-response models remain planned or blocked.

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.

One structured hurdle route is also fitted as diagnostic-only: an unlabelled q=1 relatedness intercept in hu, supplied as either a covariance matrix K or precision matrix Q:

drmTMB(
  bf(count ~ habitat, sigma ~ treatment,
    hu ~ relmat(1 | id, K = K)),
  family = truncated_nbinom2(),
  data = dat
)

Treat this as a point-estimate route only. Structured phylo(), spatial(), or animal() hurdle effects, structured hurdle slopes, labels, intervals, and coverage remain planned.

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 stats::binomial(link = "logit") when the response is an event indicator or a count of successes out of known trials with ordinary binomial sampling variation. Write the response as y01 for 0/1 data or cbind(successes, failures), where trials_i = successes_i + failures_i; do not use weights = trials, successes / trials, or cbind(successes, trials).

$$\begin{aligned} y_i \mid n_i, \mu_i &\sim \operatorname{Binomial}(n_i, \mu_i),\\ \operatorname{logit}(\mu_i) &= \beta_0 + \beta_1 \text{habitat}_i,\\ E[y_i / n_i] &= \mu_i,\\ \operatorname{Var}(y_i / n_i) &= \mu_i(1 - \mu_i) / n_i. \end{aligned}$$

drmTMB(
  bf(cbind(successes, failures) ~ habitat),
  family = stats::binomial(link = "logit"),
  data = dat
)

For plain binomial models, mu_i is the fitted event probability. fitted(fit) returns probabilities; multiply by successes_i + failures_i when the scientific summary is the expected number of successes. Ordinary mu random intercepts and independent numeric slopes are fitted first slices. Only the exact independent-slope domain in capability cell mc-0061 is inference-ready with caveats; correlated or labelled slopes remain unsupported. Factor responses, proportions with trial weights, structured effects, bivariate responses, and engine = "julia" remain unsupported.

Use beta_binomial() when the response is still counted successes out of known trials but the data show extra-binomial variation. It uses the same cbind(successes, failures) denominator syntax and adds a modelled sigma parameter for among-row probability variation.

$$\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 plain event-probability example and 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 ordinal scale or discrimination formulas after the direction of interpretation is documented, plus richer bounded-response random effects only after their fixed-effect likelihoods and diagnostics are stable. 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.