---
title: "Formula grammar"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Formula grammar}
  %\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;
    min-width: 720px;
    max-width: 100%;
    overflow-x: auto;
  }
}
</style>

This developer article documents the supported `drm_formula()` grammar,
including response blocks, distributional-parameter formulae, and validation
rules. `bf()` remains a short alias. It assumes you have already written a
first formula from a family tutorial and now need the complete grammar and
status table, not another worked example.

The design reference is `docs/design/01-formula-grammar.md`.

Fixed-effect formulas use base R's ordinary formula machinery. Transformations
and interaction expansions such as `poly(x1, 2)`, `I(x1^2)`, `x1 * x2`, and
`(x1 + x2 + x3)^2` work for implemented fixed-effect distributional-parameter
formulas. For ecological and evolutionary examples, second-order terms are
usually the practical default; third-order polynomials should be treated as
exceptional and justified by the question and data.

## Current status map

Use this table before copying syntax into an analysis. "Implemented" means the
syntax is parsed and has a fitted likelihood path for the family named in the
notes. "Reserved" means the syntax is part of the public design vocabulary but
should currently give a clear unsupported-feature error rather than a fit.
"Planned" means the syntax is shown only to explain the roadmap.

In this table, "coscale" means a model for residual correlation, currently
`rho12` in two-response Gaussian models.

If the family itself is missing from this table, the gap is not a formula
question -- see [Adding distribution families](adding-families.html) for how a
new distribution enters the package.

