Understanding a drmTMB fit with symbolizer

symbolizer does not just print equations for drmTMB; it makes a fitted distributional model auditable by connecting the R formula, symbolic model, parameter scales, assumptions, coefficient readings, and next diagnostic steps.

This vignette is the drmTMB-first tour. It assumes you have a drmTMB fit on the table and you want to know exactly what it claims, what it assumes, and what to look at next. If you want the package-wide tour first, read vignette("symbolizer"). The richer worked example with factors, transformations, and a multi-level random intercept now lives in Section 6 of this same vignette.

1. Why distributional models need a structured story

drmTMB fits Gaussian location-scale models, bivariate Gaussian location-scale-correlation models, and a family of non-Gaussian distributional models (Student-t, lognormal, Gamma, beta, beta-binomial, Poisson, nbinom2, truncated nbinom2, with optional zero-inflation / hurdle / cumulative-logit components). A univariate location-scale fit has two linear predictors: one for the mean mu, one for the residual standard deviation sigma. Both depend on covariates. Both have a link function. Both have their own coefficient vector.

That structure makes the question what does this fit actually say? genuinely hard to answer in one line. The mean coefficient on temperature lives on the response scale. The sigma coefficient on temperature lives on the log scale. The residual SD is no longer a single number; it is a vector with one entry per observation. A reader who only sees summary(fit) has to keep all of that in their head while reading the table.

symbolizer returns a symbolized_model: a structured object that connects the R formula to a symbolic equation, to assumptions, to per-coefficient readings on every relevant scale, and to a short list of diagnostic next steps. Every renderer in the package consumes that same object. No prose is generated by a language model at runtime; every phrase comes from a template in inst/extdata/.

Takeaway. A drmTMB fit has more structure than a single regression table. The symbolized_model is what you read instead of holding it all in your head.

2. From `drm_formula()` to equations

A small fit, by way of example. body_mass depends on temperature on both submodels:

library(symbolizer)
library(drmTMB)
#> 
#> Attaching package: 'drmTMB'
#> The following object is masked from 'package:base':
#> 
#>     beta

set.seed(20260523)
n <- 120
temperature <- runif(n, 10, 25)
dat <- data.frame(
  body_mass   = rnorm(n, 30 + 0.4 * temperature, exp(0.5 + 0.08 * temperature)),
  temperature = temperature
)

fit <- drmTMB(
  drm_formula(body_mass ~ temperature, sigma ~ temperature),
  family = gaussian(),
  data   = dat
)

The two arguments to drm_formula() map directly to two linear predictors. body_mass ~ temperature defines the mu submodel; sigma ~ temperature defines the sigma submodel. symbolize() reads the fit and returns the structured object:

sym <- symbolize(
  fit,
  symbols = c(body_mass = "W_i", temperature = "T_i"),
  units   = c(body_mass = "g",   temperature = "C"),
  context = "Gaussian location-scale body-size fit"
)

equations(sym) exposes the model as one row per renderable block. The notation = "both" argument controls how the result prints; both notations always live in the object.

equations(sym, notation = "both")

Index form:

\begin{aligned} W_i \mid \mu_i,\, \sigma_i \sim \mathrm{Normal}(\mu_i,\, \sigma_i^2) \\ \mu_i = \beta_{0} + \beta_{1} \, T_i \\ \log(\sigma_i) = \gamma_{0} + \gamma_{1} \, T_i \end{aligned}

Matrix form:

\begin{aligned} \mathbf{w} \mid \boldsymbol{\mu},\, \boldsymbol{\sigma} \sim \mathcal{N}(\boldsymbol{\mu},\, \mathrm{diag}(\boldsymbol{\sigma}^2)) \\ \boldsymbol{\mu} = \mathbf{X} \boldsymbol{\beta} \\ \log(\boldsymbol{\sigma}) = \mathbf{X}_{\sigma} \boldsymbol{\gamma} \end{aligned}

In math, the index form is

W_i \mid \mu_i,\, \sigma_i \sim \mathrm{Normal}(\mu_i,\, \sigma_i^2), \mu_i = \beta_0 + \beta_1\, T_i, \qquad \log(\sigma_i) = \gamma_0 + \gamma_1\, T_i,

and the matrix form is

\mathbf{w} \mid \boldsymbol{\mu},\, \boldsymbol{\sigma} \sim \mathcal{N}(\boldsymbol{\mu},\, \mathrm{diag}(\boldsymbol{\sigma}^2)), \boldsymbol{\mu} = \mathbf{X}\,\boldsymbol{\beta}, \qquad \log(\boldsymbol{\sigma}) = \mathbf{Z}\,\boldsymbol{\gamma}.

Each piece of the R call produced one piece of math. The left half of drm_formula() (body_mass ~ temperature) became the mu line; the right half (sigma ~ temperature) became the log-sigma line; the family = gaussian() choice became the distribution line and pinned the identity link on mu.

Takeaway. Each argument to drm_formula() produces one line of math. The mapping is one-to-one; nothing is implicit.

3. What sigma means in a location-scale model

In an introductory regression, sigma is a scalar — the residual SD, estimated once for the whole dataset. In a drmTMB location-scale fit it is a vector: each observation has its own predicted SD, given by a linear predictor on the log scale.

That choice is not cosmetic. It changes the units of every coefficient in the sigma submodel. drmTMB locks in the log link on sigma, so

\log(\sigma_i) = \gamma_0 + \gamma_1 T_i \quad\Longleftrightarrow\quad \sigma_i = \exp(\gamma_0)\,\exp(\gamma_1 T_i).

Two readings follow.

