---
title: "Choosing response families"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Choosing response families}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

<style>
@media (max-width: 575.98px) {
  main h1 {
    hyphens: none;
    overflow-wrap: normal;
    word-break: normal;
  }

  main table {
    display: block;
    max-width: 100%;
    overflow-x: auto;
  }
}
</style>

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?](which-scale.html) first.

For current fit and reporting boundaries, use [Can I fit and report this
model?](capability-and-limits.html).

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 remain planned |
| Two heavy-tailed continuous responses | `biv_student()` | `mu1`, `mu2`, `sigma1`, `sigma2`, shared `nu`, `rho12` | two Student-t scales and scatter/residual correlation | exact fixed-effect complete-pair source slice; at finite `nu`, zero `rho12` is uncorrelated but not independent; no random effects, scale/shape/correlation predictors, intervals, smoke, recovery, or capability claim |
| 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 |
| Two positive multiplicative responses | `biv_lognormal()` | `mu1`, `mu2`, `sigma1`, `sigma2`, `rho12` | log-response SDs and log-residual correlation | exact fixed-effect complete-pair slice; `rho12` is not raw-scale correlation or latent `eta`; direct profile is primary only within the tested fixed-effect n/rho coverage ladder, while random effects, sigma/rho predictors, and wider claims remain deferred |
| 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; exact ordinary `zoi` q1 intercept/same-raw-symbol slope plus `coi` q1 intercept/same-raw-symbol slope routes are fitted at point-fit recovery only; broader atom 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](https://glmmtmb.github.io/glmmTMB/reference/sigma.glmmTMB.html)
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}
\]

```r
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}
\]

```r
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, and additional
phylogenetic terms are planned later. A separate narrow `biv_student()` route
fits two complete responses with one shared intercept-only `nu`; it is
source-tested only and does not carry the univariate recovery evidence. 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}
\]

```r
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}
\]

```r
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 extensions beyond the
separate fixed-effect complete-pair slice 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}
\]

```r
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}
\]

```r
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}
\]

```r
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](proportion-beta-binomial.html) 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}
\]

```r
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 fitted at recovery grade. Exact ordinary
`zoi ~ 1 + (1 | id)`, same-raw-symbol `zoi ~ x + (0 + x | id)`,
`coi ~ 1 + (1 | id)`, and same-raw-symbol
`coi ~ x + (0 + x | id)` q1 routes are also fitted at point-fit recovery
only. For both `coi` routes, retained population-level recovery used `M = 64`
and 50 observations per group; sparse observed atoms or weak boundary-row
predictor spread can weaken individual conditional modes, so inspect both
before interpreting them.
Their limitations are reported as warnings: direct profiles, intervals,
coverage, joint atom effects, transformed or mismatched atom slopes, structured
effects, known covariance, denominator syntax, and bivariate bounded-response
models remain unavailable.

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}
\]

```r
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}
\]

```r
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}
\]

```r
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](count-nbinom2.html) 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}
\]

```r
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}
\]

```r
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}
\]

```r
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`:

```r
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}
\]

```r
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}
\]

```r
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](proportion-beta-binomial.html) 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.

```r
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
)

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

check_drm(fit_ordinal)
summary(fit_ordinal)
```

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.

Make that direction visible with an explicit prediction grid. The first table
is the fitted latent location shift with its ordinary Wald interval. The second
is a table of fitted category probabilities, built with the exported
`fitted_distribution()` interface:

```r
ordinal_grid <- prediction_grid(
  fit_ordinal,
  focal = "habitat",
  at = list(habitat = c(-1, 0, 1))
)

predict_parameters(
  fit_ordinal,
  newdata = ordinal_grid,
  dpar = "mu",
  conf.int = TRUE
)

ordinal_distribution <- fitted_distribution(
  fit_ordinal,
  newdata = ordinal_grid
)
ordinal_probability <- sapply(1:3, function(k) {
  ordinal_distribution$d(rep(k, nrow(ordinal_grid)))
})
colnames(ordinal_probability) <- levels(dat$score)
cbind(ordinal_grid, ordinal_probability)
```

For this fixed-effect example, the first table is the appropriate uncertainty
route for the location coefficient. The category-probability table is a fitted
description, not a category-probability interval. `profile_targets()` shows
that raw `ordinal:theta_ord:*` entries are internal optimisation coordinates;
they are not a public response-scale cutpoint estimand or a reportable cutpoint
interval on current `main`.

```r
profile_targets(fit_ordinal)[, c(
  "parm", "estimate", "scale", "profile_ready", "profile_note"
)]
```

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:

```r
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.