| Syntax | Current status | Notes |
|---|---|---|
| `drm_formula()` and `bf()` | Implemented | `drm_formula()` is the explicit constructor; `bf()` is a short alias. |
| `y ~ x1 + x2`, `sigma ~ x1` | Implemented | Univariate Gaussian location-scale model. |
| `y ~ x1 + x2`, `sigma ~ x1`, `nu ~ x2` | Implemented | Fixed-effect univariate Student-t location-scale-shape model. |
| `y ~ x1 + (1 | id) + (0 + x1 | id)`, `sigma ~ x2`, `nu ~ 1`, `family = student()` | Implemented first slice | Ordinary Student-t `mu` random intercepts and independent numeric slopes enter the identity-location predictor. A single unlabelled `spatial(1 + x | id, coords = coords)` one-slope structured effect on `mu` is recovery grade -- trust the point estimate, not the interval. The intercept-only spatial `mu` route and exact `nu ~ phylo(1 | id, tree = tree)` structured-intercept route are diagnostic grade only. Correlated slopes, labelled covariance blocks, `sigma` random effects, other `nu` random effects, other structured providers (`animal()`/`relmat()` on `mu`; `spatial()`/`animal()`/`relmat()` on `nu`), additional structured slopes and labelled q2/q4 blocks, and known covariance remain planned. |
| `bf(mu1 = y1 ~ x1, mu2 = y2 ~ x2, sigma1 = ~ 1, sigma2 = ~ 1, nu = ~ 1, rho12 = ~ 1)`, `family = biv_student()` | Implemented source slice | Exact shared-`nu` bivariate Student-t for complete finite pairs. `sigma1`/`sigma2` are scales; finite-`nu` zero `rho12` is uncorrelated but not independent. Only `mu1`/`mu2` may have predictors; random/structured effects, separate marginal `nu`, intervals, smoke, recovery, and capability claims are deferred. |
| `y ~ x1 + (1 | id) + (0 + x1 | id)`, `sigma ~ x2`, `nu ~ x3`, `family = skew_normal()` | Implemented first slice | Skew-normal ordinary `mu` random intercepts and independent numeric slopes enter the response-mean predictor. The exact independent-slope ledger cell is inference-ready with caveats for true SD 0.50 and M>=16. Public `sigma` is response SD and `nu` is residual slant. Correlated or labelled slopes, random effects outside `mu`, structured effects, bivariate routes, and broader calibration remain planned. |
| `y ~ x1 + x2`, `sigma ~ x1`, `family = lognormal()` | Implemented | Fixed-effect univariate lognormal model for positive responses; `mu` and `sigma` are on the log-response scale. |
| `y ~ x1 + (1 | id) + (0 + x1 | id)`, `sigma ~ x2`, `family = lognormal()` | Implemented first slice | Ordinary lognormal `mu` random intercepts and independent numeric slopes enter the log-response location. The exact unlabelled `(1 + x1 | id)` correlated intercept-slope block has point-fit recovery evidence. Arc 3a also admits recovery-grade unlabelled q1 `phylo(1 | id, tree = tree)` and `relmat(1 | id, K/Q = ...)` intercepts in `mu`. Labelled correlated slopes, other structured effects, known covariance, and bivariate lognormal extensions beyond the separate fixed-effect complete-pair slice remain planned. |
| `y ~ x1`, `sigma ~ x2 + (1 | id)`, `family = lognormal()` | Implemented first slice | An ordinary lognormal `sigma` random intercept enters log-`sigma`; only its exact ledger domain is inference-ready with caveats. This scale route cannot be combined with a `mu` random effect, and `sigma` slopes or labels remain unsupported. |
| `y ~ x1 + x2`, `sigma ~ x1`, `family = Gamma(link = "log")` | Implemented | Fixed-effect univariate Gamma mean-CV model for positive responses; `mu` is the response mean and `sigma` is the coefficient of variation. |
| `y ~ x1 + (1 | id) + (0 + x1 | id)`, `sigma ~ x2`, `family = Gamma(link = "log")` | Implemented first slice | Ordinary Gamma `mu` random intercepts and independent numeric slopes enter the log-mean predictor. A single unlabelled `relmat(1 | id, K = K)` or `relmat(1 + x | id, K = K)` structured intercept or one-slope on `mu` is implemented at recovery grade. Arc 3a also admits a recovery-grade unlabelled q1 `phylo(1 | id, tree = tree)` intercept. Correlated or labelled slopes, new-route slopes, other structured providers (`spatial()`/`animal()`), known covariance, and bivariate or mixed Gamma models remain planned. |
| `y ~ x1`, `sigma ~ x2 + (1 | id)`, `family = Gamma(link = "log")` | Implemented first slice | An ordinary Gamma `sigma` random intercept enters log-CV. Its ML-Laplace profile interval is inference-ready with caveats only at true SD 0.40, `n_each=12`, and `M >= 32` (`M=16` borderline). This scale route cannot be combined with a `mu` random effect; `sigma` slopes, labels, REML, and `supported` remain outside the claim. |
| `y ~ x1 + x2 + (1 | id) + (0 + x1 | id)`, `sigma ~ x1`, `nu ~ 1`, `family = tweedie()` | Implemented first slice | Tweedie model for non-negative semicontinuous responses with exact zeros; ordinary `mu` random intercepts and independent numeric slopes are fitted. The exact independent-slope ledger cell is inference-ready with caveats for true SD 0.50 and M>=16. Public `sigma` maps to `phi = sigma^2`, and `nu` is constrained to `1 < nu < 2`. Predictor-dependent `nu`, distributional random effects beyond ordinary `mu`, structured effects, zero-inflation aliases, and hurdle aliases remain planned. |
| `y ~ x1 + x2`, `sigma ~ x1`, `family = beta()` | Implemented | Fixed-effect beta mean-scale model for strict continuous proportions in `(0, 1)`; public `sigma` maps internally to `phi = 1 / sigma^2`. |
| `y ~ x1 + (1 | id) + (0 + x1 | id)`, `sigma ~ x2`, `family = beta()` | Implemented first slice | Ordinary beta `mu` random intercepts and independent numeric slopes enter the logit-mean predictor for strict `(0, 1)` responses. Row-specific structured exceptions are the recovery-grade `animal()` routes and the exact unlabelled q1 `phylo(1 | species, tree = tree)` intercept in `mu` documented below; family `sigma` remains fixed-effect-only for the phylogenetic cell. Correlated slopes, labelled covariance blocks, ordinary `sigma` random effects, exact 0/1 boundary mass, `zoi`/`coi`, `spatial()`/`relmat()`, phylogenetic slopes or scale-side structure, labelled q2/q4 or joint `mu`+`sigma` structured blocks, known covariance, beta-binomial denominators, and bivariate or mixed bounded-response models remain planned. |
| `y ~ x1 + (1 | id) + (0 + x1 | id)`, `sigma ~ x2`, `zoi ~ x3`, `coi ~ x4`, `family = zero_one_beta()` | Implemented first slice | Zero-one beta model for continuous proportions on `[0, 1]` with exact structural boundary mass and ordinary `mu` random intercepts or independent numeric slopes. The exact independent-slope ledger cell is inference-ready with caveats for true SD 0.50 and M>=16 under its generator-qualified evidence; it does not establish a general boundary-mass calibration claim. `mu` and `sigma` describe the interior beta component, `zoi` is exact-boundary probability, `coi` is the probability of an exact one conditional on the boundary, and `fitted()` includes boundary mass. Exact ordinary ML point-fit-only exceptions are one q1 intercept in `sigma` or `zoi`, one slope-only `sigma` effect `(0 + x | id)`, and one `zoi` slope when its fixed and random terms use the same raw symbol, such as `zoi ~ x + (0 + x | id)`. The two `sigma` routes now admit direct profiling, so `confint(method = "profile")` reaches their SD target; that makes an interval computable, not calibrated, and they still carry no interval, coverage, or inference claim. The `zoi` routes' direct targets are not profile-ready and they have no interval, coverage, or inference claim. `coi` random effects, other `zoi` slope forms, correlated or labelled covariance, structured effects, denominators, and bivariate or mixed bounded-response models remain planned. |
| `y01 ~ x1 + x2 + (1 | id) + (0 + x1 | id)`, `family = stats::binomial(link = "logit")` | Implemented first slice | Plain Bernoulli response model for explicit 0/1 event data. `mu` is the event probability and there is no public `sigma`; ordinary `mu` random intercepts and independent numeric slopes are fitted. Logit, probit, and complementary log-log links are available. The exact independent-slope ledger domain is inference-ready with caveats, while the exact unlabelled `(1 + x1 | id)` correlated block has point-fit recovery evidence. Factor responses, proportions with trial weights, labelled correlated slopes, structured effects beyond the exact phylogenetic intercept gate, bivariate responses, and non-logit Julia bridge routes remain unavailable. |
| `cbind(success, failure) ~ x1 + x2 + (1 | id) + (0 + x1 | id)`, `family = stats::binomial(link = "logit")` | Implemented first slice | Plain binomial response model for counted events out of known trials, with `trials = success + failure`, fixed effects plus ordinary `mu` random intercepts and independent slopes, and ordinary binomial sampling variation. Use `beta_binomial()` instead when extra-binomial variation is part of the scientific model. |
| `cbind(success, failure) ~ x1 + (1 | id) + (0 + x1 | id)`, `sigma ~ x2`, `family = beta_binomial()` | Implemented first slice | Ordinary beta-binomial `mu` random intercepts and independent numeric slopes enter the logit success-probability predictor for counted successes out of known trials. Correlated slopes, labelled covariance blocks, `sigma` random effects, `zoi`/`coi`, structured effects, known covariance, and bivariate or mixed bounded-response models remain planned. |
| `y ~ x1 + x2`, `family = poisson(link = "log")` | Implemented | Fixed-effect univariate Poisson mean model for non-negative integer counts. |
| `y ~ x1 + (1 | id) + (0 + x1 | id)`, `family = poisson(link = "log")` | Implemented first slice | Ordinary Poisson `mu` random intercepts and independent numeric slopes on the log-mean scale. The exact unlabelled `(1 + x1 | id)` correlated intercept-slope block has point-fit recovery evidence. Labelled covariance blocks, zero-inflated Poisson random effects beyond the exact diagnostic-only q=1 `zi ~ spatial(1 | id, coords = coords)` intercept, and cross-parameter covariance remain planned. |
| `y ~ x1 + offset(log(exposure))`, `family = poisson(link = "log")` | Implemented | Exposure/rate Poisson model using standard R `offset()` syntax in the `mu` formula. |
| `y ~ x1 + x2`, `zi ~ x1`, `family = poisson(link = "log")` | Implemented | Fixed-effect zero-inflated Poisson model; `mu` is the conditional count mean and `zi` is the structural-zero probability. |
| `y ~ x1 + x2`, `sigma ~ x1`, `family = nbinom2()` | Implemented | Fixed-effect univariate negative-binomial 2 model for overdispersed counts; `sigma` is an overdispersion scale. |
| `y ~ x1 + (1 | id) + (0 + x1 | id)`, `sigma ~ x2`, `family = nbinom2()` | Implemented first slice | Ordinary NB2 `mu` random intercepts and independent numeric slopes on the log-mean scale. The exact unlabelled `(1 + x1 | id)` correlated intercept-slope block has point-fit recovery evidence. Labelled covariance blocks, zero-inflated NB2 random effects outside the exact fixed-`zi` spatial structured-`mu` local-fit gate, and joint `mu`/`sigma` random effects remain planned. |
| `y ~ x1`, `sigma ~ x2 + (1 | id)`, `family = nbinom2()` | Implemented first slice | Ordinary NB2 grouped overdispersion random intercepts on the log-`sigma` scale. Exact q1 structured `sigma` intercept-plus-one-slope routes for `phylo()`, `spatial()`, `animal()`, and `relmat()` are also fitted at recovery grade. Ordinary NB2 `sigma` slopes, labelled covariance blocks, joint `mu`/`sigma` random effects, zero-inflated/truncated/hurdle scale random effects, richer structured sigma blocks, and structured-sigma intervals/coverage remain planned. |
| `y ~ x1 + offset(log(exposure))`, `sigma ~ x2`, `family = nbinom2()` | Implemented | Exposure/rate NB2 model; the offset enters the `mu` linear predictor and `sigma` remains overdispersion. |
| `y ~ x1 + x2`, `sigma ~ x1`, `zi ~ x2`, `family = nbinom2()` | Implemented | Fixed-effect zero-inflated NB2 model; `mu` and `sigma` describe the conditional count component and `zi` is structural-zero probability. |
| `count ~ x`, `sigma ~ 1 + (1 \\| group)`, `zi ~ 1`, `family = nbinom2()` | Exact point-fit-only control | Complete-response ML ZINB2 IID sigma-intercept control: fixed-effect `mu` and intercept-only `zi`; the group field enters log-`sigma`. Its direct SD target is visible but not profile-ready. Sigma predictors/slopes/labels, `mu`/`zi` random effects, missing responses, REML, intervals, and coverage remain outside this gate. |
| `count ~ x`, `sigma ~ phylo_interaction(1 \\| plant:pollinator, tree1 = plant_tree, tree2 = pollinator_tree)`, `zi ~ 1`, `family = nbinom2()` | Point-fit recovery only | Complete-response ML q1 ZINB2 route only: the structured pair field enters log-`sigma`, while `mu` and `zi` remain fixed-effect. Its direct SD target is visible but not profile-ready. A retained IID ZINB sigma-intercept control and source-bound integrated structured recovery support point-fit recovery only; no profile, interval, coverage, or inference claim follows. |
| `y ~ x1 + spatial(1 | site, coords = coords)`, `sigma ~ 1`, `zi ~ 1`, `family = nbinom2()` | Implemented diagnostic-only gate | Row-specific Q-Series v1.0 fit/extractor evidence for fixed-covariance spatial structured `mu` with fixed zero inflation. It exposes the structured SD through `summary(fit)$parameters`, conditional deviations through `ranef(fit, "spatial_mu")$terms`, and a direct `profile_targets()` row, but it does not establish point-estimate recovery; intervals, coverage, `inference_ready`, `supported`, labelled q2/q4 count covariance, zero-inflated NB2 structured `zi`, structured `sigma`, REML, AI-REML, and bridge support remain unavailable. |
| `y ~ x1 + x2`, `sigma ~ x1`, `family = truncated_nbinom2()` | Implemented | Fixed-effect zero-truncated NB2 model for positive counts; `mu` and `sigma` describe the untruncated count component and `fitted()` is the positive-count mean. |
| `y ~ x1 + (1 | id) + (0 + x1 | id)`, `sigma ~ x2`, `family = truncated_nbinom2()` | Implemented first slice | Ordinary zero-truncated NB2 `mu` random intercepts and independent numeric slopes enter the log-mean predictor while `sigma` remains fixed-effect overdispersion. One exact q1 `hu ~ relmat(K/Q)` intercept is diagnostic-only when the hurdle route is active; ordinary count-side `mu` random effects then remain blocked. Correlated slopes, labelled covariance blocks, hurdle-side random effects beyond that exact gate, `sigma` random effects, structured effects on `mu`/`sigma`, and bivariate count models remain planned. |
| `y ~ x1 + x2`, `sigma ~ x1`, `hu ~ x2`, `family = truncated_nbinom2()` | Implemented | Fixed-effect hurdle NB2 model; `hu` is the hurdle-zero probability and nonzero counts come from the zero-truncated NB2 component. A single unlabelled `hu ~ relmat(1 | id, K = K)` or `hu ~ relmat(1 | id, Q = Q)` q=1 structured intercept on the hurdle side is also implemented as diagnostic-only (local-fit only; no coverage or intervals), inherited by the `hurdle_nbinom2` alias. Structured `phylo()`/`spatial()`/`animal()`, structured slopes, labelled q2/q4 blocks, and structured `mu`/`sigma` remain planned. |
| `(1 | id)` in `mu` | Implemented | Ordinary Gaussian location random intercept. |
| `(0 + x1 | id)` in `mu` | Implemented | Independent numeric location random slope. |
| `(1 + x1 | id)` in `mu` | Implemented | Correlated ordinary intercept-slope block. |
| `(1 + x1 | p | id)` in `mu` | Implemented for naming | `p` labels the group-level covariance block for output and future matching. |
| `(1 | id)` in `sigma` | Implemented | Residual-scale random intercept in the univariate Gaussian `sigma` formula. |
| `(1 | p | id)` in both `mu` and `sigma` | Implemented | First univariate location-scale covariance slice; `p` labels one mean-scale random-intercept covariance block. |
| `sd(id) ~ x_group` | Implemented | Models the SD of an unlabelled Gaussian `mu` random intercept; predictors must be group-level. |
| `meta_V(V = V)` | Implemented | Preferred spelling for known diagonal, block-diagonal, or dense sampling covariance in a Gaussian model; bivariate Gaussian known `V` uses a complete-row `2n` by `2n` row-paired matrix. |
| `meta_known_V(V = V)` | Deprecated compatibility alias | Warns and then uses the same additive known-covariance path. New code should use `meta_V(V = V)`. |
| `mu1`, `mu2`, `sigma1`, `sigma2`, `rho12` | Implemented for fixed effects | Bivariate Gaussian location-coscale model with predictor-dependent residual correlation. |
| `(1 | p | id)` in both bivariate `mu1` and `mu2` | Implemented | First bivariate group-level covariance slice; `p` labels a `mu1`/`mu2` random-intercept covariance block. |
| `(1 | p | id)` in both bivariate `sigma1` and `sigma2` | Implemented | First bivariate residual-scale covariance slice; `p` labels a `sigma1`/`sigma2` random-intercept covariance block on the log-`sigma` scale. |
| `(0 + x | p | id)` in both bivariate `sigma1` and `sigma2` | Implemented first slice | `p` labels a `sigma1`/`sigma2` scale-slope block on the log-`sigma` scale, reported as `cor(sigma1:x,sigma2:x | p | id)`. |
| `(0 + x | p | id)` in same-response bivariate `mu1` and `sigma1`, or `mu2` and `sigma2` | Implemented first slice | Matching slope-only terms create one response-specific mean-scale-slope covariance row, such as `cor(mu1:x,sigma1:x | p | id)`. Cross-response slope labels remain planned. |
| `(1 | p | id)` in same-response bivariate `mu1` and `sigma1`, with optional independent `(1 | q | id)` in `mu2` and `sigma2` | Implemented first slice | Each labelled random-intercept pair creates its own response-specific mean-scale covariance block; residual `rho12` stays separate. |
| `(1 | p | id)` in all four bivariate `mu1`, `mu2`, `sigma1`, and `sigma2` formulas | Implemented first slice | One ordinary q=4 location-scale covariance block reports all six latent random-effect correlations. |
| `(1 + x | p | id)` in all four bivariate `mu1`, `mu2`, `sigma1`, and `sigma2` formulas | Implemented first diagnostic slice | One ordinary q8 endpoint block estimates eight endpoint SDs and reports 28 derived latent correlations. The Phase 18 smoke/recovery artifact lane reports bias, RMSE, MCSE, and interval unavailability; it is not coverage or power evidence. |
| `sd1(id) ~ x_group`, `sd2(id) ~ x_group` | Implemented | Family B direct SD formulas for labelled bivariate location random intercepts; not for residual-scale random effects. |
| `family = c(gaussian(), gaussian())` | Implemented | Routes to the bivariate Gaussian engine; mixed composed families are planned. |
| `biv_gaussian()` | Implemented legacy helper | Kept for compatibility while composed families become the public direction. |
| `mvbind(y1, y2) ~ x1` | Implemented | Shorthand for identical bivariate location formulas; explicit `mu1`/`mu2` remains preferred for different predictors. |
| `phylo(1 | species, tree = tree)` in `mu` and/or `sigma` | Implemented | Intercept-only univariate Gaussian phylogenetic location and residual-scale effect; matching `mu`/`sigma` terms estimate one latent structured correlation. Requires an ultrametric tree with branch lengths. |
| `y ~ x + phylo(1 | species, tree = tree)`, `sigma ~ z`, `family = beta()` | Implemented exact q1 recovery slice | One unlabelled intercept-only phylogenetic location effect in Beta `mu`; family `sigma` remains fixed-effect-only and controls `phi = sigma^(-2)`. Point recovery passed only in the exact tested `g = 1024, m = 4` regime. This is not `g >= 1024`; slopes, labels/q2+, phylogeny in family `sigma`, direct latent-`sd()` regression, REML, intervals, and coverage remain deferred. |
| `animal(1 | id, pedigree = pedigree)`, `animal(1 | id, A = A)`, or `animal(1 | id, Ainv = Ainv)` in `mu` and/or `sigma` | Implemented first slice | Univariate Gaussian animal-model random intercept from a dense pedigree-derived additive relationship matrix, precomputed additive relatedness, or inverse-relatedness matrix. Matching `mu`/`sigma` terms estimate one latent structured correlation. One numeric `mu` slope is fitted; the exact A-matrix q1 `sigma` one-slope route is also fitted and inference-ready with caveats under raw log-SD Wald-z intervals. Pedigree/Ainv bridge marshalling, sparse large-pedigree construction, additional multiple or labelled structured-slope layouts outside the exact fitted ledger cells, structured slope correlations, predictor-dependent `corpair()`, and direct-SD grammar remain planned. |
| matching labelled `animal(1 | p | id, pedigree = pedigree)`, `animal(1 | p | id, A = A)`, or `animal(1 | p | id, Ainv = Ainv)` in bivariate `mu1` and `mu2` | Implemented first q=2 slice | Estimates two latent animal-model location scales and one animal mean-mean correlation from the same pedigree-derived or precomputed matrix. Individual marginal SDs also include `sqrt(A[i, i])`. |
| matching labelled `animal(1 | p | id, pedigree = pedigree)`, `animal(1 | p | id, A = A)`, or `animal(1 | p | id, Ainv = Ainv)` in all four bivariate `mu1`, `mu2`, `sigma1`, and `sigma2` formulas | Implemented first q=4 slice | Estimates four latent animal-model endpoint scales and six constant latent location-location, location-scale, and scale-scale correlations from the same matrix. Individual marginal SDs also include `sqrt(A[i, i])`. |
| `relmat(1 | id, K = K)` or `relmat(1 | id, Q = Q)` in `mu` and/or `sigma` | Implemented first slice | Lower-level univariate Gaussian random intercept from a user-supplied latent relatedness covariance or precision matrix. Matching `mu`/`sigma` terms estimate one latent structured correlation. Keep observation-level known sampling covariance in `meta_V(V = V)`. |
| matching labelled `relmat(1 | p | id, K = K)` or `relmat(1 | p | id, Q = Q)` in bivariate `mu1` and `mu2` | Implemented first q=2 ML slice; exact supplied-`K` intercept-only REML cell at `point_fit_recovery` | Estimates two latent relatedness location scales and one mean-mean correlation from the same user-supplied matrix. Native `REML = TRUE` requires `biv_gaussian()`, the identical named `K`, group ordering, and label in both location formulas, intercept-only `sigma1`, `sigma2`, and `rho12`, complete pairs, unit weights, and no other random, scale-side, `meta_V()`, direct-SD, or `corpair()` layer. `Q`, slopes, q4+, intervals, and coverage remain outside this REML cell. Level-specific marginal SDs also include the known diagonal multiplier. |
| matching labelled `relmat(1 | p | id, K = K)` or `relmat(1 | p | id, Q = Q)` in all four bivariate `mu1`, `mu2`, `sigma1`, and `sigma2` formulas | Implemented first q=4 slice | Estimates four latent relatedness endpoint scales and six constant latent location-location, location-scale, and scale-scale correlations from the same matrix. Level-specific marginal SDs also include the known diagonal multiplier. |
| matching `phylo(1 | species, tree = tree)` in bivariate `mu1` and `mu2` | Implemented first slice | Estimates two phylogenetic location SDs and one phylogenetic mean-mean correlation. |
| labelled `phylo(1 | p | species, tree = tree)` in matching bivariate `mu1` and `mu2` | Implemented | The label is preserved in SD, correlation, `corpairs()`, and profile-target names for the phylogenetic mean-mean path. |
| labelled `phylo(1 | p | species, tree = tree)` in all four bivariate `mu1`, `mu2`, `sigma1`, and `sigma2` formulas | Implemented first slice | One constant q=4 phylogenetic location-scale block estimates four endpoint SDs and six latent phylogenetic correlations. Partial, unlabelled, mismatched, and slope forms remain rejected. |
| `count ~ x + provider(1 + x | p | group, input = object)`, `family = poisson(link = "log")` | Exact labelled q2 point-fit slice | Ordinary non-zero-inflated Poisson admits exactly one labelled intercept-plus-one-numeric-slope `mu` block for `phylo(tree = tree)`, `spatial(coords = coords)`, `animal(Ainv = Ainv)`, or `relmat(Q = Q)`. It estimates two latent SDs and their intercept--slope correlation. Its three direct targets (both SDs and the correlation) now admit direct profiling, so `confint(method = "profile")` reaches them; that makes an interval computable, not calibrated, and the retained evidence remains point-fit-only with no interval or coverage claim. Pure structured slopes, multiple structured slopes, q4+, zero-inflated structure, simultaneous providers, ordinary count random effects, alternate provider inputs, intervals, and coverage remain outside this slice. |
| `count ~ x + phylo(1 | species, tree = tree)` or `count ~ x + relmat(1 | id, Q = Q)`, `sigma ~ z`, `family = nbinom2()` | Implemented first slice | Ordinary non-zero-inflated NB2 q=1 structured `mu` intercept on the log-mean scale for one of `phylo()`, `spatial()`, `animal()`, or `relmat()`, while `sigma` remains fixed-effect overdispersion unless using a separate fitted sigma route; unlabelled intercept-plus-one-slope `mu` terms are also fitted for those four providers. The separate exact labelled NB2 q2 exception is `phylo(1 + x | p | species, tree = tree)` only, at point-fit-only status. Exact q1 structured `sigma` intercept-plus-one-slope routes for those four providers are recovery grade. One exact crossed `mu ~ spatial(1 | site, coords = coords) + relmat(1 | id, Q = Q)` route is recovery-only. Pure, multiple, or other labelled structured count slopes, q4+ count blocks, zero-inflated NB2 structure outside the exact local-fit gate, richer or labelled structured sigma, structured-sigma intervals/coverage, simultaneous structured types beyond that exact crossed gate, and joint `mu`/`sigma` random effects remain planned. |
| `count ~ x + phylo_interaction(1 | partner1:partner2, tree1 = tree1, tree2 = tree2)` or Gaussian `y ~ x + phylo_interaction(...)` | Implemented first slice | Single q=1 location random intercept for a two-partner phylogenetic interaction in univariate Gaussian `mu` and ordinary Poisson/NB2 `mu`. Internally this builds a sparse Kronecker precision from the two augmented phylogenetic precisions. `relmat(1 | pair, Q = Q_pair)` remains the lower-level escape hatch. Binary/Bernoulli incidence models and additive models that combine partner main phylogenies plus `phylo_interaction()` remain planned. |
| `y ~ x1 + x2 + (1 | id) + (0 + x1 | id)`, `family = cumulative_logit()` | Implemented first slice | Univariate ordinal model for ordered scores with cutpoints; ordinary `mu` random intercepts and independent numeric slopes are fitted. The public `mc-0227` slope uses ML-Laplace and is point-fit-recovery only. Its retained M>=80 O3 campaign concerns a package-private estimator and grants no public interval or reporting permission. A single unlabelled `phylo(1 | species, tree = tree)` structured `mu` intercept is diagnostic-only (fit/extractor feasibility, not point-estimate recovery or coverage). Correlated or labelled slopes, other structured routes (`spatial()`/`animal()`/`relmat()`), and scale/discrimination formulas are planned. |
| `cbind(successes, failures) ~ x1`, `family = beta_binomial()` | Implemented | Fixed-effect denominator-aware model for success counts with known trial totals; `sigma` is extra-binomial variation. |
| `phylo(1 + x1 | species, tree = tree)` | Implemented first one-slope slice | Univariate Gaussian `mu` path with independent phylogenetic intercept and slope fields; ordinary Poisson and NB2 also fit one unlabelled q1 intercept-plus-one-slope phylogenetic `mu` route at recovery/source-test grade. Pure, multiple, or labelled count slopes and labelled structured slope blocks such as `phylo(1 + x | p | species, tree = tree)` remain rejected until their covariance and recovery contracts are designed. |
| `animal(1 + x | id, pedigree = ped)` | Implemented first one-slope slice | Univariate Gaussian `mu` path with independent animal-model intercept and slope fields; additional multiple-slope layouts and slope correlations remain planned outside the exact fitted bivariate ledger cells. A labelled block in this univariate q1 route is rejected; use only the exact labelled bivariate layouts recorded as fitted in the ledger. |
| `spatial(1 | site, coords = coords)` | Implemented first slice | Univariate Gaussian `mu` spatial random intercept using a fixed dense coordinate covariance foundation, separate from the mesh/SPDE route. |
| `spatial(1 | site, mesh = mesh)` with `sigma ~ 1` | Exact fixed-domain `n = 128, 256` designs at `point_fit_recovery` | Univariate Gaussian ML `mu` intercept using a fixed-kappa projected mesh field `A_st %*% omega`. The reported GMRF field scale is not a range or uniform marginal SD. The retained `n = 64` recovery rung failed; intervals and coverage remain unclaimed. |
| `spatial(1 + x | site, coords = coords)` | Implemented first one-slope slice | Univariate Gaussian `mu` spatial random intercept plus one numeric spatial slope. The two coordinate fields are independent and have separate SDs. A labelled block in this univariate q1 route is rejected; use only the exact labelled bivariate layouts recorded as fitted in the ledger. |
| `spatial(1 | p | site, coords = coords)` in matching bivariate `mu1` and `mu2` formulas | Implemented q=2 slice; exact intercept-only native REML cell at `inference_ready_with_caveats` for tested M/H ring configurations | Bivariate Gaussian location covariance for two coordinate-spatial random intercept fields. `corpairs(level = "spatial")` reports the spatial mean-mean row separately from residual `rho12`. Native `REML = TRUE` requires matching labels; intercept-only `sigma1`, `sigma2`, and `rho12`; complete response pairs; unit weights; no known `meta_V()` covariance; and no additional ordinary random, direct-SD, or `corpair()` layer. Joint endpoint-profile calibration passes at exact M = 36 sites x 3 observations and H = 36 x 8 under the baseline ring geometry; L = 12 x 3 fails. Slopes, range estimation, other geometries or information configurations, and the `supported` tier remain planned. |
| `spatial(1 | p | site, coords = coords)` in all four bivariate `mu1`, `mu2`, `sigma1`, and `sigma2` formulas | Implemented first q=4 slice | One constant coordinate-spatial location-scale block estimates four endpoint SDs and six latent spatial correlations. Partial, unlabelled, mismatched, and slope forms remain rejected. |
| `relmat(1 + x | id, K = K)` | Implemented first one-slope slice | Univariate Gaussian `mu` path with independent lower-level relatedness intercept and slope fields; additional multiple-slope layouts and slope correlations remain planned outside the exact fitted bivariate ledger cells. A labelled block in this univariate q1 route is rejected; use only the exact labelled bivariate layouts recorded as fitted in the ledger. |
| `corpair(id, level = "group", block = "p", from = "mu1", to = "mu2") ~ x_group` | Implemented | Predictor-dependent ordinary q=2 location-location latent random-effect correlation regression for matching labelled `mu1`/`mu2` random intercepts. Predictors must be constant within `id`. |
| `corpair(species, level = "phylogenetic", block = "p", from = "mu1", to = "mu2") ~ ecology` | Implemented | Predictor-dependent phylogenetic q=2 location-location latent random-effect correlation regression for matching labelled `mu1`/`mu2` `phylo()` terms. Predictors must be constant within `species`. Location-scale, scale-scale, q=4, and spatial `corpair()` regressions remain planned. |
| Matching slope-only `(0 + x | p | id)` in bivariate `mu1` and `mu2` | Implemented first bivariate slope slice | Estimates the ordinary slope1-slope2 plasticity-syndrome correlation without intercept-slope correlations; broader bivariate slope blocks remain planned. |
| Matching one-slope intercept-plus-slope `(1 + x | p | id)` blocks in bivariate `mu1` and `mu2` | Implemented first q=4 location slice | Estimates two location-intercept SDs, two location-slope SDs, and six group-level latent correlations; the SDs are direct profile targets and the correlations are derived-unavailable interval rows. |
| Matching two-slope `(1 + x + z | p | id)` blocks in bivariate `mu1` and `mu2` | Implemented q=6 location slice with smoke artifact routing | Estimates six location SDs and 15 group-level latent correlations; the SDs are direct profile targets and the correlations are derived-unavailable interval rows. |
| p8 endpoint blocks, q8 variants beyond the first matching all-four ordinary endpoint slice, labelled or cross-formula residual-scale slope blocks, broader bivariate residual-scale blocks beyond the matching q2 `sigma1`/`sigma2` scale-slope slice, non-Gaussian cross-parameter covariance, predictor-dependent q=4 phylogenetic/spatial correlations, or `rho12` random effects | Planned | These need a larger covariance parameterization, simulation, and naming checks before fitting. Unlabelled ordinary univariate correlated `sigma` blocks are already fitted. Intercept-slope `corpair()` rows are distant-future; broader slope-correlation regressions need coefficient-aware syntax. |
| `corpairs(fit)` | Implemented for current fitted correlations | Reports residual `rho12`, ordinary and bivariate group-level covariance rows including ordinary q=4 blocks, bivariate phylogenetic mean-mean or q=4 rows, bivariate coordinate-spatial q=2 or constant q=4 rows, and the animal/`relmat()` q=2 or constant q=4 rows. |