\gamma_0 is \log(\sigma) at the reference (here, T_i = 0). \exp(\gamma_0) is the residual SD itself at that reference value.
\gamma_1 is the change in \log(\sigma) per one-unit change in temperature. The biological reading is multiplicative: a one-unit change in temperature multiplies the residual SD by \exp(\gamma_1). If \gamma_1 > 0, hotter sites are more variable; if \gamma_1 < 0, they are less variable.

parameter_interpretation() reads each coefficient on every relevant scale at once:

parameter_interpretation(sym)

submodel	term_label	coefficient_role	estimate	95% CI	link_scale_reading	natural_scale_reading	variance_scale_reading	biological_reading
mu	(Intercept)	intercept	31.3	27.1, 35.4 *	Expected body_mass at the reference	Expected body_mass for the reference case	—	Baseline body_mass in the reference condition
mu	temperature	slope	0.311	0.0337, 0.588 *	Linear change in expected body_mass per unit of temperature	Expected body_mass changes by 0.311 per unit of temperature	—	A unit change in temperature shifts the expected body_mass by 0.311
sigma	(Intercept)	intercept	0.222	-0.285, 0.729	Log residual SD at the reference (SD = exp(0.222))	Residual SD = exp(0.222) at the reference	Residual variance = exp(2*0.222)	Baseline level of unexplained individual variation in body_mass
sigma	temperature	slope	0.0940	0.0662, 0.122 *	Log residual SD changes by 0.0940 per unit of temperature	Residual SD multiplied by exp(0.0940) per unit of temperature	Residual variance multiplied by exp(2*0.0940) per unit	A unit change in temperature multiplies the unexplained variability of body_mass by exp(0.0940)

Rows marked * have a 95% confidence interval that excludes zero (CI method: wald).

The variance_scale_reading column is the one that does the work for sigma coefficients: it states the multiplicative effect on the residual SD itself, with the fitted value of \exp(\hat{\gamma}_1) already substituted in.

Takeaway. Each observation has its own predicted sigma. The sigma coefficients live on the log scale; exp() of the slope is the multiplicative factor on the residual SD.

4. What is and isn’t assumed

A reader auditing the fit needs to know which assumptions are in the formula, which follow automatically from the parameterisation, and which are still the user’s responsibility. assumption_table() splits them into three buckets:

assumption_table(sym)

assumption	expression	biological meaning	status
conditional_distribution	W_i \mid \mu_i,\, \sigma_i \sim \mathrm{Normal}(\mu_i,\, \sigma_i^2)	body_mass varies normally around its expected value	explicit
linear_predictor	\mu_i = \beta_0 + \sum_k \beta_k X_{ki}	Expected body_mass is a linear combination of the mean-model predictors	explicit
linear_predictor	\log(\sigma_i) = \gamma_0 + \sum_k \gamma_k Z_{ki}	Log residual SD of body_mass is a linear combination of the scale-model predictors	explicit
independence	W_i \perp W_j \mid X \text{ for } i \ne j	Observations are conditionally independent given the predictors	follows from the formula
positivity	\sigma_i > 0	Residual SD is constrained positive via the log link	follows from the formula
no_missing_at_random	—	Observations are assumed not missing in a way that depends on the unobserved response	your responsibility

The three statuses map onto three honest questions.

stated — the formula declares this. The Gaussian distribution on W_i, the identity link on mu, the log link on sigma, and the linear-predictor structure of each submodel are all stated: they are what drm_formula() and family = gaussian() mean. There is nothing for the user to verify; this is what the model is.
implied — these follow from the parameterisation, not from data. The log link on sigma forces \sigma_i > 0 automatically; conditional independence given the predictors is the default of any fixed-effects regression. The package surfaces these so they are not invisible, but they are not separate empirical claims.
not_checked — these are the user’s responsibility. Missing-at- random, for example, cannot be inferred from the fit; symbolizer states it as an assumption the analyst still owns.

Three things symbolizer deliberately does not check. Residual normality on mu. Homoscedasticity within the sigma model (the model allows sigma to vary; whether the fitted sigma actually absorbs all the heteroscedasticity is a diagnostic question, not a symbolic one). Missing-at-random. The first two belong in section 6 (diagnostics); the third is documented in the assumption table itself.

Takeaway. assumption_table() is the audit lens. Stated and implied are guarantees of the parameterisation; not_checked is on the analyst.

5. Random intercepts and the assumption swap

When a random intercept enters the mu submodel, the independence statement automatically updates. symbolizer swaps the unconditional “observations are conditionally independent given the predictors” row for the conditional “observations are conditionally independent given the predictors and the random effects” row.

A small fit demonstrates:

set.seed(20260523)
n <- 120
n_sites <- 6
site <- factor(rep(letters[seq_len(n_sites)], length.out = n))
temperature <- runif(n, 10, 25)
site_re <- rnorm(n_sites, sd = 1.0)
dat_re <- data.frame(
  body_mass   = rnorm(
    n,
    30 + 0.4 * temperature + site_re[as.integer(site)],
    exp(0.5 + 0.05 * temperature)
  ),
  temperature = temperature,
  site        = site
)

fit_re <- drmTMB(
  drm_formula(body_mass ~ temperature + (1 | site), sigma ~ temperature),
  family = gaussian(),
  data   = dat_re
)

sym_re <- symbolize(
  fit_re,
  symbols = c(body_mass = "W_i", temperature = "T_i"),
  units   = c(body_mass = "g",   temperature = "C")
)

The mu equation picks up a + u_{\mathrm{site}(i)} term, and a new random-effect distribution row appears in the equation block:

equations(sym_re, notation = "both")

