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, 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.
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 for structural exact 0/1 outcomes |
probability-component fixed effects, with the exact relatedness hurdle route used only to confirm fit/extractor feasibility | random effects or structured dependence in zi,
zoi, or coi, 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 |
| 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 |
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 check the roadmap 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-0227, mc-0464, mc-0539, and
mc-0575 slopes plus mc-0382 lognormal
sigma-intercept domain are inference-ready with caveats under their
separate designs |
| 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, 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 |
sdpars$mu, sdpars$sigma, marker-specific
ranef() blocks, 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, 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 |
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 |
Implementation map |
| 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() 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
|
Implementation map |
| 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) |
sdpars$mu, ranef(fit, "mu")
|
groups differ in the x slope |
(0 + x | id) |
sdpars$mu, profile_targets(fit)
|
| group baselines and slopes are correlated |
(1 + x | id) or (1 + x | p | id)
|
corpars$mu,
corpairs(fit, class = "mean-slope")
|
| groups differ in residual scale | sigma ~ z + (1 | id) |
sdpars$sigma, sigma(fit)
|
| groups differ in residual-scale slopes | sigma ~ z + (0 + w | id) |
sdpars$sigma, 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, in the future, mesh/SPDE
structure. 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 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 now flags near-boundary corpars$phylo
values 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 interval gate remains blocked, while mesh/SPDE, 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 check the roadmap or design notes 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 |
corpars$mu or group-level rows from
corpairs(fit)
|
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"),
ranef(fit, "spatial_sigma"), sdpars$mu,
sdpars$sigma, 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:
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 |