## Offsets

Every univariate family accepts a standard R `offset()` term in its `mu`
formula. The offset is a known constant added to the location linear predictor,
so what it means depends on the family's link.

For **log-link** families — Poisson, `nbinom2()` and its zero-inflated path,
`Gamma(link = "log")`, and `tweedie()` — an offset is the familiar exposure or
effort model. Use it when counts or amounts accumulate over unequal observation
windows; the coefficients then describe a rate:

```r
drmTMB(bf(mu = count ~ habitat + offset(log(trap_nights))),
       data = surveys, family = poisson(link = "log"))
```

For **identity-link** families — Gaussian, `student()`, `lognormal()`, and
`skew_normal()` — an offset is a known additive shift of the mean. Use it when
part of the response is already known and should not cost a parameter, such as
a baseline you want to subtract.

For **logit-link** families — Bernoulli/binomial, `beta_binomial()`, `beta()`,
`zero_one_beta()`, and `cumulative_logit()` — an offset is a known shift in
log-odds. This is a **calibration** term, not an exposure: writing
`offset(log(effort))` here does not give you a rate model, and the value you
supply is on the log-odds scale. Two family-specific details matter: a
`zero_one_beta()` offset moves only the interior beta component and leaves the
`zoi` and `coi` boundary probabilities alone, and a `cumulative_logit()` offset
shifts the latent location relative to fixed cutpoints rather than an intercept.