Index form:

\begin{aligned} W_i \mid \mu_i,\, \sigma_i \sim \mathrm{Normal}(\mu_i,\, \sigma_i^2) \\ \mu_i = \beta_{0} + \beta_{1} \, T_i + u_{site(i)} \\ \log(\sigma_i) = \gamma_{0} + \gamma_{1} \, T_i \\ u_{site} \sim \mathcal{N}(0,\, \sigma_{site}^2) \end{aligned}

Matrix form:

\begin{aligned} \mathbf{w} \mid \boldsymbol{\mu},\, \boldsymbol{\sigma} \sim \mathcal{N}(\boldsymbol{\mu},\, \mathrm{diag}(\boldsymbol{\sigma}^2)) \\ \boldsymbol{\mu} = \mathbf{X} \boldsymbol{\beta} + \mathbf{u} \\ \log(\boldsymbol{\sigma}) = \mathbf{X}_{\sigma} \boldsymbol{\gamma} \\ \mathbf{u}_{site} \sim \mathcal{N}(\mathbf{0},\, \sigma_{site}^2 \mathbf{I}_{6}) \end{aligned}

The assumption table swaps the independence row. Side by side:

# Without random effects
assumption_table(sym)[assumption_table(sym)$assumption %in%
                       c("independence", "independence_given_random_effects"), ]

assumption	expression	biological meaning	status
independence	W_i \perp W_j \mid X \text{ for } i \ne j	Observations are conditionally independent given the predictors	follows from the formula

# With (1 | site)
assumption_table(sym_re)[assumption_table(sym_re)$assumption %in%
                          c("independence", "independence_given_random_effects"), ]

assumption	expression	biological meaning	status
independence_given_random_effects	W_i \perp W_j \mid X\, \mathbf{u} \text{ for } i \ne j	Observations are conditionally independent given the predictors and the random effects	explicit

The structured object also exposes the random-effect structure as tidy tibbles, ready for downstream renderers:

sym_re$random_effects

submodel	term	group	levels	random effect	between-group SD
mu	`(1 \| site)`	site	6	u_{site(i)}	\sigma_{site}

sym_re$variance_components

parameter	sd_estimate	var_estimate
sigma_site_0	1.05	1.10

sd(site) is the between-site standard deviation; the symbol table reports it as \sigma_{\mathrm{site}}.

Takeaway. Adding (1 | group) automatically reworks the independence assumption from unconditional to conditional. The audit trail is built into the template, not into prose.

See the three views of this RE fit

The same as_html_three_views() widget that anchors the ladder article works on every symbolized_model. Calling it on sym_re lets a reader walk from the per-observation form y_i = \beta_0 + \beta_1 T_i + u_{\mathrm{site}(i)} + \varepsilon_i through the matrix abstraction \boldsymbol{\mu} = \mathbf{X}\boldsymbol{\beta} + \mathbf{Z}\mathbf{u} to the actual numeric arrays — with both the \mathbf{X}\boldsymbol{\beta} fixed-effect block and the \mathbf{Z}\mathbf{u} random-effect block populated from the data:

as_html_three_views(sym_re)

Skip three-views widget

What happens for each observation i – the per-individual reading.

Each observation is normally distributed around a mean that may shift with the fixed-effect predictors and a group offset, with a residual SD that may also shift with its own predictors – both location and spread are modeled.

Coefficient reading. On the response scale, \hat\beta is the additive change in the mean of the response for a one-unit increase in the predictor (identity link – no back-transformation needed).

\begin{aligned} W_i \mid \mu_i,\, \sigma_i & \sim \mathrm{Normal}(\mu_i,\, \sigma_i^2) \\ \mu_i & = \beta_{0} + \beta_{1} \, T_i + u_{site(i)} \\ \log(\sigma_i) & = \gamma_{0} + \gamma_{1} \, T_i \\ u_{site} & \sim \mathcal{N}(0,\, \sigma_{site}^2) \end{aligned}

where:

W_i — response variable \mathbb{R}^{120}
T_i — continuous predictor column of X (length 120)
\mu_i — conditional mu of body_mass \mathbb{R}^{120}
\sigma_i — residual standard deviation \mathbb{R}^{120}
\beta_{0}, \beta_{1} — mu submodel coefficients \mathbb{R}^{2}
\gamma_{0}, \gamma_{1} — sigma submodel coefficients \mathbb{R}^{2}
u_{site(i)} — random intercept by site scalar; \mathbb{R}^{6} in matrix form
\sigma_{site} — between-site standard deviation scalar

Where does the variation live? Where the variation lives – each row is one variance component (shown as a share of the total when a single total variance is defined).

site: variance = 1.10

ICC: ICC not available on this scale yet. (the residual SD varies across observations (a location-scale model), so there is no single within-group variance to form the ICC.)

Point estimates only; uncertainty not shown.

The same matrix equation, with your actual numbers stacked inside the brackets – what the computer multiplies. Showing first 5 and last 2 rows of n = 120.

For observation i = 1 of your data:

\begin{aligned} w_{1} &= \hat\beta_{0} + \hat\beta_{1}\,\mathrm{temperature}_{1} + \hat{u}_{\mathrm{site},\,\mathrm{a}} + \hat\varepsilon_{1} &\quad(\text{response equation, one row of the model}) \\ \hat\mu_{1} &= 29 + 0.458 \times 24.7 + (0.482) \approx 40.7 &\quad(\text{predicted mean} = \text{linear predictor}) \\ w_{1} &= \underbrace{40.7}_{\textstyle\,\hat\mu_{1}\,\text{(predicted)}\,} \;+\; \underbrace{(-0.526)}_{\textstyle\,\hat\varepsilon_{1}\,\text{(residual)}\,} &\quad(\text{observed} = \text{predicted mean} + \text{residual}) \end{aligned}

