What can I fit today?

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.

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.

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
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
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.
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
Fixed-effect one-response families	Stable, with first count random-effect slices	Gaussian, Student-t, lognormal, Gamma, beta, beta-binomial, Poisson, NB2, zero-inflated Poisson, zero-inflated NB2, truncated NB2, hurdle NB2, cumulative-logit ordinal location, and ordinary Poisson/NB2 mu random intercepts plus independent slopes	ordinal scale/discrimination, zero-one-inflated beta or ordered beta, non-Gaussian sigma and shape random effects, correlated count slopes, zero-inflated count random effects, and random effects for most other non-Gaussian families	Wald fixed-effect intervals by default when sdreport() is computed; direct profile targets appear in profile_targets(), including Poisson and NB2 random-effect SDs through log_sd_mu
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, and multiple independent sigma slopes	correlated residual-scale slopes, coefficient-specific sd() slope models, and treating large q blocks as routine	check_drm() reports replication, weak-slope, boundary, and Hessian diagnostics; q > 2 mu block SDs and independent sigma slope SDs are direct profile targets, while q > 2 mu 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	Stable, with dense-storage guardrails	Gaussian meta_V(V = V) for diagonal, dense, and row-paired bivariate known sampling covariance, with 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
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
Ordinary bivariate covariance and corpairs()	First slice	matching labelled random intercepts in mu1/mu2, sigma1/sigma2, same-response mu/sigma, all-four q=4 blocks, and q=2 corpair(..., level = "group") ~ x	bivariate random slopes, full cross-parameter slope covariance, and slope1-slope2 plasticity-syndrome correlations	constant q=2 SD/correlation targets are direct; predictor-dependent corpair() values use newdata; q=4 rows are derived and mark derived intervals unavailable
Phylogenetic structured effects	First slices	univariate mu, bivariate mu1/mu2, labelled q=4 location-scale blocks, sd_phylo*() direct-SD surfaces, and q=2 phylogenetic corpair() regression	phylogenetic slopes, standalone or partial phylogenetic scale terms, structured rho12, and predictor-dependent q=4 phylogenetic correlations	direct phylogenetic SD and constant q=2 correlation targets are profile-ready; predictor-dependent corpair() values use newdata; q=4 correlations are derived-only for intervals
Coordinate spatial structured effects	First slice	spatial(1 | site, coords = coords) and one numeric spatial(1 + x | site, coords = coords) slope in univariate Gaussian mu	mesh/SPDE, multiple slopes, slope correlations, spatial sigma, bivariate spatial covariance, and spatial corpair()	sdpars$mu, ranef("spatial_mu"), profile_targets(), and check_drm() expose the fitted coordinate fields
Animal and lower-level relatedness markers	Planned markers	animal() and relmat() are exported, documented, and parsed so examples and roadmap text use the intended public grammar	fitted pedigree-derived animal models, direct A or Ainv relatedness inputs, low-level K or Q matrix fits, and their diagnostics	no fitted likelihood or interval targets yet; dense K/A examples will be small-to-moderate, while scalable claims require sparse Q/Ainv precision evidence
Profile intervals and diagnostics	First slice	fixed effects, direct SD/correlation targets, row-specific sigma, sigma1, sigma2, rho12, and fitted q=2 corpair() values	automatic intervals for every derived covariance summary and public bootstrap intervals	conf.status, profile.boundary, and profile.message expose interval status; q=4 rows use derived_interval_unavailable, and bootstrap requests error before interval work begins
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, phylogenetic slopes, mesh/SPDE, spatial corpair(), bivariate random slopes, 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
Do groups differ in expected response?	ordinary mu random intercepts and one-slope blocks	When variance carries signal
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
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 a strict continuous proportion or successes out of known trials?	fixed-effect beta() or beta_binomial()	Proportions and success rates
Does the response have known sampling variance or covariance?	Gaussian meta_known_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
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

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)
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 are the phylogenetic and coordinate-spatial Gaussian pieces. The animal and relmat() markers are exported and documented as the planned public grammar, but drmTMB() still rejects them before fitting until their likelihood, diagnostics, profile targets, and recovery tests exist.

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 residual sigma remains the within-observation SD. 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() slice is univariate only: spatial(1 | site, coords = coords) fits a spatial location intercept, and spatial(1 + x | site, coords = coords) fits independent spatial intercept and numeric slope fields. Mesh/SPDE spatial fields, spatial bivariate covariance blocks, spatial direct-SD models, multiple spatial slopes, slope correlations, 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, use the documented planned names in notes and design examples, but do not treat them as runnable analysis code yet:

# Planned, not fitted yet:
drmTMB(
  drm_formula(
    body_size ~ sex + animal(1 | individual, pedigree = pedigree),
    sigma ~ cohort
  ),
  family = gaussian(),
  data = dat
)

The fitted fallback is an ordinary Gaussian random-effect model:

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 until the animal-model likelihood, diagnostics, profile targets, and recovery tests are implemented.

Use relmat(1 | id, K = K) or relmat(1 | id, Q = Q) only for future 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 future latent random-effect relatedness or precision matrices.

For fitted 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 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
correlated residual-scale slope blocks in sigma	needs covariance recovery for scale intercept-slope blocks before public use
larger cross-formula labelled covariance blocks linking mu and sigma	extends beyond the first fitted random-intercept mean-scale correlations
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
phylogenetic slopes, predictor-dependent q=4 phylogenetic correlations, and richer phylogenetic covariance blocks	the first bivariate mu1/mu2 phylogenetic location slice, the q=2 predictor-dependent location-location slice, and the first constant q=4 phylogenetic location-scale block are fitted
richer spatial random effects	first univariate Gaussian coordinate-based mu intercept and one-slope paths implemented; mesh/SPDE, spatial scale, multiple slopes, slope correlations, bivariate spatial covariance, and spatial corpair() paths remain planned
animal models and generic relatedness matrices	animal() and relmat() are documented markers only until pedigree or known-matrix likelihoods, diagnostics, profile targets, and recovery tests are implemented
skew-normal residual asymmetry	skew_normal() with nu ~ ... is planned only; use Gaussian or fitted Student-t sensitivity checks until density, normal-limit, skew-recovery, and false-positive heteroscedasticity tests exist
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"), sdpars$mu, and check_drm() for the first univariate coordinate slice	coordinate-structured site intercept or one numeric slope deviations in one response; spatial correlations are 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; use confint() on the explicit cor:phylo: target when a direct profile interval is worth the extra compute.

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
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
post-fit checks, predictions, residuals, and simulation	Checking and using fitted models