---
title: "What can I fit today?"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{What can I fit today?}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

<style>
@media (max-width: 575.98px) {
  main table {
    display: block;
    max-width: 100%;
    overflow-x: auto;
  }
}
</style>

This guide is the user-facing status map for `drmTMB`. It is not a worked
tutorial. Use it when you want to know whether a model is implemented, which
article to read next, and which syntax belongs to the current package rather
than the roadmap.

If you have not fit a model yet, begin with
[Distributional regression with drmTMB](drmTMB.html), then return here when you
need the implemented-versus-planned boundary.

The main design rule is one formula per distributional parameter. For one
response, the most common parameters are the location `mu`, the residual scale
`sigma`, and family-specific shape or zero-inflation parameters. For two
responses, `mu1`, `mu2`, `sigma1`, `sigma2`, and `rho12` describe the two
means, two residual scales, and residual response-response correlation.
Equivalently, `mu` is the location term, `sigma` is the scale term,
family-specific parameters such as Student-t `nu` are shape terms, and
coscale means a residual-correlation parameter such as bivariate `rho12`.

Use this guide as the landing page when you are deciding which article to read
next. The quick rule is: use the worked tutorials for a fitted analysis, use
the model guides for vocabulary and implemented-versus-planned status, and use
the reference pages when you already know the function name.
For a finer table of random-effect types, dependence layers, q, random slopes,
`corpairs()`, `zi`, and `hu`, use the
[implementation map](implementation-map.html).

Start from the scientific question, then move to syntax and interpretation.
This keeps planned features from looking like fitted routes.

| If your question is... | Fit now | Interpret as | Planned neighbour |
| --- | --- | --- | --- |
| Do predictors change the expected response? | `mu ~ x` or `mu1 = y1 ~ x`, `mu2 = y2 ~ x` | location effects on the response mean or expected value | richer mixed-response bivariate families |
| Do predictors change residual variability? | `sigma ~ z`; ordinary NB2, lognormal, and Gamma can add `sigma ~ z + (1 | id)` | residual scale on the family-specific scale, often `log(sigma)` | non-Gaussian `sigma` slopes and random effects outside those three intercept gates |
| Do groups vary in average response or slopes? | `(1 | id)`, `(0 + x | id)`, or supported `(1 + x | id)` routes | ordinary group-level random effects in `mu` | correlated non-Gaussian slopes and broader p8/q8 endpoint variants |
| Do related species, sites, animals, or known matrices structure `mu`? | `phylo()`, coordinate `spatial()`, `animal()`, or `relmat()` in fitted Gaussian routes; ordinary Poisson and NB2 also support one q=1 structured `mu` intercept-plus-one-slope term through the same marker families | latent structured deviations, separate from residual `rho12` | pure, multiple, or labelled count structured slopes and richer non-Gaussian structured routes beyond the exact ordinary Poisson/NB2 q=1 gate |
| Do two responses remain correlated after their means and scales are modelled? | `rho12 ~ x` | residual coscale, not group or structured covariance | random effects in `rho12` |
| Do extra zeros, exact boundaries, or hurdle probability vary with predictors? | fixed-effect `zi ~ x` or `hu ~ x` where supported; a truncated-NB2 hurdle can use one diagnostic-only q=1 `hu ~ relmat(1 | id, K/Q = ...)` intercept; `zero_one_beta()` uses fixed-effect `zoi ~ x` and `coi ~ x`, plus exact point-fit-only ordinary `zoi ~ 1 + (1 | id)`, same-raw-symbol `zoi ~ x + (0 + x | id)`, `coi ~ 1 + (1 | id)`, or same-raw-symbol `coi ~ x + (0 + x | id)` q1 routes | probability-component fixed effects; the exact zero-one-beta atom gates support point recovery, not profiles, intervals, or coverage; sparse observed atoms or weak boundary-row predictor spread can weaken individual `coi` modes | random effects or structured dependence in `zi`, transformed/mismatched or broader atom slopes, joint atom effects, and `hu` routes beyond the exact q=1 `relmat()` intercept |

The map hierarchy is: this page helps users choose a fitted route, the
implementation map checks fitted-versus-planned boundaries, the source map
points contributors to code and tests, and the validation-debt register records
the evidence behind advertised surfaces.

The Reference index is also part of the workflow. Post-fit helpers are grouped
so that a reader can move from fitted models to summaries, predictions,
uncertainty, and plots without guessing which functions are user-facing:

| Reference group | Current user-facing functions | Use when... |
| --- | --- | --- |
| Model fitting and post-fit tools | `drmTMB()`, `check_drm()`, `summary()`, `confint()`, `profile_targets()`, `prediction_grid()`, `predict_parameters()`, `marginal_parameters()`, `fixef()`, `ranef()`, `sigma()`, `rho12()`, `corpairs()` | You need to fit, check, summarize, predict, or extract fitted quantities. |
| Visualization | `plot_parameter_surface()`, `plot_corpairs()` | You already have a `predict_parameters()` or `corpairs()` table and want an optional `ggplot2` display of fitted estimates and any explicit interval columns. |

Planned plotting helpers stay out of the Reference index until they are
exported, documented, tested, and connected to a stable data contract.

For Phase 17 visualization decisions, start from the table or extractor that
matches the question before choosing a plot:

| If you want to show... | Start from | Current display route | Check before styling |
| --- | --- | --- | --- |
| observed responses | the original analysis data | draw raw points or summaries on the observed-response scale | do not place raw response points on scale, correlation, or random-effect SD axes |
| fitted `mu`, `sigma`, shape, or zero-inflation surfaces | `prediction_grid()` then `predict_parameters(conf.int = TRUE)` when Wald fixed-effect intervals are needed | `plot_parameter_surface()` for estimate surfaces and explicit interval bands, or your own `ggplot2` layers | `dpar`, `type`, `conf.status`, `conf.level`, and `interval_source` columns |
| fitted random-effect SD surfaces | `prediction_grid()` then `predict_parameters(..., dpar = "sd(group)")`, optionally followed by `marginal_parameters()` | table first; custom plot only after residual `sigma` and random-effect SDs are separated | component `random-effect-sd-model`, group-level predictors, and interval status |
| adjusted or empirical marginal summaries | `prediction_grid(..., margin = "empirical")` then `marginal_parameters()` | table first; custom plot only after the averaging rule is clear | focal variables, `by`, margin metadata, and weights rule |
| residual, group-level, phylogenetic, or spatial correlations | `rho12()` or `corpairs()` | `plot_corpairs()` for explicit `corpairs()` tables; custom plots for planned layers | correlation layer, `conf.status`, `interval_source`, and whether the row is direct or derived |
| fitted intervals | `predict_parameters(conf.int = TRUE)`, `confint()`, or `summary(conf.int = TRUE)` | interval table first; ribbons are drawn only from explicit interval columns | `conf.status`, `conf.level`, `profile.boundary`, `profile.message`, and `interval_source` |
| diagnostics before interpretation | `check_drm()` | diagnostic table first; diagnostic plots remain planned | convergence, Hessian, boundary, design, and retained-object status |