Stacking the same response equation for all n = 120 observations:

\underbrace{\begin{bmatrix} 40.2 \\ 41.4 \\ 44.2 \\ 33.8 \\ 38 \\ \vdots \\ 34.8 \\ 36.2 \end{bmatrix}}_{\textstyle\,\mathbf{w}_{\,120 \times 1}\;\text{(observed)}\,} \;=\; \underbrace{\begin{bmatrix} 1 & 24.7 \\ 1 & 18.9 \\ 1 & 20.6 \\ 1 & 16.3 \\ 1 & 21.1 \\ \vdots & \vdots \\ 1 & 11.9 \\ 1 & 18.9 \end{bmatrix}}_{\textstyle\,\mathbf{X}_{\,120 \times 2}\,}\, \underbrace{\begin{bmatrix} 29 \\ 0.458 \end{bmatrix}}_{\textstyle\,\hat{\boldsymbol{\beta}}_{\,2 \times 1}\;\text{(estimated)}\,} \;+\; \underbrace{\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix}}_{\textstyle\,\mathbf{Z}_{\text{site},\,\,120 \times 6}\,}\, \underbrace{\begin{bmatrix} 0.482 \\ 1.47 \\ -0.889 \\ -0.823 \\ -0.45 \\ 0.207 \end{bmatrix}}_{\textstyle\,\hat{\mathbf{u}}_{\text{site},\,6 \times 1}\;\text{(BLUP)}\,} \;+\; \underbrace{\begin{bmatrix} -0.526 \\ 2.37 \\ 6.67 \\ -1.77 \\ -0.109 \\ \vdots \\ 0.81 \\ -1.6 \end{bmatrix}}_{\textstyle\,\hat{\boldsymbol{\varepsilon}}_{\,120 \times 1}\;\text{(residual)}\,}

Left: observed vector \mathbf{w}. Middle: the prediction \mathbf{X}\hat{\boldsymbol{\beta}} + \mathbf{Z}\hat{\mathbf{u}} = \hat{\boldsymbol{\mu}}. Right: the residual vector \hat{\boldsymbol{\varepsilon}} = \mathbf{w} - \hat{\boldsymbol{\mu}}. Every row of this matrix equation is one of the response-equation rows from the worked row above.

Partial pooling. The random-effect estimates shown here (the BLUPs) are partially pooled: each group’s estimate is shrunk toward zero by an amount that grows when the group has little data and shrinks when the between-group variance is large. Groups with the least data are pulled hardest toward the overall mean.

And the \sigma submodel (no observed counterpart – \sigma’s job is to describe the spread of \hat{\boldsymbol{\varepsilon}}). For the same observation i = 1:

\begin{aligned} \log\hat\sigma_{1} &= \hat\gamma_{0} + \hat\gamma_{1}\,\mathrm{temperature}_{1} &\quad(\text{sigma submodel for observation 1, log link}) \\ \log\hat\sigma_{1} &= 0.578 + 0.0422 \times 24.7 = 1.62 &\quad(\text{with your numbers}) \\ \hat\sigma_{1} &= \exp(1.62) \approx 5.05 &\quad(\text{predicted residual SD for observation 1}) \end{aligned}

Stacking the same log-link equation for all n = 120 observations:

\log\!\underbrace{\begin{bmatrix} 5.05 \\ 3.95 \\ 4.25 \\ 3.55 \\ 4.34 \\ \vdots \\ 2.95 \\ 3.95 \end{bmatrix}}_{\textstyle\,\boldsymbol{\sigma}_{\,120 \times 1}\,} \;=\; \underbrace{\begin{bmatrix} 1 & 24.7 \\ 1 & 18.9 \\ 1 & 20.6 \\ 1 & 16.3 \\ 1 & 21.1 \\ \vdots & \vdots \\ 1 & 11.9 \\ 1 & 18.9 \end{bmatrix}}_{\textstyle\,\mathbf{X}_{\sigma,\,120 \times 2}\,}\, \underbrace{\begin{bmatrix} 0.578 \\ 0.0422 \end{bmatrix}}_{\textstyle\,\boldsymbol{\gamma}_{\,2 \times 1}\,}

6. A richer worked example: four predictor roles in one fit

The previous sections walked the basics on a simple two-predictor fit. A real model an ecologist would run usually has more going on: continuous predictors, log-transformed predictors, factor contrasts, and a random intercept on the mean. This section shows that the same symbolize() surfaces scale up — with two educator-facing surfaces (formula_bridge and notation_bridge) that get more interesting when the fit is more interesting.

set.seed(20260523)
nr           <- 200
temperature  <- runif(nr, 8, 30)
food         <- rlnorm(nr, meanlog = 1.5, sdlog = 0.6)
sex          <- factor(sample(c("female", "male"), nr, replace = TRUE))
site         <- factor(sample(letters[1:8], nr, replace = TRUE))

site_re <- rnorm(8, sd = 1.5)
mu_r    <- 5 + 0.5 * temperature + 2 * log(food) +
             3 * (sex == "male") + site_re[site]
sg_r    <- exp(0.3 + 0.05 * temperature + 0.1 * (sex == "male"))
dat_r   <- data.frame(
  body_mass = rnorm(nr, mu_r, sg_r),
  temperature = temperature, food = food, sex = sex, site = site
)

