Choosing response families

Choose a response family from how you measured the outcome: a continuous measurement, a count, a proportion, or an ordered score. Then decide whether predictors should explain its average, its variability, or another feature such as tail weight or the probability of a zero. 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.

For current fit and reporting boundaries, use Can I fit and report this model?.

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.

Start with one measurement question

Suppose you measured growth as a continuous response and want to know whether two habitats differ in both average growth and its predictability. Because the measurements can take any real value and roughly symmetric residual variation is a sensible starting point, begin with gaussian(). This family makes the two model components easy to state: mu is expected growth and sigma is the residual standard deviation around that expectation.

The following example simulates a small data set and fits one formula for each component:

set.seed(13)
n <- 120
growth_dat <- data.frame(
  habitat = factor(rep(c("forest", "grassland"), each = n / 2))
)
growth_dat$growth <- rnorm(
  n,
  mean = ifelse(growth_dat$habitat == "grassland", 1.6, 1),
  sd = ifelse(growth_dat$habitat == "grassland", 0.9, 0.6)
)

fit_growth <- drmTMB(
  bf(growth ~ habitat, sigma ~ habitat),
  family = gaussian(),
  data = growth_dat
)

check_drm(fit_growth)
#> ok: 15; notes: 0; warnings: 0; errors: 0

The mu habitat coefficient is the fitted difference in expected growth. The sigma habitat coefficient is on the log residual-SD scale, so exponentiating it gives the fitted ratio of residual SDs:

sd_ratio <- exp(coef(fit_growth, "sigma")[["habitatgrassland"]])
round(sd_ratio, 2)
#> [1] 1.5

A ratio above 1 means that growth varies more among grassland observations after accounting for their different mean. A ratio below 1 would mean that grassland growth is more predictable. Next, run check_drm(fit_growth) as shown above and continue to When variance carries signal, Part 1 for confidence intervals, predictions, and reporting. If the response was not continuous and roughly symmetric, use the matrix below to choose a different family before adapting the formula.

At a glance

Start from the measurement process, not from the name of the distribution. The table below summarizes the available 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 can be fitted; see below for the narrower checks on structured effects
Two heavy-tailed continuous responses	biv_student()	mu1, mu2, sigma1, sigma2, shared nu, rho12	two Student-t scales and scatter/residual correlation	complete pairs and fixed effects only, with constant scales, shape, and correlation; simulation validation and intervals are unavailable; at finite nu, zero rho12 does not imply independence
Continuous with asymmetric residuals	skew_normal()	mu, sigma, nu	response standard deviation; nu is residual slant	ordinary mu random intercepts/slopes can be fitted; slope-interval evidence covers a true SD of 0.50 and at least 16 groups; sigma/nu random effects and structured effects are unavailable
Positive multiplicative response	lognormal()	mu, sigma	SD on the log-response scale	fit mu random effects or a sigma random intercept separately; scale-interval checks cover only the tested designs and missed the true value slightly too often
Two positive multiplicative responses	biv_lognormal()	mu1, mu2, sigma1, sigma2, rho12	log-response SDs and log-residual correlation	complete pairs and fixed effects only, with constant scales and correlation; use direct profile intervals within the tested sample sizes and correlations; rho12 is on the log scale, distinct from latent association eta
Positive response with mean-CV interpretation	Gamma(link = "log")	mu, sigma	coefficient of variation	fit mu random effects or a sigma random intercept separately; see below for the scale-intercept interval checks; sigma slopes and bivariate models are unavailable
Non-negative semicontinuous response with exact zeros	tweedie()	mu, sigma, nu	public scale with phi = sigma^2	ordinary mu random intercepts/slopes can be fitted; slope-interval evidence covers a true SD of 0.50 and at least 16 groups; keep nu ~ 1 and other parameters free of random effects
Continuous proportion in (0, 1)	beta_family()	mu, sigma	public scale mapped to beta precision	ordinary mu random intercepts/slopes and the tested single-term animal() effects in mu or sigma can be fitted; use zero_one_beta() for structural exact 0/1 values
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/slopes can be fitted; the separate boundary-probability random effects described below have point-estimate checks but no validated intervals
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	ordinary mu random intercepts/slopes have point-estimate checks; sigma random effects and bivariate models are unavailable
Count baseline or rate with exposure	poisson(link = "log")	mu, optional zi	no modelled sigma	ordinary and structured mu effects are available without zero inflation; consult the count tutorial before adding spatial effects to a zero-inflated model
Overdispersed count or rate	nbinom2()	mu, sigma, optional zi	extra-Poisson dispersion scale	ordinary mu effects, a grouped sigma intercept, and separate structured mu or sigma effects are available without zero inflation; see the count tutorial for the permitted formulas
Positive count, zeros absent by design	truncated_nbinom2()	mu, sigma	NB2 dispersion for the untruncated component	ordinary mu random intercepts/slopes can be fitted without hu; sigma random effects and bivariate models are unavailable
Count with a separate zero process	truncated_nbinom2() plus hu ~	mu, sigma, hu	NB2 dispersion for nonzero counts	start with fixed effects; the relatedness effect in hu described below has been checked for fitting only
Ordered categories	cumulative_logit()	mu, cutpoints	fixed latent logistic scale	ordinary mu random intercepts/slopes can be fitted; the single phylogenetic intercept has been checked for fitting only; other structured effects and scale/discrimination formulas are unavailable
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

