Skip to contents

Equations via symbolizer

Source: vignettes/equations-via-symbolizer.Rmd
equations-via-symbolizer.Rmd

drmSEM produces a graph of distributional models: each endogenous node is one drmTMB fit, and paths into mu / sigma / nu / zi / hu / rho12 are first-class. Reporting that kind of model in a paper requires writing the per-node equations, families, links, and assumptions — error-prone work that also goes stale when the code changes.

symbolizer renders the equations from the fitted object itself. Its drmTMB extractor covers exactly the families and components drmSEM uses, so wiring it up is a thin bridge: walk the SEM in topological order, symbolize each node’s drmTMB fit, and collate.

A two-node example 

sem <- drm_sem(
  size      = drm_node(drmTMB::bf(size ~ temp + habitat, sigma ~ temp),
                       family = stats::gaussian()),
  abundance = drm_node(drmTMB::bf(abundance ~ size + temp, zi ~ habitat),
                       family = drmTMB::nbinom2()),
  data = dat
)
#> ℹ Fitting node "size"
#> ✔ Fitting node "abundance" [121ms]
#> 
#> ℹ Fitting node "abundance"
#> ✔ Fitting node "abundance" [93ms]
#>

Render 

sym <- symbolizer::symbolize(sem) # walks sem$records, symbolizes each node
sym                          # prints a summary header
#> <symbolized_drm_sem: 2 nodes, 6 edges>
#> Nodes (topological order): "size" and "abundance"
#> Use `symbolizer::as_latex()` / `symbolizer::equations()` /
#> `symbolizer::assumption_table()` to render.

symbolizer::as_latex(sym)    # publication-quality equations, node-by-node
#> [1] "%% Node: size\n\\begin{aligned}\n\\mathrm{size}_i \\mid \\mu_i,\\, \\sigma_i & \\sim \\mathrm{Normal}(\\mu_i,\\, \\sigma_i^2) \\\\\n\\mu_i & = \\beta_{0} + \\beta_{1} \\, \\mathrm{temp}_i + \\beta_{2} \\, [habitat = \\mathrm{B}] \\\\\n\\log(\\sigma_i) & = \\gamma_{0} + \\gamma_{1} \\, \\mathrm{temp}_i\n\\end{aligned}\n\n%% Node: abundance\n\\begin{aligned}\n\\mathrm{abundance}_i \\mid \\mu_i,\\, \\sigma_i & \\sim \\mathrm{NegBin}(\\mu_i,\\, \\mathrm{size} = \\exp(\\sigma_i)); \\mathrm{Var}(\\mathrm{abundance}_i) = \\mu_i + \\mu_i^2 / \\exp(\\sigma_i) \\\\\n\\log(\\mu_i) & = \\beta_{0} + \\beta_{1} \\, \\mathrm{size}_i + \\beta_{2} \\, \\mathrm{temp}_i \\\\\n\\mathrm{logit}(\\pi_{\\mathrm{zi},\\, i}) & = \\alpha_{0} + \\alpha_{1} \\, [habitat = \\mathrm{B}]\n\\end{aligned}"
symbolizer::equations(sym)   # tidy tibble: one row per component, keyed by node

sizei∣μi,σi∼Normal(μi,σi2)μi=β0+β1tempi+β2[habitat=B]log⁡(σi)=γ0+γ1tempiabundancei∣μi,σi∼NegBin(μi,size=exp⁡(σi));Var(abundancei)=μi+μi2/exp⁡(σi)log⁡(μi)=β0+β1sizei+β2tempilogit(πzi,i)=α0+α1[habitat=B] \begin{aligned} \mathrm{size}_i \mid \mu_i,\, \sigma_i \sim \mathrm{Normal}(\mu_i,\, \sigma_i^2) \\ \mu_i = \beta_{0} + \beta_{1} \, \mathrm{temp}_i + \beta_{2} \, [habitat = \mathrm{B}] \\ \log(\sigma_i) = \gamma_{0} + \gamma_{1} \, \mathrm{temp}_i \\ \mathrm{abundance}_i \mid \mu_i,\, \sigma_i \sim \mathrm{NegBin}(\mu_i,\, \mathrm{size} = \exp(\sigma_i)); \mathrm{Var}(\mathrm{abundance}_i) = \mu_i + \mu_i^2 / \exp(\sigma_i) \\ \log(\mu_i) = \beta_{0} + \beta_{1} \, \mathrm{size}_i + \beta_{2} \, \mathrm{temp}_i \\ \mathrm{logit}(\pi_{\mathrm{zi},\, i}) = \alpha_{0} + \alpha_{1} \, [habitat = \mathrm{B}] \end{aligned} 