fit_rich <- drmTMB(
  drm_formula(
    body_mass ~ temperature + log(food) + sex + (1 | site),
    sigma ~ temperature + sex
  ),
  family = gaussian(),
  data   = dat_r
)
sym_rich <- symbolize(
  fit_rich,
  symbols = c(body_mass = "W_i", temperature = "T_i",
              food      = "F_i", sex         = "S_i"),
  units   = c(body_mass = "g",   temperature = "C", food = "g/day"),
  context = "structured body-size location-scale model"
)

formula_bridge() is the educator-facing translation table from R formula syntax to mathematics. Two rows (one per submodel), in both notations:

formula_bridge(sym_rich)

submodel	R syntax	meaning	math (index)	math (matrix)
mu	`body_mass ~ temperature + log(food) + sex + (1 \| site)`	Expected body_mass is a linear function of the mean-model predictors	\mu_i = \beta_{0} + \beta_{1} \, T_i + \beta_{2} \, \mathrm{log}(F_i) + \beta_{3} \, [sex = \mathrm{male}] + u_{site(i)}	\boldsymbol{\mu} = \mathbf{X} \boldsymbol{\beta} + \mathbf{u}
sigma	`sigma ~ temperature + sex`	Log residual SD of body_mass is a linear function of the scale-model predictors	\log(\sigma_i) = \gamma_{0} + \gamma_{1} \, T_i + \gamma_{2} \, [sex = \mathrm{male}]	\log(\boldsymbol{\sigma}) = \mathbf{X}_{\sigma} \boldsymbol{\gamma}

For the mu submodel, the R syntax body_mass ~ temperature + log(food) + sex + (1 | site) maps to the indexed equation in the mathematics column and to \boldsymbol{\mu} = \mathbf{X} \boldsymbol{\beta} + \mathbf{u} in the mathematics_matrix column. The random-intercept piece appears as + \mathbf{u} in the matrix form and as + u_{\mathrm{site}(i)} in the indexed form. A reader can compare the two columns and learn the translation without doing it by hand.

notation_bridge() is the deeper teaching surface: every symbol and every equation pairs its index form with its matrix form and tags both shapes:

notation_bridge(sym_rich)

concept	index	matrix	shape	concrete
conditional_distribution	W_i \mid \mu_i,\, \sigma_i \sim \mathrm{Normal}(\mu_i,\, \sigma_i^2)	\mathbf{w} \mid \boldsymbol{\mu},\, \boldsymbol{\sigma} \sim \mathcal{N}(\boldsymbol{\mu},\, \mathrm{diag}(\boldsymbol{\sigma}^2))	\mathbb{R}^n	\mathbb{R}^{200}
mu_linear_predictor	\mu_i = \beta_{0} + \beta_{1} \, T_i + \beta_{2} \, \mathrm{log}(F_i) + \beta_{3} \, [sex = \mathrm{male}] + u_{site(i)}	\boldsymbol{\mu} = \mathbf{X} \boldsymbol{\beta} + \mathbf{u}	\mathbb{R}^n	\mathbb{R}^{200}
sigma_linear_predictor	\log(\sigma_i) = \gamma_{0} + \gamma_{1} \, T_i + \gamma_{2} \, [sex = \mathrm{male}]	\log(\boldsymbol{\sigma}) = \mathbf{X}_{\sigma} \boldsymbol{\gamma}	\mathbb{R}^n	\mathbb{R}^{200}
body_mass	W_i	\mathbf{w}	\mathbb{R}^n	\mathbb{R}^{200}
parameter	\mu_i	\boldsymbol{\mu}	\mathbb{R}^n	\mathbb{R}^{200}
parameter	\sigma_i	\boldsymbol{\sigma}	\mathbb{R}^n	\mathbb{R}^{200}
coefficient	\beta_{0}, \beta_{1}, \beta_{2}, \beta_{3}	\boldsymbol{\beta}	\mathbb{R}^{p_\mu}	\mathbb{R}^{4}
coefficient	\gamma_{0}, \gamma_{1}, \gamma_{2}	\boldsymbol{\gamma}	\mathbb{R}^{p_\sigma}	\mathbb{R}^{3}
design_matrix	—	\mathbf{X}	\mathbb{R}^{n \times p_\mu}	\mathbb{R}^{200 \times 4}
design_matrix	—	\mathbf{X}_{\sigma}	\mathbb{R}^{n \times p_\sigma}	\mathbb{R}^{200 \times 3}
site	u_{site(i)}	\mathbf{u}_{site}	$scalar; \mathbb{R}^{G_{site}} in matrix form$	$scalar; \mathbb{R}^{8} in matrix form$
variance_component	\sigma_{site}	\sigma_{site}	scalar	scalar

The dimension column carries the shape rule (e.g. \mathbb{R}^n, \mathbb{R}^{n \times p_\mu}). The dimension_concrete column plugs in the actual numbers from this fit (e.g. \mathbb{R}^{200}, \mathbb{R}^{200 \times 4}). A reader learning matrix notation can stare at this table and see, in one place, that the \boldsymbol{\mu} vector has 200 entries, the \boldsymbol{\beta} vector has four, and \mathbf{X} is the 200 \times 4 matrix that connects them.

Takeaway. The same renderers (equations, formula_bridge, notation_bridge) scale from a two-predictor toy fit to a four-role mixed model. Abstract shape rules on one side, concrete numbers from this fit on the other.

7. Two responses at once: bivariate Gaussian

When you measure two responses on the same subject — paired growth and reproduction, two morphometric traits, a behaviour and its physiological correlate — a univariate fit fits each response in isolation and ignores the residual correlation between them. drmTMB’s biv_gaussian() family models both responses jointly, with a residual correlation \rho_{12} that itself can depend on covariates.