What still rejects an offset, and why:

| Route | Reason |
| --- | --- |
| `truncated_nbinom2()` and its hurdle path | Both renormalise the observed mean over a restricted support, so doubling exposure would not double the reported mean the way it does for plain Poisson or NB2. Extending them needs a decision about whether the offset targets the latent untruncated rate or the observed mean. |
| Every bivariate family | There is no per-response offset contract yet. |
| Gaussian sufficient-statistic aggregation | Aggregation reduces rows to group summaries, so a row-level offset has no aggregated counterpart; fit the row-level model instead. |
| `sigma`, `zi`, `hu`, `zoi`, `coi`, `nu`, `rho12` | Offsets are restricted to the location formula. |

If you need a rate model for a route in that first row, put the exposure in the
formula as a log-scale covariate and interpret the estimated coefficient instead
of fixing it at one.

Current implemented univariate Gaussian random-effect forms include:

```r
drm_formula(y ~ x + (1 | id), sigma ~ z)
drm_formula(y ~ x + (0 + x | id), sigma ~ z)
drm_formula(y ~ x + (1 + x | id), sigma ~ z)
drm_formula(y ~ x + (1 + x | p | id), sigma ~ z)
drm_formula(y ~ x + (1 | id), sigma ~ z + (1 | id) + (0 + w | id))
drm_formula(y ~ x, sigma ~ z + (1 + w | id))
drm_formula(y ~ x + (1 | p | id), sigma ~ z + (1 | p | id))
drm_formula(y ~ x + (1 | id), sigma ~ z, sd(id) ~ x_group)
drm_formula(y ~ x + phylo(1 | species, tree = tree), sigma ~ z)
drm_formula(y ~ x + animal(1 | id, Ainv = Ainv), sigma ~ z)
drm_formula(y ~ x + relmat(1 | line, Q = Q), sigma ~ z)
drm_formula(
  y ~ x + (1 | id) + (1 | site),
  sigma ~ z,
  sd(id) ~ x_group,
  sd(site) ~ site_type
)
```