For repeated observations within sites, individuals, or species, consider a random intercept such as (1 | site) in mu. An independent numeric slope, such as (0 + effort | site), lets the effect of effort vary among sites. The ordinary mu models with point-estimate simulation checks include Student-t, skew-normal, lognormal, Gamma, Tweedie, beta, zero-one beta, beta-binomial, binomial, Poisson, NB2, truncated NB2, and cumulative-logit models. A successful fit does not by itself establish reliable confidence intervals. Check Can I fit and report this model? for the evidence on your particular model and interval method.

For dependence described by a tree, coordinates, a pedigree, or a relatedness matrix, use phylo(), spatial(), animal(), or relmat(), respectively. Poisson and NB2 models without zero inflation allow one such term with an intercept and one numeric slope in mu. NB2 also allows these terms separately in sigma. The count tutorial shows the formulas and the narrower choices for zero-inflated and hurdle models. 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.

For non-Gaussian models with phylogenetic, spatial, pedigree, or relatedness effects, use only the family-and-term combinations listed above. Some combinations have been checked only for whether they fit and return results; others also have simulation evidence for point estimates. Neither statement automatically supports confidence intervals. Before reporting an interval, check Can I fit and report this model?. Do not combine listed terms or move them to another model component unless that combination is explicitly supported. Fixed-effect Wald intervals are available where the fixed-coefficient covariance is available; profile_targets() lists the parameters that can currently be profiled.

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_family()	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_family()	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 with simulation checks on point estimates; those checks do not validate intervals. For scale and shape effects, keep the fixed-effect formulas shown here unless using the specific structured model described above. A separate biv_student() model fits two complete responses with one shared intercept-only nu; it is checked in source tests only, without simulation validation. 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. Ordinary mu random intercepts and independent numeric slopes are available. Keep sigma and nu free of random effects; known sampling covariance, structured effects, and bivariate responses are unavailable 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 model for positive responses when multiplicative variation is biologically plausible. You may include a random intercept or independent slopes in the mean model, or a random intercept in the sigma model, but not both in the same fit. Treat intervals for the sigma random intercept cautiously: in the tested simulations, they were slightly less reliable than their nominal confidence level. Phylogenetic terms, scale slopes, known sampling covariance, and bivariate lognormal models beyond the supported fixed-effect complete-pair case are not yet available.

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. Simulation checks support profile intervals for the tested design with true random-effect SD 0.40, 12 observations per group, and at least 32 groups; results at 16 groups were borderline. These checks used maximum likelihood with a Laplace approximation. Fit it separately from the ordinary mu random-effect route; combining mu and sigma random effects is rejected. These results do not validate scale slopes, labelled covariance blocks, REML, or other study designs.

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. Ordinary mu random intercepts and independent numeric slopes are available. Keep nu ~ 1; random effects in other parameters, structured effects, bivariate responses, and separate zero-inflation or hurdle formulas are unavailable.

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_family(),
  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 have simulation checks on point estimates. You can also fit each of zoi ~ 1 + (1 | id), zoi ~ x + (0 + x | id), coi ~ 1 + (1 | id), and coi ~ x + (0 + x | id) separately. For a slope, use the same untransformed predictor in the fixed and random terms. Only point estimates have been validated. For both coi models, the population-level simulation checks used 64 groups and 50 observations per group; sparse observed atoms or weak boundary-row predictor spread can weaken individual conditional modes, so inspect both before interpreting the group-specific estimates. 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}$$

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.

To explore whether related groups have similar zero probabilities, you can fit one 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
)

This model has been checked for fitting and extracting results only; its point-estimate accuracy and intervals have not been validated. Use the fixed-effect hurdle example above for the introductory workflow. Other structured hurdle effects and slopes are unavailable.

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 can be fitted. Interval checks for independent slopes cover only specific simulation designs; consult Can I fit and report this model? before reporting an interval. 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
)

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:

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, report uncertainty for the habitat coefficient from the first table. Treat the category probabilities as fitted descriptions, not probability intervals. Do not report intervals for the individual cutpoints: drmTMB does not yet provide a public, interpretable cutpoint-interval method.

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

Keep ordinal scale fixed: sigma ~ predictors, discrimination formulas, and bivariate ordinal models are unavailable. Use the category probabilities above to describe how the fitted response changes across habitats. COM-Poisson and generalized Poisson models are also unavailable; for overdispersed counts, use the NB2 example instead.

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 for both c() and list() spellings. To begin analysing such data, fit each response separately using its appropriate family; those separate fits do not estimate their association.