set.seed(20260524)
nb <- 80
x1 <- rnorm(nb)
x2 <- rnorm(nb)
# Cholesky factorisation for two correlated standard normals (rho = 0.6),
# no extra package needed:
e1 <- rnorm(nb)
e2 <- 0.6 * e1 + sqrt(1 - 0.6^2) * rnorm(nb)
dat_biv <- data.frame(
  growth = 30 + 1.5 * x1 + e1,
  repro  = 10 + 0.8 * x2 + e2,
  x1     = x1,
  x2     = x2
)

fit_biv <- drmTMB(
  drm_formula(
    mu1    = growth ~ x1,
    mu2    = repro  ~ x2,
    sigma1 = ~ 1,
    sigma2 = ~ 1,
    rho12  = ~ 1
  ),
  family = biv_gaussian(),
  data   = dat_biv
)

sym_biv <- symbolize(fit_biv)

The symbolized_model carries five submodels instead of two. Each appears in formula_bridge() with its own R syntax, statistical meaning, and symbolic form:

formula_bridge(sym_biv)

submodel	R syntax	meaning	math (index)	math (matrix)
mu1	`growth ~ x1`	Expected growth is a linear function of the first mean-model predictors	\mu_{1i} = \beta_{1,0} + \beta_{1,1} \, \mathrm{x1}_i	\boldsymbol{\mu}_{1} = \mathbf{X}_{1} \boldsymbol{\beta}_{1}
mu2	`repro ~ x2`	Expected repro is a linear function of the second mean-model predictors	\mu_{2i} = \beta_{2,0} + \beta_{2,1} \, \mathrm{x2}_i	\boldsymbol{\mu}_{2} = \mathbf{X}_{2} \boldsymbol{\beta}_{2}
sigma1	`~1`	Log residual SD of growth is a linear function of the first scale-model predictors	\log(\sigma_{1i}) = \gamma_{1,0}	\log(\boldsymbol{\sigma}_{1}) = \mathbf{X}_{\sigma_{1}} \boldsymbol{\gamma}_{1}
sigma2	`~1`	Log residual SD of repro is a linear function of the second scale-model predictors	\log(\sigma_{2i}) = \gamma_{2,0}	\log(\boldsymbol{\sigma}_{2}) = \mathbf{X}_{\sigma_{2}} \boldsymbol{\gamma}_{2}
rho12	`~1`	Fisher-z residual correlation between growth and repro is a linear function of the correlation-model predictors	\mathrm{tanh}^{-1}(\rho_{12,i}) = \rho_{0}	\mathrm{tanh}^{-1}(\boldsymbol{\rho}_{12}) = \mathbf{W} \boldsymbol{\rho}

mu1, mu2 are the two mean predictors; sigma1, sigma2 are the two log-SD predictors (one residual SD per response); rho12 is the residual correlation on the Fisher-z scale (so the linear predictor can range over all of \mathbb{R} while the correlation stays in (-1, 1)).

The joint conditional distribution is rendered in both index and matrix forms:

equations(sym_biv, notation = "both")

Index form:

\begin{aligned} (Y_{1i}, Y_{2i}) \mid \mu_{1i},\, \mu_{2i},\, \sigma_{1i},\, \sigma_{2i},\, \rho_{12,i} \sim \mathcal{N}_2\!\left((\mu_{1i}, \mu_{2i}),\, \Sigma_i\right) \\ \mu_{1i} = \beta_{1,0} + \beta_{1,1} \, \mathrm{x1}_i \\ \mu_{2i} = \beta_{2,0} + \beta_{2,1} \, \mathrm{x2}_i \\ \log(\sigma_{1i}) = \gamma_{1,0} \\ \log(\sigma_{2i}) = \gamma_{2,0} \\ \mathrm{tanh}^{-1}(\rho_{12,i}) = \rho_{0} \\ \Sigma_i = \begin{pmatrix}\sigma_{1i}^2 & \rho_{12,i}\sigma_{1i}\sigma_{2i} \\ \rho_{12,i}\sigma_{1i}\sigma_{2i} & \sigma_{2i}^2\end{pmatrix} \end{aligned}

Matrix form:

\begin{aligned} \mathbf{Y} \mid \mathbf{M},\, \mathbf{S},\, \boldsymbol{\rho}_{12} \sim \mathcal{N}_2(\mathbf{M},\, \boldsymbol{\Sigma}) \\ \boldsymbol{\mu}_{1} = \mathbf{X}_{1} \boldsymbol{\beta}_{1} \\ \boldsymbol{\mu}_{2} = \mathbf{X}_{2} \boldsymbol{\beta}_{2} \\ \log(\boldsymbol{\sigma}_{1}) = \mathbf{X}_{\sigma_{1}} \boldsymbol{\gamma}_{1} \\ \log(\boldsymbol{\sigma}_{2}) = \mathbf{X}_{\sigma_{2}} \boldsymbol{\gamma}_{2} \\ \mathrm{tanh}^{-1}(\boldsymbol{\rho}_{12}) = \mathbf{W} \boldsymbol{\rho} \\ \boldsymbol{\Sigma} = \mathrm{diag}(\boldsymbol{\sigma}_{1}) \mathbf{R}(\boldsymbol{\rho}_{12}) \mathrm{diag}(\boldsymbol{\sigma}_{2}) \end{aligned}

The matrix form is the cleaner reading: \mathbf{Y} is the n \times 2 response matrix, \mathbf{M} stacks the per-observation mean vectors, and \boldsymbol{\Sigma} is the per-observation 2 \times 2 covariance matrix with \sigma_{1,i}\,\sigma_{2,i}\,\rho_{12,i} on the off-diagonal.