Current implemented Student-t form:

```r
drm_formula(y ~ x, sigma ~ z, nu ~ 1)
drm_formula(y ~ x, sigma ~ z, nu ~ x)
drm_formula(y ~ x + (1 | id) + (0 + x | id), sigma ~ z, nu ~ 1)
```

Student-t ordinary `mu` random intercepts and independent numeric slopes are a
first location slice, joined by a `spatial(1 | id, coords = coords)` structured
`mu` intercept (diagnostic grade only; it confirms fit/extractor feasibility but
does not establish point-estimate recovery)
and a `phylo(1 | id, tree = tree)` structured `nu` intercept (diagnostic grade
only). Correlated Student-t slopes, `sigma` random effects, other `nu` random
effects, known sampling covariance, other structured providers on `mu`/`nu`,
and richer bivariate Student-t models are later phases. The separate
`biv_student()` source slice admits fixed-effect `mu1`/`mu2` with intercept-only
Student-t scales, one shared `nu`, and intercept-only `rho12`; it does not
inherit the univariate recovery or interval evidence.

Future skew-normal and skew-t formulas should also start with fixed-effect
residual-shape terms such as `nu ~ x`. A later latent-effect spelling such as
`skew(id) ~ x_group` would ask a different question about the distribution of
group effects and should not be taught as an alias for residual `nu`.

