Use this guide to choose a route from your scientific question. It is not a tutorial: each linked guide gives the data setup, runnable example, and interpretation. If you have not fit a model yet, begin with the route that matches your data below.
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.
Choose a route
| If your data and question are… | Eligible route | Key limitation | Next guide |
|---|---|---|---|
| One response, and predictors may change its average or residual variability | Gaussian location-scale regression, then a family matched to the response | A family changes the meaning of mu, sigma,
and any extra parameters. |
Distributional regression, then Choosing response families |
| A continuous trait or response measured across related species, with a tree | Gaussian phylo() model |
You need enough replication to distinguish tree-structured and residual variation. | Phylogenetic mixed models |
| Effect sizes with known sampling variances or covariance | Gaussian meta_V(V = V) model |
It is not a response family for raw observations; check diagnostics before using uncertainty. | Mean effects and residual heterogeneity |
| Two responses measured together | Bivariate Gaussian model with response-specific means and scales |
rho12 is residual coupling, not a correlation among
groups or species. |
Changing residual coupling with rho12 |
| Counts, proportions, excess zeros, or robust continuous data | An eligible fixed-effect response-family route | Do not add a random or structured term merely because the neighbouring family accepts one. | Choosing response families |
| A tree, coordinates, pedigree, or known relationship matrix structures observations | A Gaussian structured-effect route | Structured deviations are distinct from residual variation and residual correlation. | Structural dependence overview |
Before reporting an estimate or interval, run
check_drm() and read Can I fit and report this model?.
If your model combines a family, random effect, or structured term not
shown in its worked guide, check that page before extending the formula:
nearby examples are not automatically interchangeable.
After fitting, start with check_drm(), then use
summary(), sigma(), rho12(), or
corpairs() according to the parameter you set out to
estimate. Use an interval only when the fitted object reports a usable
status for that specific quantity.
| 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 |
| residual variability, predictability, extra heterogeneity, or dispersion |
sigma or sigma1 / sigma2
|
Which scale are you modelling? |
| among-group variation in expected response | random-effect SDs or sd(group)
|
Which scale are you modelling? |
| residual trait coupling after the means and residual SDs are modelled | rho12 |
Changing residual coupling with rho12 |
| pedigree, tree, coordinate, or known-matrix dependence |
animal(), phylo(), spatial(),
or relmat() depending on the source |
Structural dependence overview |
| fitted interval status or boundary diagnostics |
profile_targets(), conf.status,
profile.boundary, profile.message
|
Checking and using fitted models |
Checking a less common model
Use this detailed table only when the first routes above do not cover your analysis. Read it from left to right: choose the model feature, check the documented form that is available now, then read the reporting caveat before using an interval. A model that fits is not automatically a model whose intervals are ready to report. This is not a menu of features to combine.
| Model feature | Scope | Documented form available now | Not yet supported | What to check before reporting |
|---|---|---|---|---|
| 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 tested forms 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 forms,
shape random effects beyond the exact Student-t phylogenetic form,
zero-one beta distributional random effects beyond ordinary
mu, zero-inflated count random effects beyond the exact
forms, 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 model, 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; selected independent mu slopes
for binomial, skew-normal, Tweedie, and zero-one beta, plus the
lognormal sigma random intercept, have assessed intervals
only within their separately tested designs and limitations; see the reporting guide. The public
cumulative-logit independent mu slope supports
point-estimate recovery under ML-Laplace only; an interval study using
an estimator unavailable through drmTMB() does not justify
reporting its public intervals |
| 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 and tested; check uncertainty for the specific model | 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 the available evidence does not establish their
confidence-interval 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 basic fitting
checks, the first q8 all-endpoint block with basic fitting and
parameter-recovery checks, 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 checks, 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 have parameter-recovery studies or basic fitting checks, with
direct log_sd_phylo targets, without evidence supporting
interval reporting or coverage |
| 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 forms |
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 forms lack evidence
supporting parameter recovery, interval reporting, or coverage; the
spatial sigma-slope intervals remain unsuitable for reporting; q=4
spatial correlations are derived-only and currently have basic
result-extraction and diagnostic checks 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 forms, 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:
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 |
| Do groups differ in expected response? | ordinary mu random intercepts and one-slope blocks |
When variance carries signal, Part 1 |
| Do groups differ in residual SD? | sigma ~ x + (1 | group) |
Which scale are you modelling? |
| Do groups differ in residual-scale slopes? |
sigma ~ x + (0 + w | group) for independent Gaussian
residual-scale slopes |
Which scale are you modelling? |
| 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? |
| Is the response continuous, count, bounded, robust, or zero-inflated? | fixed-effect non-Gaussian families | Choosing response families |
| Do repeated groups need non-Gaussian mean random effects? | selected ordinary non-Gaussian mu random intercepts and
independent numeric slopes |
Can I fit and report this model? |
| Are the large residuals changing the location and scale conclusions? | fixed-effect Student-t mu, sigma, and
nu formulas |
Robust continuous responses |
| 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_family() 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 |
| Does the response have known sampling variance or covariance? | Gaussian meta_V(V = V)
|
Mean effects and residual heterogeneity |
| Do species means remain similar after shared ancestry is included? |
phylo(1 | species, tree = tree) in univariate Gaussian
mu
|
Structural dependence overview |
| Are count means phylogenetically structured? |
phylo(1 | species, tree = tree) in ordinary Poisson or
ordinary NB2 mu
|
Can I fit and report this model? |
| 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 |
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:
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:
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.
drmTMB(
drm_formula(
y ~ x1 + phylo(1 | species, tree = tree),
sigma ~ x1
),
family = gaussian(),
data = dat
)
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 intervals remain unsuitable for reporting. 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 forms, 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:
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:
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 intervals remain unsuitable for reporting, while mesh/SPDE
beyond the fixed-kappa Gaussian mu intercept, additional
multiple or labelled Gaussian-slope layouts outside the exact documented
fitted forms, 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 documented fitted
forms, 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 numerical and simulation checks; 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? 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 forms, recovery-grade NB2 q1 spatial
sigma, Student-t spatial mu, Poisson spatial
zi, fixed-zi Poisson spatial mu,
and fixed-zi NB2 spatial mu forms 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:
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:
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:
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 |
| a worked Gaussian mean-variance analysis | When variance carries signal, Part 1 |
| a scale vocabulary check | Which scale are you modelling? |
| a robust continuous-response example | Robust continuous responses |
| a fixed-effect count example with overdispersion or structural zeros | Count abundance and extra zeros |
| a two-response residual-correlation example | Changing residual coupling with rho12 |
| known sampling variances or covariance | Mean effects and residual heterogeneity |
| implemented phylogenetic random intercepts or coordinate-spatial intercept/slope terms | Structural dependence overview |
| post-fit checks, predictions, residuals, and simulation | Checking and using fitted models |