The rho12 coefficient gets the same multi-scale treatment as every other coefficient. The reading row points the reader at the Fisher-z scale (where the coefficient lives), the natural correlation scale (tanh() of the link), and the biological meaning (residual co-variation between the two responses):

rho_rows <- parameter_interpretation(sym_biv)
rho_rows[rho_rows$submodel == "rho12", ]

submodel	term_label	coefficient_role	estimate	95% CI	link_scale_reading	natural_scale_reading	variance_scale_reading	biological_reading
rho12	(Intercept)	intercept	0.674	0.454, 0.894 *	Fisher-z residual correlation at the reference: tanh^{-1}(0.674)	Residual correlation rho12 = tanh(0.674) at the reference	Residual covariance at the reference is tanh(0.674) * sigma_1 * sigma_2	Two responses share co-variation with baseline correlation tanh(0.674)

Rows marked * have a 95% confidence interval that excludes zero (CI method: wald).

Takeaway. The same symbolized_model interface scales from two submodels to five. formula_bridge shows the new submodels, the joint distribution exposes \boldsymbol{\Sigma}, and parameter_interpretation gives the rho12 row the same multi-scale reading as every other coefficient.

8. What to inspect next

The point of the structured object is not just to render equations; it is to tell the reader what to look at next. model_card(sym) packages the full teaching bundle in one object, including a class- specific list of extraction calls. For a drmTMB fit, the extraction calls are ordered diagnostics-first:

mc <- model_card(sym)
mc$extraction_calls
#>                                                 Convergence and diagnostic checks 
#>                                                          "drmTMB::check_drm(fit)" 
#>                                                                         Residuals 
#>                                                                  "residuals(fit)" 
#>                                                                     Fitted values 
#>                                                                    "predict(fit)" 
#>                                                                Fixed effects (mu) 
#>                                               "drmTMB::fixef(fit, dpar = \"mu\")" 
#>                                                             Fixed effects (sigma) 
#>                                            "drmTMB::fixef(fit, dpar = \"sigma\")" 
#>                                                            Random effects (BLUPs) 
#>                                                              "drmTMB::ranef(fit)" 
#>                                                                    Symbolic story 
#>                                                      "symbolizer::symbolize(fit)" 
#>                                                                   Teaching bundle 
#>                              "symbolizer::model_card(symbolizer::symbolize(fit))" 
#>                                                 Group means (alongside contrasts) 
#>                                                    "symbolizer::group_means(sym)" 
#>                                         Per-group slopes (alongside interactions) 
#> "symbolizer::group_slopes(sym, continuous = <one of your continuous predictors>)" 
#>                                          Pairwise contrasts (which levels differ) 
#>                    "symbolizer::group_contrasts(sym, by = <one of your factors>)"

The ordering is intentional. A reader should know whether the fit is trustworthy before they interpret a single coefficient.

drmTMB::check_drm(fit) — convergence and identifiability. Did the optimiser converge? Are any parameters at the boundary? Are the standard errors finite?
residuals(fit) and predict(fit) — residual diagnostics. Plot residuals against the mu predictors (any structure left in the mean?) and absolute residuals against the sigma predictors (did the sigma submodel absorb the heteroscedasticity?).
drmTMB::ranef(fit) — random-effect BLUPs, when there are any. The distribution of BLUPs is worth a histogram or a QQ plot.
drmTMB::fixef(fit, dpar = "mu") and drmTMB::fixef(fit, dpar = "sigma") — the actual coefficient tables, by submodel. Read these after the diagnostics, not before.

model_card() also ships a set of recommended-plot recipes, in the same diagnostics-first order:

mc$recommended_plots
#>                                                                                     Residuals vs fitted 
#>                                                          "Pearson residuals vs fitted; check for trend" 
#>                                                                            Residuals on the sigma scale 
#> "abs residuals vs scale-model predictor; check that the sigma submodel captures the heteroscedasticity" 
#>                                                                                       Mean vs predictor 
#>                                     "plot fitted mu against the mean-model predictor; overlay raw data" 
#>                                                                                      Sigma vs predictor 
#>                                     "plot fitted sigma against the scale-model predictor; on log scale" 
#>                                                                              Random-effect distribution 
#>                                                                              "histogram or QQ of BLUPs"

These are text recipes only – symbolizer does not draw the plots itself. The recipes name the diagnostic plot the reader should make next, in their own plotting framework of choice.

Takeaway. Diagnostics before coefficients. The extraction-calls list is the order a reader should work through the fit.

9. Drafting the Methods paragraph

Once you’ve inspected the fit and you’re ready to write the paper, methods_text(sym) returns a draft Methods-section paragraph that matches what your model actually does. The prose is composed from a CSV template plus slots filled from the symbolized_model — there’s no LLM in the loop, so every phrase traces to the template row.

methods_text(sym)

We fit a Gaussian location-scale model to 120 observations of body_mass (units: g) using the drmTMB package (version 0.1.4) in R. The mean of body_mass was modelled as a linear function of temperature. The residual standard deviation was modelled on the log scale as a function of temperature. Parameters were estimated by maximum likelihood via the Template Model Builder framework. We report 95% Wald confidence intervals on the fixed-effect coefficients.

Reminders:

This is a draft. Confirm every claim against your fit before pasting.
Package version 0.1.4 was read from the local install; verify it matches what you actually fit with.

Two things to notice in the output:

The paragraph names the response, the predictors per submodel, the random-effect group (when present), the package and version, the number of observations, and the CI method. Each of those traces back to a slot you can inspect on methods_text(sym)$slots.
The print method appends Reminders before you paste — a short checklist that nudges you to confirm the package version, soften or sharpen the prose, and reword the random-effect sentence if you used a structure more complex than (1 | group).

The same call works for the richer fit from Section 6 (a 4-role location-scale model with a random intercept):

methods_text(sym_rich)

We fit a Gaussian location-scale model to 200 observations of body_mass (units: g) using the drmTMB package (version 0.1.4) in R. The mean of body_mass was modelled as a linear function of temperature, food, and sex. The residual standard deviation was modelled on the log scale as a function of temperature and sex. A random intercept for site was included in the model. Parameters were estimated by maximum likelihood via the Template Model Builder framework. We report 95% Wald confidence intervals on the fixed-effect coefficients.

Reminders:

This is a draft. Confirm every claim against your fit before pasting.
Package version 0.1.4 was read from the local install; verify it matches what you actually fit with.
Random-effect structure described here; if you used a more complex RE structure (random slopes, nested, crossed), reword that sentence.

and for the bivariate-Gaussian fit from Section 7:

methods_text(sym_biv)

We fit a bivariate Gaussian location-scale-correlation model to 80 paired observations of growth and repro using the drmTMB package (version 0.1.4) in R. The means of the two responses were modelled as linear functions of x1 and x2 respectively. The residual standard deviations of the two responses were modelled on the log scale, with both intercept-only. The residual correlation between the two responses was modelled with an intercept only on the Fisher-z linear predictor. Parameters were estimated by maximum likelihood via the Template Model Builder framework. We report 95% Wald confidence intervals on the fixed-effect coefficients.

Reminders:

This is a draft. Confirm every claim against your fit before pasting.
Package version 0.1.4 was read from the local install; verify it matches what you actually fit with.

The biv_gaussian template carries dedicated slots for mu1, mu2, sigma1, sigma2, and rho12 so the paragraph names both responses and reports the Fisher-z scale on the residual correlation.

Takeaway. methods_text(sym) produces a draft, not a final paragraph. Treat the output as starting material to edit; the reminders block is the editorial checklist.

10. The capability gate

symbolize() does not silently dispatch on every drmTMB fit. The capability registry gates each (class, family, component) tuple with one of five status words: Stable, First slice, Opt-in control, Planned or reserved, Unsupported or blocked.

caps <- symbolizer_capabilities()
caps[caps$class == "drmTMB", ]
#> # A tibble: 34 × 7
#>    class  family            component      status           since notes lives_in
#>    <chr>  <chr>             <chr>          <chr>            <chr> <chr> <chr>   
#>  1 drmTMB gaussian          mu             Stable           0.1.0 Univ… symboli…
#>  2 drmTMB gaussian          sigma          Stable           0.1.0 Univ… symboli…
#>  3 drmTMB gaussian          random_effects First slice      0.3.1 Rand… symboli…
#>  4 drmTMB gaussian          zi             Planned or rese… NA    Zero… symboli…
#>  5 drmTMB gaussian          hu             Planned or rese… NA    Hurd… symboli…
#>  6 drmTMB poisson           zi             First slice      0.4.0 Zero… symboli…
#>  7 drmTMB nbinom2           zi             First slice      0.4.0 Zero… symboli…
#>  8 drmTMB truncated_nbinom2 hu             First slice      0.4.0 Hurd… symboli…
#>  9 drmTMB gaussian          rho12          Planned or rese… NA    Biva… symboli…
#> 10 drmTMB student           mu             First slice      0.2.2 Stud… symboli…
#> # ℹ 24 more rows

The relevant rows are:

drmTMB / gaussian / mu — Stable
drmTMB / gaussian / sigma — Stable
drmTMB / gaussian / random_effects — First slice (intercept-only (1 | group))
drmTMB / biv_gaussian / mu1, mu2, sigma1, sigma2, rho12 — First slice (see Section 6)
Student-t, lognormal, Gamma, beta, beta_binomial, Poisson, nbinom2, truncated_nbinom2, cumulative_logit (with their zero-inflation and hurdle submodels) — First slice (the caps table above is the authoritative list)

If you call symbolize() on an unsupported tuple, the package errors with a friendly “Today we can read…” message that names the row in the registry. The registry, not the prose, decides what symbolize() accepts; check the status word before designing a workflow around a tuple.

Takeaway. The capability registry is the single source of truth. For drmTMB it covers the Gaussian location-scale family (including random intercepts and random slopes), the full bivariate Gaussian surface (mu1, mu2, sigma1, sigma2, rho12), the Student-t / lognormal / Gamma / beta / beta_binomial / Poisson / nbinom2 / truncated_nbinom2 / cumulative_logit families, and the zero-inflation and hurdle submodels. See vignette("symbolizer-roadmap") for the full matrix.

11. Where to read next

vignette("symbolizer") — the package-wide Get Started tour.
vignette("symbolizer-factors") — the categorical pedagogy guide: dummies, contrasts, interactions, and six common pitfalls.
vignette("symbolizer-gllvm") — the parallel story for latent-variable models in gllvmTMB.
vignette("symbolizer-compare") — comparing two symbolized models with compare_symbolic().

1. Why distributional models need a structured story

2. From drm_formula() to equations

3. What sigma means in a location-scale model

4. What is and isn’t assumed

5. Random intercepts and the assumption swap

See the three views of this RE fit

6. A richer worked example: four predictor roles in one fit

7. Two responses at once: bivariate Gaussian

8. What to inspect next

9. Drafting the Methods paragraph

10. The capability gate

11. Where to read next

2. From `drm_formula()` to equations