symbolizer::assumption_table(sym) # families, links, deferred components
assumption expression biological meaning status
conditional_distribution sizei∣μi,σi∼Normal(μi,σi2)\mathrm{size}_i \mid \mu_i,\, \sigma_i \sim \mathrm{Normal}(\mu_i,\, \sigma_i^2) size varies normally around its expected value explicit
linear_predictor μi=β0+∑kβkXki\mu_i = \beta_0 + \sum_k \beta_k X_{ki} Expected size is a linear combination of the mean-model predictors explicit
linear_predictor log⁡(σi)=γ0+∑kγkZki\log(\sigma_i) = \gamma_0 + \sum_k \gamma_k Z_{ki} Log residual SD of size is a linear combination of the scale-model predictors explicit
independence sizei⟂sizej∣X for i≠j\mathrm{size}_i \perp \mathrm{size}_j \mid X \text{ for } i \ne j Observations are conditionally independent given the predictors follows from the formula
positivity σi>0\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
conditional_distribution abundancei∣μiσi∼NegBin(μisize=exp⁡(σi))\mathrm{abundance}_i \mid \mu_i\, \sigma_i \sim \mathrm{NegBin}(\mu_i\, \mathrm{size}=\exp(\sigma_i)) abundance is a non-negative count with overdispersion; mean = mu_i; Var(abundance) = mu_i + mu_i^2 / exp(sigma_i) explicit
non_negative_integer_response abundancei∈{0$|1|2||linearpredictor|\mathrm{abundance}_i \in \{0\$ |1\ |2\ | |linear_predictor |(_i) = _0 + _k k X{ki}$ log link on mu; mu_i is the expected count of abundance; a unit change in a predictor multiplies the expected count by exp(beta) explicit
linear_predictor log⁡(σi)=γ0+∑kγkZki\log(\sigma_i) = \gamma_0 + \sum_k \gamma_k Z_{ki} sigma_i parameterises overdispersion on the log scale; small exp(sigma_i) (close to zero) means highly overdispersed, very large exp(sigma_i) approaches Poisson explicit
independence abundancei⟂abundancej∣X for i≠j\mathrm{abundance}_i \perp \mathrm{abundance}_j \mid X \text{ for } i \ne j Observations are conditionally independent given the predictors follows from the formula
positivity σi>0 unconstrained on log scale\sigma_i > 0 \text{ unconstrained on log scale} σi\sigma_i is unconstrained on the log scale; exp⁡(σi)\exp(\sigma_i) is the size parameter (>0> 0) follows from the formula
no_missing_at_random — Observations are assumed not missing in a way that depends on the unobserved response your responsibility
zero_inflation_mixture P(abundancei=0∣𝐗i)=πzii+(1−πzii)⋅NegBin(0;μiexp⁡(σi))P(\mathrm{abundance}_i = 0 \mid \mathbf{X}_i) = \pi_{\mathrm{zi}\, i} + (1 - \pi_{\mathrm{zi}\, i}) \cdot \mathrm{NegBin}(0; \mu_i\, \exp(\sigma_i)) Mixture of a structural-zero component and the count distribution explicit
linear_predictor logit(πzi$|i)=α0+∑kαkVki|logitlinkonthezero−inflationprobability||hurdletwopart|\mathrm{logit}(\pi_{\mathrm{zi}\$ |i}) = \alpha_0 + \sum_k \alpha_k V_{ki} |logit link on the zero-inflation probability | |hurdle_two_part |P(i = 0) = {, i}; P(_i = k _i > 0) = ^{+}(k; _i, (_i))$ Hurdle structure on nbinom2: logit for zero vs non-zero, then truncated nbinom2 for positives explicit
linear_predictor $(_{$ i}) = _0 + _k k V{ki} logit link on the hurdle probability

The collator labels each block with ## Node: <name> so the reader can see exactly which equations belong to which response. sigma and zi sub-models appear next to their parent mu, with the link function spelled out — which makes the component-labelled grammar of drmSEM concrete instead of abstract.

Scope

symbolize.drm_sem() covers the model specification layer (per-node distributional equations, families, links, assumptions). The SEM layer — the DAG, [dsep()], [fisher_c()], and the simulation-based effect calculus — stays with drmSEM itself, because the distribution-mediated effect row of [indirect_effects()] has no closed-form symbolic representation by construction. Symbolizer correctly stops at the per-node equations and lets drmSEM do the propagation.

See also

References