| If the scientific phrase is... | It usually maps to... | Start here |
| --- | --- | --- |
| average response, expected trait value, or treatment effect | `mu` or `mu1` / `mu2` | [When variance carries signal, Part 1](location-scale.html) |
| residual variability, predictability, extra heterogeneity, or dispersion | `sigma` or `sigma1` / `sigma2` | [Which scale are you modelling?](which-scale.html) |
| among-group variation in expected response | random-effect SDs or `sd(group)` | [Which scale are you modelling?](which-scale.html) |
| residual trait coupling after the means and residual SDs are modelled | `rho12` | [Changing residual coupling with rho12](bivariate-coscale.html) |
| pedigree, tree, coordinate, or known-matrix dependence | `animal()`, `phylo()`, `spatial()`, or `relmat()` depending on the source | [Structural dependence overview](structural-dependence.html) |
| fitted interval status or boundary diagnostics | `profile_targets()`, `conf.status`, `profile.boundary`, `profile.message` | [Checking and using fitted models](https://itchyshin.github.io/drmTMB/articles/model-workflow.html) |

## Stable-core matrix

Use this matrix before treating a syntax pattern as available. "Stable" means a
routine fitted surface with tests and user-facing documentation. "First slice"
means the model fits, but only inside the stated boundary. "Opt-in control"
means the feature is a hardening or large-data path, not a promise that every
neighbouring model scales the same way. The evidence and debt ledger behind the
matrix lives in `docs/design/34-validation-debt-register.md`.

Read the status words this way:

| Status word | Meaning for a fitted analysis |
|---|---|
| Stable | Routine fitted path with tests, diagnostics or interval status, and a reader-facing example or guide. |
| Implemented, source-tested; ledger tier unregistered | The syntax and focused source tests exist, but no capability-ledger cell currently authorizes a stronger reader-facing maturity tier. |
| First slice | Fitted and tested, but intentionally narrow; stay inside the named formula, family, and data-shape boundary. |
| Opt-in control | Available for hardening, scalability, or memory control, but not a modelling guarantee for neighbouring surfaces. |
| Planned or reserved | Public grammar or roadmap wording may exist, but `drmTMB()` should reject it or treat it as design-only until likelihood, tests, docs, and after-task evidence land. |
| Unsupported or blocked | Do not use as analysis syntax; fit the nearest implemented model or follow the named fallback in [Can I fit and report this model?](capability-and-limits.html) before interpreting a richer structure. |

| Model surface | Status | Fitted today | Not yet | Interval and diagnostic status |
|---|---|---|---|---|
| One-response families | Stable, with first non-Gaussian random-effect slices | Gaussian, Student-t, lognormal, Gamma, Tweedie, beta, zero-one beta, beta-binomial, plain binomial logit, Poisson, NB2, zero-inflated Poisson, zero-inflated NB2, truncated NB2, hurdle NB2, and cumulative-logit ordinal location are fitted. Eligible ordinary routes across Student-t, skew-normal, lognormal, Gamma, Tweedie, beta, zero-one beta, beta-binomial, binomial, Poisson, NB2, truncated NB2, and cumulative-logit fit `mu` random intercepts plus independent slopes. Additional exact gates include ordinary NB2/lognormal/Gamma log-`sigma` random intercepts, Poisson/NB2 q=1 `phylo()`/`spatial()`/`animal()`/`relmat()` `mu` intercept-plus-one-slope routes, recovery-grade Gamma-phylo and lognormal-phylo/relmat q1 `mu` intercepts, separate recovery-grade NB2 q=1 structured `sigma` routes, one crossed NB2 `mu ~ spatial() + relmat()` recovery-only route, and one truncated-NB2 q=1 `hu ~ relmat(1 | id, K/Q = ...)` diagnostic-only route | correlated or labelled non-Gaussian slopes, binomial structured effects or Julia bridge, ordinal scale/discrimination, ordered beta, Tweedie distributional random effects or predictor-dependent power, non-Gaussian `sigma` random effects outside the named gates, shape random effects beyond the exact Student-t phylogenetic gate, zero-one beta distributional random effects beyond ordinary `mu`, zero-inflated count random effects beyond the exact gates, hurdle random effects beyond the exact q=1 `hu ~ relmat()` intercept, pure, multiple, or labelled structured count slopes, simultaneous structured count types beyond the exact crossed NB2 gate, and unsupported neighbouring family/parameter routes | Wald fixed-effect intervals by default when `sdreport()` is computed; binomial and Tweedie fixed-effect rows use the fixed-effect interval path; direct profile targets appear in `profile_targets()`, including binomial, Poisson, NB2, Student-t, zero-truncated NB2, lognormal, Gamma, beta, and beta-binomial `mu` random-effect SDs through `log_sd_mu`, NB2/lognormal/Gamma `sigma` random-intercept SDs through `log_sd_sigma`, and Poisson/NB2 structured SDs through `log_sd_phylo`; the exact `mc-0061`, `mc-0464`, `mc-0539`, and `mc-0575` slopes plus `mc-0382` lognormal sigma-intercept domain are inference-ready with caveats under their separate designs; public `mc-0227` remains ML-Laplace point-fit-recovery only, and its O3 campaign is internal and non-reportable |
| Gaussian ordinary random effects | Stable, with advanced q > 2 blocks | `mu` intercepts, independent slopes, one-slope correlated blocks, ordinary q > 2 numeric multi-slope `mu` blocks, `sigma` intercepts, independent slopes, and unlabelled correlated intercept-slope or multi-slope blocks on log-`sigma` | labelled residual-scale or cross-formula `mu`-`sigma` slope covariance, coefficient-specific `sd()` slope models, and treating large q blocks as routine | `check_drm()` reports replication, weak-slope, boundary, and Hessian diagnostics; q=3 recovery and q=4 output-contract checks cover ordinary multi-slope paths; q > 2 block SDs are direct profile targets, while q > 2 correlations are derived-unavailable for direct profiling |
| Random-effect scale models | First slice | `sd(group) ~ x_group` for unlabelled Gaussian `mu` random intercepts | slope-specific `sd(id, dpar = "mu", coef = "x") ~ ...` and residual-scale direct-SD models | SD-surface coefficients are direct; row-specific group SD summaries are derived |
| Known sampling covariance | Implemented, source-tested; ledger tier unregistered | Gaussian `meta_V(V = V)` for diagonal, dense, and row-paired bivariate known sampling covariance, with deprecated `meta_known_V(V = V)` as a compatibility alias | non-Gaussian known covariance, sparse known covariance, broad dense-`V` scalability claims, and dense known `V` with non-unit likelihood weights | `check_drm()` reports dense full `V` as a note with dimension, density, size, rank, and conditioning; fixed effects and response-scale residual summaries use the usual interval routes only when Hessian diagnostics are clean. Some `meta_V(V = V)` fits with predictor-dependent `sigma` can return plausible point estimates while reporting `pdHess = FALSE`; treat their Wald SEs and intervals as unreliable until a profile, bootstrap, or simpler `sigma` model supports the target |
| Bivariate Gaussian residual coupling | Stable | fixed-effect `mu1`, `mu2`, `sigma1`, `sigma2`, and predictor-dependent residual `rho12` | random effects in `rho12` and mixed composed families | `rho12()` reports response-scale residual correlation; `confint(..., parm = "rho12", newdata = ...)` profiles supplied rows, but current ledger evidence does not certify their CI coverage |
| Ordinary bivariate covariance and `corpairs()` | First slice | matching labelled random intercepts in `mu1`/`mu2`, `sigma1`/`sigma2`, same-response `mu`/`sigma`, all-four q=4 intercept blocks, matching slope-only `mu1`/`mu2`, same-response `mu`/`sigma`, and `sigma1`/`sigma2` blocks, matching q4/q6 `mu1`/`mu2` location blocks with smoke artifact routing, the first q8 all-endpoint block with smoke/recovery artifact routing, and q=2 `corpair(..., level = "group") ~ x` | broader p8/q8 endpoint variants and predictor-dependent slope `corpair()` regressions | constant q=2 SD/correlation targets, the slope-slope row, the same-response mean-scale-slope row, and the scale-slope row are direct; q > 2 location SDs and q8 endpoint SDs are direct targets; predictor-dependent `corpair()` values use `newdata`; q > 2 correlations, including q8 endpoint correlations, are derived and mark derived intervals unavailable; q8 has diagnostic artifacts, not coverage or power evidence |
| Phylogenetic structured effects | First slices | Gaussian univariate `mu` and `sigma` intercepts with optional matching `mu`/`sigma` correlation, one numeric `mu` slope, the exact q1 `sigma` one-slope route, Gaussian bivariate `mu1`/`mu2`, labelled q=4 location-scale blocks, `sd_phylo*()` direct-SD surfaces, q=2 phylogenetic `corpair()` regression, ordinary Poisson/NB2 q=1 `mu` intercept-plus-one-slope routes, and a separate recovery-grade NB2 q1 structured `sigma` intercept-plus-one-slope route | pure, multiple, or labelled Poisson/NB2 phylogenetic slopes, zero-inflated phylogenetic effects, multiple or labelled Gaussian phylogenetic slopes, slope correlations, direct-SD formulas combined with structured `sigma`, structured `rho12`, and predictor-dependent q=4 phylogenetic correlations | the exact Gaussian q1 sigma one-slope route is inference-ready with caveats under raw log-SD Wald-z intervals; direct phylogenetic SD and constant q=2 correlation targets can be mechanically profile-ready, but profile availability alone is diagnostic and does not validate an interval; predictor-dependent `corpair()` values use `newdata` for diagnostic profiles; full q=4 correlations are derived-only, while block-diagonal q=4 fallback correlations are direct targets but still need fit-specific profile diagnostics; the Poisson/NB2 q=1 routes are recovery-grade or smoke-level with direct `log_sd_phylo` targets, not interval/coverage promotions |
| Coordinate spatial structured effects | First slices | `spatial(1 | site, coords = coords)` in univariate Gaussian `mu` and/or `sigma`, matching univariate `mu`/`sigma` correlation, one numeric `mu` slope, a q1 `sigma` one-slope point-fit/extractor route, matching bivariate `mu1`/`mu2` terms, matching all-four q=4 spatial terms, ordinary Poisson/NB2 q=1 `mu` intercept-plus-one-slope routes, a separate recovery-grade NB2 q1 structured `sigma` intercept-plus-one-slope route, the exact diagnostic-only Poisson `zi ~ spatial()` intercept, the exact diagnostic-only fixed-`zi` Poisson `mu ~ spatial()` intercept, and the exact diagnostic-only fixed-`zi` NB2 `mu ~ spatial()` intercept | mesh/SPDE beyond the fixed-kappa Gaussian `mu` intercept, multiple or labelled slopes, spatial sigma-slope intervals and coverage, slope correlations, spatial direct-SD, spatial `corpair()`, pure, multiple, or labelled count spatial slopes, labelled count covariance, and zero-inflated spatial effects outside the exact Poisson `zi`, fixed-`zi` Poisson `mu`, and fixed-`zi` NB2 `mu` gates | `summary()$parameters`, marker-specific `ranef()$terms`, `profile_targets()`, `check_drm()`, `corpairs(level = "spatial")`, and `summary()$covariance` expose the fitted coordinate fields; both fixed-`zi` spatial-`mu` gates have no recovery, interval, or coverage promotion; the spatial sigma-slope interval gate remains blocked; q=4 spatial correlations are derived-only and currently have extractor/diagnostic smoke rather than formal coverage evidence |
| Animal and lower-level relatedness markers | First slice | univariate Gaussian `mu` and `sigma` random intercepts with `animal(1 | id, pedigree/A/Ainv = ...)` or `relmat(1 | id, K/Q = ...)`, optional matching `mu`/`sigma` correlation, one numeric `mu` slope, the exact A-matrix animal and K/Q relmat q1 `sigma` one-slope routes, matching labelled bivariate q=2 `mu1`/`mu2` location covariance, constant all-four q=4 location-scale blocks, ordinary Poisson/NB2 q=1 `mu` intercept-plus-one-slope routes, and separate recovery-grade NB2 q1 structured `sigma` intercept-plus-one-slope routes | pedigree/Ainv bridge marshalling, sparse large-pedigree precision construction, multiple or labelled structured slopes, slope correlations, predictor-dependent `corpair()` regression, pure, multiple, or labelled count animal/`relmat()` slopes or labelled covariance, and generic direct-SD grammar | the exact animal/relmat sigma slopes are inference-ready with caveats under raw log-SD Wald-z intervals; direct SD and constant q=2 correlation targets can be mechanically profile-ready, but profile availability alone is diagnostic and does not validate an interval; q=4 correlation rows are derived-only; `corpairs()`, `summary()$covariance`, `profile_targets()`, and `check_drm()` report the fitted structured rows; dense `pedigree`, `K`, and `A` examples are small-to-moderate, while scalable claims should prefer sparse `Q`/`Ainv` evidence |
| Profile intervals and diagnostics | First slice | fast Wald intervals for fixed effects and direct scale, SD, and correlation targets; selected direct profile targets; row-specific `sigma`, `sigma1`, `sigma2`, `rho12`, and fitted q=2 `corpair()` values; and `confint(..., method = "bootstrap")` simulate/refit intervals for direct targets | automatic intervals for every derived covariance summary and bootstrap support outside direct `confint()` targets | `conf.status`, `profile.boundary`, and `profile.message` expose interval status; availability is not a coverage certification, especially for bivariate fixed-effect `rho12` rows. SD Wald intervals use the fitted log-SD scale, correlation Wald intervals use a guarded Fisher-z/atanh scale, `profile_precision = "fast"` gives a quicker first-pass profile, q=4 rows use `derived_interval_unavailable`, and bootstrap rows report successful and failed refits |
| Large-data fit controls | Opt-in control | memory-light fitted objects, optional `se = FALSE` standard-error skipping, sparse fixed-effect `mu` matrices, and Gaussian sufficient-statistic aggregation | broad random-effect, structured-effect, non-Gaussian, bivariate, and known-covariance scalability claims | `check_drm()` reports `sdreport`, sparse-design, and aggregation diagnostics where fitted |
| Reserved or planned neighbours | Reserved/rejected or design-only | coefficient-specific `sd()` slopes, random effects in `rho12`, shape random effects, ID-level skewness such as future `skew(id) ~ x`, multiple or labelled phylogenetic slopes, non-Gaussian phylogenetic slopes outside the exact unlabelled Poisson/NB2 q1 intercept-plus-one-slope gates, phylogenetic slope correlations, mesh/SPDE beyond the fixed-kappa Gaussian intercept, spatial `corpair()`, broader or mismatched residual-scale/location-scale endpoint bivariate slope covariance beyond the fitted exact q2 routes, and mixed composed families | runnable analysis syntax | planned-feature errors should appear before fitting; no interval target is advertised |

The mesh exception is deliberately narrow: `bf(y ~ spatial(1 | site, mesh = mesh),
sigma ~ 1)` is an ML, fixed-kappa, univariate Gaussian `mu` intercept on
projected coordinates. It reports a raw GMRF field scale at local-fit level;
all other mesh/SPDE entries in this model map remain deferred.

## One response

The most mature path is Gaussian location-scale regression. Use it when the
question is about the expected response and residual predictability:

```r
drmTMB(
  drm_formula(y ~ x1, sigma ~ x2),
  family = gaussian(),
  data = dat
)
```

Here `y ~ x1` models the conditional mean `mu`, while `sigma ~ x2` models the
residual standard deviation. In Gaussian models, `sigma` coefficients are on a
log-SD scale, so exponentiating a coefficient gives a residual SD ratio.

Implemented one-response surfaces include:

| Question | Implemented surface | Read next |
|---|---|---|
| Do predictors change the mean and residual SD? | `drm_formula(y ~ x, sigma ~ z)` | [When variance carries signal, Part 1](location-scale.html) |
| Do groups differ in expected response? | ordinary `mu` random intercepts and one-slope blocks | [When variance carries signal, Part 1](location-scale.html) |
| Do groups differ in residual SD? | `sigma ~ x + (1 | group)` | [Which scale are you modelling?](which-scale.html) |
| Do groups differ in residual-scale slopes? | `sigma ~ x + (0 + w | group)` for independent Gaussian residual-scale slopes | [Which scale are you modelling?](which-scale.html) |
| Does a group-level predictor change among-group SD in the mean model? | `sd(group) ~ x_group` for unlabelled Gaussian `mu` random intercepts | [Which scale are you modelling?](which-scale.html) |
| Is the response continuous, count, bounded, robust, or zero-inflated? | fixed-effect non-Gaussian families | [Choosing response families](distribution-families.html) |
| Do repeated groups need non-Gaussian mean random effects? | selected ordinary non-Gaussian `mu` random intercepts and independent numeric slopes | [Implementation map](implementation-map.html) |
| Are the large residuals changing the location and scale conclusions? | fixed-effect Student-t `mu`, `sigma`, and `nu` formulas | [Robust continuous responses](robust-student.html) |
| Is the response an event probability, successes out of known trials, a strict continuous proportion, or a structural-boundary continuous proportion? | `stats::binomial(link = "logit")` for ordinary binomial event probabilities, `beta_binomial()` for overdispersed success counts, `beta()` for strict continuous proportions, and `zero_one_beta()` for exact structural 0/1 values; each has ordinary `mu` random intercepts and independent numeric slopes; exact binomial and zero-one-beta slope cells have design-specific inference-ready-with-caveats evidence | [Proportions and success rates](proportion-beta-binomial.html) |
| Does the response have known sampling variance or covariance? | Gaussian `meta_V(V = V)` | [Mean effects and residual heterogeneity](meta-analysis.html) |
| Do species means remain similar after shared ancestry is included? | `phylo(1 | species, tree = tree)` in univariate Gaussian `mu` | [Structural dependence overview](structural-dependence.html) |
| Are count means phylogenetically structured? | `phylo(1 | species, tree = tree)` in ordinary Poisson or ordinary NB2 `mu` | [Implementation map](implementation-map.html) |
| Do nearby sites share smooth location deviations or smooth slope differences? | `spatial(1 | site, coords = coords)` or one numeric `spatial(1 + x | site, coords = coords)` slope in univariate Gaussian `mu` | [Structural dependence overview](structural-dependence.html) |

The random-effect scale terms are deliberately separate. `sigma ~ x` models
within-observation residual SD. `sd(group) ~ x_group` models the standard
deviation of a group-level mean effect. Those two scales can answer different
biological questions even when they use similar predictors.

For ordinary grouped Gaussian models, read the random-effect layer in this
order:

| Model phrase | Syntax | Read with |
|---|---|---|
| groups differ in baseline `mu` | `(1 | id)` | `summary(fit)$parameters`, `ranef(fit, "mu")$terms` |
| groups differ in the `x` slope | `(0 + x | id)` | `summary(fit)$parameters`, `profile_targets(fit)` |
| group baselines and slopes are correlated | `(1 + x | id)` or `(1 + x | p | id)` | `summary(fit)$covariance`, `corpairs(fit, class = "mean-slope")` |
| groups differ in residual scale | `sigma ~ z + (1 | id)` | `summary(fit)$parameters`, `sigma(fit)` |
| groups differ in residual-scale slopes | `sigma ~ z + (0 + w | id)` | `summary(fit)$parameters`, `profile_targets(fit)`, `sigma(fit)` |
| group-level predictors change among-group `mu` SD | `sd(id) ~ x_group` | `coef(fit, "sd(id)")`, `prediction_grid()`, `predict_parameters(..., dpar = "sd(id)")`, `marginal_parameters()` |

The last column is the output contract. It keeps among-group `mu` variation,
residual `sigma`, group-level random-effect correlations, and residual `rho12`
as different model layers.

## Two responses

The implemented bivariate Gaussian path estimates response-specific means,
response-specific residual SDs, and residual coupling `rho12`:

```r
drmTMB(
  drm_formula(
    mu1 = y1 ~ x1 + x2,
    mu2 = y2 ~ x1,
    sigma1 = ~ x1 + x2,
    sigma2 = ~ x1,
    rho12 = ~ x1 + x2
  ),
  family = c(gaussian(), gaussian()),
  data = dat
)
```

`rho12` is the residual correlation between the two responses after the model
has already accounted for the predictors in `mu1`, `mu2`, `sigma1`, and
`sigma2`. Extract fitted response-scale residual correlations with
`rho12(fit)`, and use `corpairs(fit)` when you want a long table that keeps
residual correlations distinct from group-level correlations.

The first bivariate group-level covariance slice is also implemented for
matching labelled random intercepts in the two location formulas:

```r
drmTMB(
  drm_formula(
    mu1 = y1 ~ x1 + (1 | p | id),
    mu2 = y2 ~ x1 + (1 | p | id),
    sigma1 = ~ x1,
    sigma2 = ~ x1,
    rho12 = ~ x1
  ),
  family = c(gaussian(), gaussian()),
  data = dat
)
```

The `p` label creates a group-level covariance block for the `mu1` and `mu2`
random intercepts. That group-level correlation is not `rho12`: it asks
whether groups with higher average response 1 also tend to have higher average
response 2. `rho12` still describes within-observation residual coupling.

## Structural dependence

Read structural-dependence syntax in biological order before reading the
implementation order. Animal models use a known pedigree or additive
relatedness matrix. Phylogenetic models use a tree or tree-derived covariance.
Spatial models use coordinates or the bounded fixed-kappa mesh/SPDE Gaussian
intercept. Some
analyses need phylogenetic and spatial layers in the same predictor, and some
advanced users need a generic known relatedness or precision matrix through
`relmat()`.

The current fitted paths cover first Gaussian structured pieces for
`phylo()`, coordinate `spatial()`, `animal()`, and `relmat()`. Each route fits
a one-response Gaussian `mu` random intercept and one numeric `mu` slope, for
example `phylo(1 + x | species, tree = tree)`,
`spatial(1 + x | site, coords = coords)`,
`animal(1 + x | id, Ainv = Ainv)`, or `relmat(1 + x | id, Q = Q)`. The same
structured families also include selected intercept-only `sigma` routes and q1
`sigma` one-slope routes, matching labelled q=2 bivariate `mu1`/`mu2`
covariance rows, and constant all-four q=4 location-scale blocks where marked.
Phylo, A-matrix animal, and K/Q relmat sigma slopes are inference-ready with
caveats; spatial remains point-fit/extractor only. Large-pedigree sparse
construction, pedigree/Ainv bridge marshalling, multiple or labelled slopes,
predictor-dependent `corpair()` regression, and generic direct-SD grammar
remain planned.

```r
drmTMB(
  drm_formula(
    y ~ x1 + phylo(1 | species, tree = tree),
    sigma ~ x1
  ),
  family = gaussian(),
  data = dat
)
```

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

Each term adds a structured random effect to the univariate Gaussian `mu`
predictor, the `sigma` predictor, or a matching univariate `mu`/`sigma` pair.
Matching intercept-only `phylo()` terms in bivariate `mu1` and `mu2` formulas
fit the first phylogenetic mean-mean correlation slice. Matching labelled
`phylo()` terms across `mu1`, `mu2`, `sigma1`, and `sigma2` also fit the first
constant q=4 phylogenetic location-scale covariance block. The coordinate-based
`spatial()` route fits `spatial(1 | site, coords = coords)` in `mu` and/or
`sigma`, `spatial(1 + x | site, coords = coords)` in `mu`, the first bivariate
q=2 `mu1`/`mu2` spatial covariance, and the constant q=4 location-scale block
when matching labelled spatial terms appear in all four endpoints. The q1
spatial `sigma` one-slope route also has point-fit/extractor evidence, but its
interval gate remains blocked. Mesh/SPDE spatial fields, spatial direct-SD
models, multiple or labelled spatial slopes, slope correlations,
non-Gaussian spatial effects outside the exact ordinary Poisson/NB2 q1 spatial
`mu` intercept-plus-one-slope, recovery-grade NB2 q1 spatial `sigma`, Student-t
spatial `mu`, Poisson spatial `zi`, fixed-`zi` Poisson spatial `mu`, and
fixed-`zi` NB2 spatial `mu` gates, and
spatial `corpair()` regressions remain planned extensions.
Run `check_drm()` before interpreting that fitted phylogenetic correlation; it
flags near-boundary phylogenetic correlation rows and tiny phylogenetic
location SDs relative to the matching residual scales.

For animal-model syntax, the first runnable route uses a precomputed
relationship matrix:

```r
drmTMB(
  drm_formula(
    body_size ~ sex + animal(1 | individual, Ainv = Ainv),
    sigma ~ cohort
  ),
  family = gaussian(),
  data = dat
)
```

The ordinary Gaussian random-effect fallback remains useful when the pedigree
or relatedness matrix is unavailable:

```r
drmTMB(
  drm_formula(
    body_size ~ sex + (1 | individual),
    sigma ~ cohort
  ),
  family = gaussian(),
  data = dat
)
```

Report that fallback as repeatability or grouped heterogeneity, not additive
genetic variance. It ignores pedigree or `Ainv` structure.

Use `relmat(1 | id, K = K)` or `relmat(1 | id, Q = Q)` for lower-level
examples where the known dependence matrix is neither naturally a pedigree, a
phylogeny, nor a coordinate-spatial surface.
If that matrix is known sampling covariance among observations or effect-size
estimates, use `meta_V(V = V)` instead; `relmat()` is for latent random-effect
relatedness or precision matrices.

For fitted full q=4 phylogenetic blocks, read the six latent covariance rows
with `corpairs(fit, level = "phylogenetic")` or
`summary(fit)$covariance`. Those rows are not residual `rho12`, and they remain
derived and unavailable for intervals. Block-diagonal q=4 fallback fits expose two constant
block correlations as direct targets, but a direct target can still fail on a
weak, one-sided, or boundary-limited profile. Predictor-dependent q=4
`corpair()` regressions remain planned.

## Planned, not implemented

These surfaces are roadmap items, not syntax to rely on in analyses yet:

| Planned surface | Why it is separate |
|---|---|
| labelled residual-scale slope blocks and larger cross-formula labelled covariance blocks linking `mu` and `sigma` | unlabelled same-`sigma` correlated slopes are fitted, but labelled or cross-formula slope covariance extends beyond the validated block |
| residual-scale bivariate random slopes and full cross-parameter covariance blocks | expands group-level covariance blocks beyond the first same-parameter and same-response random-intercept slices and the matching slope-only plus q4/q6 location `mu1`/`mu2` slices |
| multiple phylogenetic slopes, slope correlations, predictor-dependent q=4 phylogenetic correlations, and richer phylogenetic covariance blocks | the first univariate one-slope `mu` path, bivariate `mu1`/`mu2` phylogenetic location slice, q=2 predictor-dependent location-location slice, and first constant q=4 phylogenetic location-scale block are fitted |
| richer spatial random effects | first coordinate-spatial Gaussian `mu` and `sigma` intercepts, one-slope `mu`, a q1 `sigma` one-slope point-fit/extractor route, bivariate q=2 location-covariance, constant q=4 location-scale paths, ordinary Poisson/NB2 q=1 spatial `mu` intercept-plus-one-slope routes, and a separate recovery-grade NB2 q1 spatial `sigma` intercept-plus-one-slope route are implemented; the Gaussian spatial sigma-slope interval gate remains blocked, while mesh/SPDE beyond the fixed-kappa Gaussian `mu` intercept, additional multiple or labelled Gaussian-slope layouts outside the exact fitted ledger cells, slope correlations, pure, multiple, or labelled count spatial slopes, and spatial `corpair()` paths remain planned |
| animal models and generic relatedness matrices | Gaussian `mu` and `sigma` intercepts plus one numeric `mu` slope are fitted for `animal(pedigree/A/Ainv)` and `relmat(K/Q)`; the exact A-matrix animal and K/Q relmat q1 `sigma` one-slope routes are fitted and inference-ready with caveats; matching labelled `mu1`/`mu2` terms fit the first q=2 bivariate location covariance, matching all-four `mu1`/`mu2`/`sigma1`/`sigma2` terms fit constant q=4 location-scale blocks, and ordinary Poisson/NB2 q=1 `animal()` or `relmat()` `mu` intercept-plus-one-slope routes plus separate recovery-grade NB2 q1 structured `sigma` routes are fitted; pedigree/Ainv bridge marshalling, sparse large-pedigree construction, additional multiple or labelled Gaussian-slope layouts outside the exact fitted ledger cells, slope correlations, predictor-dependent `corpair()`, pure, multiple, or labelled count slopes, and direct-SD grammar remain planned |
| skew-normal random, structured, bivariate, known-covariance, `rho12`, and latent-skew extensions | fixed-effect `skew_normal()` with `nu ~ ...` is fitted and has focused source plus Phase 18 artifact tests; these richer routes need separate likelihood, extractor, diagnostic, interval, comparator, and recovery evidence |
| higher-dimensional multivariate models | belong in `gllvmTMB`, not `drmTMB` |

When unsupported syntax fails, the right next step is usually to fit the
nearest implemented fixed-effect or univariate model, then use
[Can I fit and report this model?](capability-and-limits.html) for the exact
reporting boundary and fallback before interpreting a richer covariance
structure as available.

## Correlation layers

Several model components can be called correlations, but they answer different
questions:

| Layer | Example extractor or surface | Interpretation |
|---|---|---|
| residual | `rho12(fit)` or residual rows from `corpairs(fit)` | within-observation coupling between two responses |
| ordinary group-level | group-level rows from `corpairs(fit)` or `summary(fit)$covariance` | correlation among random intercepts or slopes |
| phylogenetic | `corpairs(fit, level = "phylogenetic")` | correlation among tree-structured species mean, scale, or q=2 modelled location deviations, depending on the fitted block |
| spatial | `ranef(fit, "spatial_mu")$terms`, `ranef(fit, "spatial_sigma")$terms`, `summary(fit)$parameters`, `check_drm()`, and `corpairs(fit, level = "spatial")` for fitted univariate, bivariate q=2, and q=4 spatial rows | coordinate-structured site intercept, one numeric `mu` slope, or the q1 `sigma` one-slope point-fit/extractor route in one response; univariate location-scale, bivariate coordinate-spatial mean-mean, and constant location-scale correlations are fitted, while spatial sigma-slope intervals, predictor-dependent correlation, multiple or labelled slopes, slope-correlation, and non-Gaussian paths outside the exact ordinary Poisson/NB2 q1 spatial `mu` intercept-plus-one-slope gates, recovery-grade NB2 q1 spatial `sigma`, Student-t spatial `mu`, Poisson spatial `zi`, fixed-`zi` Poisson spatial `mu`, and fixed-`zi` NB2 spatial `mu` gates remain planned |

Keep these layers separate in writing. A residual `rho12` result is not a
phylogenetic correlation, a spatial correlation, or a personality/plasticity
correlation unless the fitted model actually contains that higher-level
covariance block.

## A practical trait protocol

Comparative trait studies often want one model that combines two trait means,
two residual SDs, residual trait coupling, individual or species differences,
and shared ancestry. That is the right scientific direction, but the fully
combined bivariate phylogenetic double-hierarchical model is still a roadmap
target. Today, fit the pieces that are implemented and report them as separate
answers. The first bivariate phylogenetic location block can now join that
staged protocol when both trait means use matching intercept-only `phylo()`
terms.

For two traits without phylogeny, use the bivariate Gaussian path:

```r
fit_biv <- drmTMB(
  drm_formula(
    mu1 = log_body_mass ~ diet + body_size + (1 | p | species),
    mu2 = log_litter_size ~ diet + body_size + (1 | p | species),
    sigma1 = ~ diet,
    sigma2 = ~ diet,
    rho12 = ~ diet
  ),
  family = c(gaussian(), gaussian()),
  data = mammals
)
```

This model can ask whether `diet` changes each trait mean, each residual SD, and
the residual association between the two traits. The shared `p` label adds the
implemented group-level `mu1`/`mu2` random-intercept covariance for `species`.
It is an ordinary species-level grouping term, not a phylogenetic covariance.

For shared ancestry in one trait, use the implemented univariate phylogenetic
path:

```r
fit_phylo <- drmTMB(
  drm_formula(
    log_body_mass ~ diet + body_size + phylo(1 | species, tree = tree),
    sigma ~ diet
  ),
  family = gaussian(),
  data = mammals
)
```

This model can ask whether species trait means remain similar after shared
ancestry is included. It does not estimate bivariate residual `rho12`, and it
does not estimate phylogenetic correlations among `mu1`, `mu2`, `sigma1`, and
`sigma2`.

For shared ancestry in two trait means, use matching bivariate `phylo()` terms:

```r
fit_biv_phylo <- drmTMB(
  drm_formula(
    mu1 = log_body_mass ~ diet + body_size + phylo(1 | species, tree = tree),
    mu2 = log_litter_size ~ diet + body_size + phylo(1 | species, tree = tree),
    sigma1 = ~ diet,
    sigma2 = ~ diet,
    rho12 = ~ diet
  ),
  family = c(gaussian(), gaussian()),
  data = mammals
)
```

This model can ask whether phylogenetic deviations in one trait mean are
associated with phylogenetic deviations in the other trait mean. The residual
`rho12` term still answers a separate within-observation question after those
mean deviations are included. Read the fitted phylogenetic layer with
`corpairs(fit_biv_phylo, level = "phylogenetic")` or
`summary(fit_biv_phylo)$covariance`. A `confint()` call on the explicit
`cor:phylo:` target can diagnose profile geometry, but this intercept-only q2
row is not interval-validated and the diagnostic profile is not a reporting
interval.

If the scientific question separates shared ancestry from a remaining
species-level association, that is an advanced identifiability comparison, not
the introductory path. A phylogenetic species effect is species-level structure
induced by the tree; it is not a residual observation-level correlation.
Adding a non-phylogenetic species block beside `phylo()` asks the model to split
two species-level covariance layers that share the same grouping factor. That
comparison should live in a dedicated identifiability example with
`check_drm()`, profile intervals, and simpler-model comparisons, rather than in
the main model map.

Use this staged interpretation until the full combined location-scale model is
implemented:

| Fitted piece | What it answers today | What it does not answer yet |
|---|---|---|
| `rho12(fit_biv)` | within-observation residual coupling between two traits | phylogenetic or species-level correlation |
| `corpairs(fit_biv, level = "group")` | ordinary species or individual random-intercept correlation | residual trait coupling or shared-ancestry correlation |
| `fit_phylo` with `phylo(1 | species, tree = tree)` | one-trait shared-ancestry mean structure | bivariate phylogenetic covariance |
| `corpairs(fit_biv_phylo, level = "phylogenetic")` | bivariate phylogenetic mean-mean correlation, or six q=4 rows when the all-four labelled block is fitted | predictor-dependent phylogenetic correlations |
| `summary(fit_biv_phylo)$covariance` | variance and covariance point summaries for the fitted phylogenetic mean-mean or q=4 layer | derived-profile intervals for q=4 correlations |

For intervals, read the status column before reading the bounds.
`conf.status = "profile"` means a profile interval was returned;
`newdata_required` means the fitted surface needs a row supplied to
`confint(..., newdata = ...)`; and `derived_interval_unavailable` means the
point estimate is not a one-parameter profile target yet. Profile rows also
carry `profile.boundary` and `profile.message` so boundary-like intervals are
not mistaken for ordinary symmetric intervals.

In developer shorthand, the q=4 phylogenetic endpoint means four distributional
endpoints: `mu1`, `mu2`, `sigma1`, and `sigma2`. It does not mean four fitted
correlations. A four-endpoint covariance has six pairwise correlations; q=4
fits should therefore be read through the six `corpairs()` rows and treated as
a larger covariance model than the mean-mean path.

## Where to go next

| If you need... | Start with... |
|---|---|
| a first fitted model | [Getting started](drmTMB.html) |
| a worked Gaussian mean-variance analysis | [When variance carries signal, Part 1](location-scale.html) |
| a scale vocabulary check | [Which scale are you modelling?](which-scale.html) |
| a robust continuous-response example | [Robust continuous responses](robust-student.html) |
| a fixed-effect count example with overdispersion or structural zeros | [Count abundance and extra zeros](count-nbinom2.html) |
| a two-response residual-correlation example | [Changing residual coupling with rho12](bivariate-coscale.html) |
| known sampling variances or covariance | [Mean effects and residual heterogeneity](meta-analysis.html) |
| implemented phylogenetic random intercepts or coordinate-spatial intercept/slope terms | [Structural dependence overview](structural-dependence.html) |
| post-fit checks, predictions, residuals, and simulation | [Checking and using fitted models](https://itchyshin.github.io/drmTMB/articles/model-workflow.html) |