Current implemented lognormal form:

```r
drm_formula(biomass ~ habitat, sigma ~ treatment)
```

paired with:

```r
family = lognormal()
```

The response must be positive and finite. `mu` and `sigma` are defined on the
log-response scale, and `fitted()` returns the arithmetic response mean
`exp(mu + sigma^2 / 2)`.

Current implemented Gamma form:

```r
drm_formula(biomass ~ habitat, sigma ~ treatment)
```

paired with:

```r
family = Gamma(link = "log")
```

Here `mu` is the response mean and `sigma` is the coefficient of variation.
Gamma `mu` random intercepts and independent numeric slopes are an implemented
first slice at recovery grade. A separate ordinary `sigma` random intercept is
also fitted at recovery grade; it cannot be combined with a `mu` random effect.
A single unlabelled `relmat(1 | id, K = K)` intercept
or one-slope structured effect on `mu` is fitted at recovery grade. Correlated
or labelled slopes, `sigma` slopes, known sampling covariance,
`phylo()`/`spatial()`/`animal()` structured terms, and bivariate or mixed Gamma
models are later phases.

The middle `p` is a group-level covariance-block label. It is not residual
`rho12`, and reserved distributional parameter names such as `rho12` should not
be used as labels. Matching `(1 | p | id)` terms in univariate Gaussian `mu`
and `sigma` fit the first mean-scale covariance block; larger cross-formula
blocks with slopes remain future work.

Current implemented bivariate form:

```r
drm_formula(
  mu1 = y1 ~ x1 + x2,
  mu2 = y2 ~ x1,
  sigma1 = ~ x1 + x2,
  sigma2 = ~ x1,
  rho12 = ~ x1 + x2
)
```

When both responses share the same location predictors, the shorter form is
also implemented:

```r
drm_formula(
  mvbind(y1, y2) ~ x1 + x2,
  sigma1 = ~ x1 + x2,
  sigma2 = ~ x1,
  rho12 = ~ x1 + x2
)
```

This expands internally to `mu1 = y1 ~ x1 + x2` and
`mu2 = y2 ~ x1 + x2`. Use the explicit `mu1` and `mu2` formulas when the two
responses need different location predictors.

The first fitted structured-effect marker examples are phylogenetic location
effects:

```r
drm_formula(y ~ x + phylo(1 | species, tree = tree), sigma ~ z)
drm_formula(y ~ x + phylo(1 + x | species, tree = tree), sigma ~ z)
```

`phylo()` uses an ultrametric branch-length tree and adds tree-structured
Gaussian random intercept and one-slope fields to the `mu` predictor. Matching
intercept-only `phylo()` terms in bivariate `mu1` and `mu2` formulas fit the first
phylogenetic mean-mean correlation slice:

```r
drm_formula(
  mu1 = y1 ~ x + phylo(1 | species, tree = tree),
  mu2 = y2 ~ x + phylo(1 | species, tree = tree),
  sigma1 = ~ 1,
  sigma2 = ~ 1,
  rho12 = ~ 1
)
```

Matching labelled `phylo()` terms across all four bivariate Gaussian endpoints
fit the constant q=4 phylogenetic location-scale block. A single shared label
estimates one full q4 tree block; using one label for `mu1`/`mu2` and another
for `sigma1`/`sigma2` estimates the block-diagonal fallback:

```r
drm_formula(
  mu1 = y1 ~ x + phylo(1 | pl | species, tree = tree),
  mu2 = y2 ~ x + phylo(1 | pl | species, tree = tree),
  sigma1 = ~ z + phylo(1 | ps | species, tree = tree),
  sigma2 = ~ z + phylo(1 | ps | species, tree = tree),
  rho12 = ~ 1
)
```

Additional multiple or labelled phylogenetic-slope layouts outside the exact
fitted ledger cells, richer residual-scale structured slopes beyond the exact
q1 routes, partial or mismatched bivariate phylogenetic scale terms,
structured effects in `rho12`, multiple animal-model slopes, multiple
lower-level `relmat()` slopes, multiple spatial slopes, and mesh-based spatial
terms are still planned:

```r
# planned; drmTMB() will currently reject these
drm_formula(y ~ x + phylo(1 + x1 + x2 | species, tree = tree), sigma ~ z)
drm_formula(y ~ x + phylo(1 + x | p | species, tree = tree), sigma ~ z)
drm_formula(y ~ x + animal(1 + x1 + x2 | id, pedigree = ped), sigma ~ z)
drm_formula(y ~ x + animal(1 + x | p | id, pedigree = ped), sigma ~ z)
drm_formula(y ~ x + relmat(1 + x1 + x2 | id, K = K), sigma ~ z)
drm_formula(y ~ x + relmat(1 + x | p | id, K = K), sigma ~ z)
drm_formula(y ~ x + spatial(1 + x | p | site, coords = coords), sigma ~ z)
drm_formula(y ~ x + spatial(1 | site, mesh = mesh), sigma ~ z)
```

The first fitted `spatial()` paths use coordinates directly for univariate
Gaussian `mu` coefficient fields:

```r
drmTMB(
  bf(y ~ x + spatial(1 | site, coords = coords), sigma ~ z),
  family = gaussian(),
  data = dat
)

drmTMB(
  bf(y ~ x + spatial(1 + x | site, coords = coords), sigma ~ z),
  family = gaussian(),
  data = dat
)
```

`coords` can have one row per site or one row per observation, but coordinates
must be constant within each site. A precomputed `mesh` is fitted only for the
fixed-kappa Gaussian intercept documented above. Mesh is not the scientific
sampling level; it is the computational scaffold for the projected SPDE/GMRF
field. These structured effects are separate from residual `rho12`.

For structured slopes, the fitted coordinate-spatial path now supports one
numeric `mu` slope as independent intercept and slope fields with separate SDs.
Additional structured slopes remain the advanced path after more recovery
evidence.
Multiple random factors should remain separate additive blocks. The first slope
models should not estimate predictor-dependent intercept-slope `corpair()` rows.
The first fitted bivariate targets are the slope1-slope2 correlation and the
q4/q6 location blocks with smoke artifact routing for the same covariates across
responses.

The first bivariate group-level covariance syntax is explicit about the
distinction between group-level covariance labels and residual `rho12`:

```r
drm_formula(
  mu1 = y1 ~ x1 + x2 + (1 | p | ID),
  mu2 = y2 ~ x1      + (1 | p | ID),
  sigma1 = ~ x1 + x2 + (1 | q | ID),
  sigma2 = ~ x1      + (1 | q | ID),
  rho12 = ~ x1 + x2
)
```

The shared `p` label fits a group-level random-intercept covariance block for
the two response means. The shared `q` label fits a separate scale-scale block
for the two residual-scale random intercepts on the log-`sigma` scale. The
same-response mean-scale bridge is also implemented for one or both responses:
a shared label in `mu1` and `sigma1` gives a response-1 mean-scale row, while a
different shared label in `mu2` and `sigma2` gives a separate response-2
mean-scale row. Reusing the same label and group in all four bivariate formulas
fits the first ordinary q=4 location-scale block:

```r
drm_formula(
  mu1 = y1 ~ x + (1 | p | ID),
  mu2 = y2 ~ x + (1 | p | ID),
  sigma1 = ~ z + (1 | p | ID),
  sigma2 = ~ z + (1 | p | ID),
  rho12 = ~ w
)
```

That q=4 block reports one location-location row, four location-scale rows,
and one scale-scale row in `corpairs()`. Fitted rows still use the older
`mean-mean` and `mean-scale` class names, but `corpairs()` accepts
`location-location` and `location-scale` as filter aliases to match the
reserved `corpair()` formula wording. The first fitted `corpair()` route is the
ordinary q=2 location-location model:
`corpair(ID, level = "group", block = "p", from = "mu1", to = "mu2") ~ w`.
The predictor `w` must be constant within `ID`, and `corpairs()` reports the
mean, range, and number of fitted group-level correlations. Later
location-scale routes should also be endpoint-specific, for example
`corpair(ID, level = "group", block = "p", from = "mu1", to = "sigma2") ~ w`,
because the four location-scale pairs in q=4 models do not have identical
interpretation. For phylogenetic life-history trade-off models, the same
grammar now also covers the fitted q=2 location-location route:
`corpair(species, level = "phylogenetic", block = "p", from = "mu1", to = "mu2") ~ ecology`.
The phylogenetic formula has an extra design gate: because the tree couples all
species, `drmTMB` needs a positive-definite covariance contract for the full
species block before a species-level predictor can change the latent
phylogenetic correlation. The selected q=2 contract uses two independent unit
phylogenetic fields and species-specific loadings, so variable `rho_l` is a
nonstationary phylogenetic covariance model rather than a residual-style
correlation row. It covers `mu1`-`mu2` only; phylogenetic location-scale and
scale-scale correlation regressions remain q=4 extensions. Multiple-slope
bivariate location blocks beyond the smoke-routed ordinary q4/q6 lanes,
bivariate spatial slope blocks, predictor-dependent q=4 phylogenetic or spatial
correlations, p8 endpoint slope covariance, q8 variants beyond the first
ordinary diagnostic lane, and `rho12` random effects remain planned.
Double-hierarchical bivariate fits should continue to report correlation pairs
in a long table, not just as unlabeled matrices. The design note
`docs/design/20-coscale-correlation-pairs.md` records the planned labels:
correlation level, group, covariance block, distributional parameters,
responses, and random-effect coefficients. The current `corpairs(fit)` helper
already uses this table shape for residual `rho12`, ordinary univariate
Gaussian `mu` random-effect correlations, and the implemented bivariate
mean-mean, scale-scale, response-specific mean-scale, ordinary q=4, and
phylogenetic mean-mean and q=4 random-intercept correlations.