---
title: "Known-matrix relatedness with relmat"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Known-matrix relatedness with relmat}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.1,
  dpi = 144,
  dev.args = list(bg = "white")
)
library(drmTMB)

relmat_guide_theme <- function() {
  ggplot2::theme_minimal(base_size = 11) +
    ggplot2::theme(
      panel.grid.minor = ggplot2::element_blank(),
      axis.title = ggplot2::element_text(colour = "grey15"),
      axis.text = ggplot2::element_text(colour = "grey25"),
      plot.title = ggplot2::element_text(
        face = "bold",
        colour = "grey10",
        margin = ggplot2::margin(b = 4)
      ),
      plot.subtitle = ggplot2::element_text(
        colour = "grey30",
        margin = ggplot2::margin(b = 8)
      ),
      plot.background = ggplot2::element_rect(fill = "white", colour = NA),
      panel.background = ggplot2::element_rect(fill = "white", colour = NA),
      legend.position = "bottom"
    )
}

relmat_eye_theme <- function() {
  relmat_guide_theme() +
    ggplot2::theme(
      panel.grid.major.y = ggplot2::element_blank(),
      legend.position = "none"
    )
}

relatedness_heatmap_data <- function(x) {
  out <- as.data.frame(as.table(x), stringsAsFactors = FALSE)
  names(out) <- c("row", "column", "relatedness")
  out$row <- factor(out$row, levels = rev(rownames(x)))
  out$column <- factor(out$column, levels = colnames(x))
  out
}

simulate_relmat_guide_data <- function(seed = 20260524) {
  set.seed(seed)
  n_line <- 24
  n_each <- 5
  line_levels <- paste0("line_", seq_len(n_line))
  line <- factor(rep(line_levels, each = n_each))
  temperature <- runif(length(line), -1, 1)
  treatment <- factor(rep(rep(c("control", "heated"), each = n_each), n_line / 2))

  K <- outer(seq_len(n_line), seq_len(n_line), function(i, j) {
    0.36^abs(i - j)
  })
  diag(K) <- diag(K) + 0.08
  dimnames(K) <- list(line_levels, line_levels)
  Q <- solve(K)

  line_effect <- as.vector(t(chol(K)) %*% rnorm(n_line)) * 0.42
  names(line_effect) <- line_levels
  sigma <- exp(-1.20 + 0.10 * (treatment == "heated"))
  seed_mass <- 2.2 + 0.50 * temperature + 0.25 * (treatment == "heated") +
    line_effect[line] + rnorm(length(line), sd = sigma)

  list(
    data = data.frame(
      seed_mass = seed_mass,
      temperature = temperature,
      treatment = treatment,
      line = line
    ),
    K = K,
    Q = Q
  )
}

simulate_relmat_q2_guide_data <- function(seed = 1) {
  set.seed(seed)
  n_observation_block <- 40
  n_line <- 10
  n_each <- 4
  observation <- factor(rep(paste0("obs_", seq_len(n_observation_block)), each = n_each))
  sex <- factor(rep(rep(c("female", "male"), each = n_each), n_observation_block / 2))
  age <- runif(length(observation), 0, 2)
  line_levels <- paste0("line_", seq_len(n_line))
  line <- factor(
    rep(rep(line_levels, each = 4), each = n_each),
    levels = line_levels
  )

  K <- outer(seq_len(n_line), seq_len(n_line), function(i, j) {
    0.25^abs(i - j)
  })
  diag(K) <- diag(K) + 0.10
  dimnames(K) <- list(line_levels, line_levels)
  Q <- solve(K)

  z_mass <- rnorm(n_line)
  z_height <- 0.55 * z_mass + sqrt(1 - 0.55^2) * rnorm(n_line)
  u_mass <- as.vector(t(chol(K)) %*% z_mass) * 0.35
  u_height <- as.vector(t(chol(K)) %*% z_height) * 0.30
  names(u_mass) <- line_levels
  names(u_height) <- line_levels

  e_mass <- rnorm(length(line))
  e_height <- 0.10 * e_mass + sqrt(1 - 0.10^2) * rnorm(length(line))
  sigma_mass <- exp(-1.0 + 0.15 * (sex == "male"))
  sigma_height <- exp(-1.1 - 0.10 * (sex == "male"))

  list(
    data = data.frame(
      line = line,
      age = age,
      sex = sex,
      seed_mass = 2.1 + 0.45 * age + 0.20 * (sex == "male") +
        u_mass[line] + sigma_mass * e_mass,
      plant_height = 1.4 + 0.30 * age - 0.12 * (sex == "male") +
        u_height[line] + sigma_height * e_height
    ),
    K = K,
    Q = Q
  )
}
```

<style>
@media (max-width: 575.98px) {
  main h1 {
    hyphens: none;
    overflow-wrap: normal;
    word-break: normal;
  }

  main table {
    display: block;
    min-width: 720px;
    max-width: 100%;
    overflow-x: auto;
  }
}
</style>

Use `relmat()` when the dependence is among latent group-level deviations and
the matrix is already known. The matrix can be a covariance, a correlation, or
an inverse covariance/precision matrix. The common case is a correlation-like
matrix among lines, strains, plots, populations, or experimental units:

```text
u_id ~ N(0, sd_relmat^2 K)
y_i  = mu_i + u_id[i] + error_i
```

If the natural input is a precision matrix, use `Q` instead:

```text
u_id ~ N(0, sd_relmat^2 Q^-1)
```

This is a higher-level random-effect structure. It is not the same as known
sampling covariance among observed effect-size estimates; that observation-level
matrix belongs to the meta-analysis route with `meta_V(V = V)`.

## When is `relmat()` the right tool?

Use `relmat()` only after the named biological or spatial route is not the
better description. It is the general fallback for a matrix you already have,
so reach for it once the [structural-dependence
overview](structural-dependence.html) has ruled out the named routes.

| Situation | Use | Why |
| --- | --- | --- |
| Genomic or marker-based similarity among inbred lines, strains, or cultivars | `relmat(1 | line, K = G)` | `G` is often a correlation or genomic relationship matrix. The fitted SD turns it into the latent covariance `sd_relmat^2 G`. |
| A lab, assay, or ecological kernel among experimental units | `relmat(1 | unit, K = K_unit)` | The matrix says which unit-level deviations should be similar, but it is not a pedigree, species tree, or coordinate-spatial surface. |
| A precomputed inverse relatedness matrix from another tool | `relmat(1 | id, Q = Q_id)` | The natural object is already precision-like, so `drmTMB` does not need to invert a dense covariance matrix. |
| A graph, network, river, areal, or custom Gaussian Markov random-field precision that was built and checked outside `drmTMB` | `relmat(1 | node, Q = Q_node)` | The scientific structure is a latent node-level field with a known precision, but it is not the current coordinate `spatial()` route. |

Do not use `relmat()` just because you have any covariance matrix. If the
matrix is known sampling covariance for observed estimates, use
`meta_V(V = V)`. If the matrix is additive relatedness from a pedigree or animal
model, use `animal()`. If it comes from an ultrametric species tree, use
`phylo()`. If coordinates define the dependence, use `spatial()`.

## Choosing `K` or `Q`

Use `K` when the matrix is easiest to interpret as relatedness, covariance, or
correlation. A correlation matrix with diagonal 1 is usually the clearest input:
the fitted `sd_relmat` then controls the variance scale of the latent
deviations. A covariance matrix with non-unit diagonal is also allowed, but then
the diagonal scale is part of the known structure.

Use `Q` when the matrix is easiest to interpret as an inverse covariance or
precision. This is common when another model, pedigree engine, graph model, or
Gaussian Markov random-field construction gives a sparse precision matrix. In
that case the latent covariance is proportional to `solve(Q)`, but the model can
work with `Q` directly.

## Two concrete examples

Suppose `G` is a marker-derived genomic relationship matrix among experimental
lines. The scientific question is whether lines that are genomically similar
also have similar unexplained deviations in seed mass:

```r
fit_grm <- drmTMB(
  bf(seed_mass ~ temperature + relmat(1 | line, K = G), sigma ~ temperature),
  data = dat,
  family = gaussian()
)
```

Here `G` can be a correlation matrix. The `sd:mu:relmat(1 | line)` row of
`summary(fit_grm)$parameters` reports the fitted structured SD scale `s`, so
the latent covariance is `s^2 G` and line `i` has
marginal SD `s sqrt(G[i, i])`. Thus `s` is a line-level marginal SD only when
the corresponding diagonal entry of `G` is one. This is not known sampling
error; it is latent among-line structure left after the fixed effects.

Now suppose `Q_river` is a positive-definite precision matrix for sites on a
river network, built outside `drmTMB` from an adjacency or flow model. The
scientific question is whether nearby or connected sites have similar latent
condition deviations:

```r
fit_river <- drmTMB(
  bf(condition ~ treatment + relmat(1 | reach, Q = Q_river), sigma ~ 1),
  data = dat,
  family = gaussian()
)
```

This is a good `Q` example because the natural object is already an inverse
covariance. If `C = solve(Q_river)`, the model estimates the fitted structured
scale `s`, the latent covariance is `s^2 C`, and reach `i` has marginal SD
`s sqrt(C[i, i])`. The known precision controls the pattern of similarity
across reaches; `s` is not a common site-level marginal SD unless `diag(C)` is
one.

## What is fitted today

| Question | Syntax | Status |
| --- | --- | --- |
| Does one Gaussian response have deviations structured by a user-supplied latent covariance or precision matrix? | `relmat(1 | id, K = K)` or `relmat(1 | id, Q = Q)` in `mu` | Fitted first slice. `K` is a covariance matrix; `Q` is a precision matrix. |
| Does one numeric fixed-effect slope also have known-matrix deviations? | `relmat(1 + x | id, K = K)` or `relmat(1 + x | id, Q = Q)` in `mu` | Fitted first univariate one-slope slice with independent intercept and slope fields. Additional structured-slope layouts outside the exact fitted bivariate ledger cells and intercept-slope correlations remain planned. |
| Does a one-response Gaussian model have known-matrix residual-scale deviations? | `relmat(1 | id, Q = Q)` in `sigma`, alone or matching the `mu` term; `relmat(1 + x | id, K = K)` for the exact one-slope route | The residual-scale intercept is fitted. The exact K/Q q1 `sigma` one-slope route is fitted and inference-ready with caveats under raw uncorrected log-SD Wald-z intervals; profile is diagnostic-only at `g = 8`. Matching univariate `mu` and `sigma` `relmat()` intercepts estimate one known-matrix mean-scale correlation. Broader K/Q bridge claims, additional multiple or labelled sigma-slope layouts outside the exact fitted bivariate ledger cells, and structured slope correlations remain planned. |
| Do two Gaussian response means share a known-matrix latent correlation? | matching labelled `relmat(1 | p | id, K = K)` or `relmat(1 | p | id, Q = Q)` terms in `mu1` and `mu2` | ML fits both representations. Native REML has `point_fit_recovery` evidence only for the exact supplied-`K`, location-intercept cell with constant residual formulas, complete pairs, unit weights, and no other model layer. `corpairs(level = "relmat")` reports the latent relatedness correlation separately from residual `rho12`. `Q` remains ML-only for this bivariate route. |
| Does the known-matrix layer link means and residual scales across two responses? | the same labelled `relmat(1 | p | id, Q = Q)` term in `mu1`, `mu2`, `sigma1`, and `sigma2` | Fitted first q=4 location-scale slice. `corpairs(level = "relmat")` reports six constant latent relatedness correlations; profile intervals for those derived correlations remain unavailable. |
| Does the known-matrix layer change multiple slopes or a predictor-dependent correlation? | examples such as `relmat(1 + x + z | id, K = K)` or `relmat()` `corpair()` formulas | Planned. These routes need matrix validation, diagnostics, recovery tests, interval targets, and examples before use. |

Use `animal()` when the matrix is additive relatedness for individual animal
models. Use `phylo()` when the matrix comes from a species tree and the tree
route is available. Use `spatial()` when the structure is induced by sampling
coordinates. Use `relmat()` for validated latent matrices that do not belong to
those named routes.

## Start with the smallest useful model

For one response, start with the fitted location-intercept route:

```r
fit_relmat <- drmTMB(
  y ~ treatment + relmat(1 | line, K = K),
  data = dat,
  family = gaussian()
)
```

For the exact Arc 1a REML route, keep `sigma ~ 1`, use an unlabelled intercept
or independent intercept-plus-one-numeric-slope shape, and set `REML = TRUE`:

```r
fit_relmat_reml <- drmTMB(
  bf(
    y ~ x + relmat(1 + x | line, K = K),
    sigma ~ 1
  ),
  data = dat,
  family = gaussian(),
  REML = TRUE
)
```

The multi-seed campaign used the `K` representation shown here, with
`n_each = 20` and exactly `M = {8, 16, 32}` relatedness levels. `Q` has
deterministic representation-parity evidence only, not multi-seed campaign
coverage. Those Arc 1a domains do not admit slope-only, labelled,
multiple-slope, scale-side, bivariate, or non-Gaussian REML routes. The
separate exact bivariate supplied-`K` exception is shown next.

For two responses under REML, build one named covariance matrix whose rows and
columns exactly match the factor levels, then use the same matrix and block
label in `mu1` and `mu2`. This complete example constructs the data and a dense
correlation matrix rather than leaving either object undefined:

```r
set.seed(20260715)
line_levels <- paste0("line_", seq_len(10))
line <- factor(rep(line_levels, each = 6), levels = line_levels)
treatment <- rep(rep(c(0, 1), each = 3), length(line_levels))
index <- seq_along(line_levels)
K <- outer(index, index, function(i, j) 0.4^abs(i - j))
dimnames(K) <- list(line_levels, line_levels)
L <- t(chol(K))
z1 <- rnorm(length(line_levels))
z2 <- 0.35 * z1 + sqrt(1 - 0.35^2) * rnorm(length(line_levels))
u1 <- setNames(as.vector(L %*% z1) * 0.8, line_levels)
u2 <- setNames(as.vector(L %*% z2) * 0.65, line_levels)
e1 <- rnorm(length(line))
e2 <- -0.2 * e1 + sqrt(1 - 0.2^2) * rnorm(length(line))
dat <- data.frame(
  trait1 = 0.3 + 0.5 * treatment + u1[line] + 0.3 * e1,
  trait2 = -0.2 - 0.25 * treatment + u2[line] + 0.35 * e2,
  treatment = treatment,
  line = line
)

fit_relmat_q2_reml <- drmTMB(
  bf(
    mu1 = trait1 ~ treatment + relmat(1 | p | line, K = K),
    mu2 = trait2 ~ treatment + relmat(1 | p | line, K = K),
    sigma1 = ~ 1,
    sigma2 = ~ 1,
    rho12 = ~ 1
  ),
  data = dat,
  family = biv_gaussian(),
  REML = TRUE
)
corpairs(fit_relmat_q2_reml, level = "relmat")
rho12(fit_relmat_q2_reml)
```

`corpairs()` reports the fitted latent known-matrix correlation among line or
unit-level location deviations. `rho12()` reports the residual correlation
between responses after fixed effects and random effects have been included.
The REML claim stops at `point_fit_recovery`: intervals and coverage remain
unvalidated.

The precision representation remains available for the established ML q2
route, not for this bivariate REML exception:

```r
Q <- solve(K)
fit_relmat_q2_ml <- drmTMB(
  bf(
    mu1 = trait1 ~ treatment + relmat(1 | p | line, Q = Q),
    mu2 = trait2 ~ treatment + relmat(1 | p | line, Q = Q),
    sigma1 = ~ 1,
    sigma2 = ~ 1,
    rho12 = ~ 1
  ),
  data = dat,
  family = biv_gaussian(),
  REML = FALSE
)
```

For reports, keep this intercept-only q=2 known-matrix correlation point-only.
The profile mechanism is callable, but interval calibration and coverage are
still planned for this row:

```r
relmat_pairs <- corpairs(fit_relmat_q2_reml, level = "relmat")
plot_corpairs(relmat_pairs)
```

That plot should be read as a latent known-matrix correlation display, not as a
residual `rho12` plot. Do not add an interval eye until this exact row has
validated interval evidence. The q=4 rows below are also point-only: current q=4
`relmat()` correlations are derived rows.

When the scientific question is whether known-matrix deviations in means and
residual scales covary, use the same labelled `relmat()` term in all four
bivariate endpoint formulas:

```r
fit_relmat_q4 <- drmTMB(
  bf(
    mu1 = trait1 ~ treatment +
      relmat(1 | p | line, Q = Q),
    mu2 = trait2 ~ treatment +
      relmat(1 | p | line, Q = Q),
    sigma1 = ~ treatment +
      relmat(1 | p | line, Q = Q),
    sigma2 = ~ treatment +
      relmat(1 | p | line, Q = Q)
  ),
  data = dat,
  family = biv_gaussian()
)
corpairs(fit_relmat_q4, level = "relmat")
```

The q=4 path is constant across levels and responses: it estimates four latent
known-matrix endpoint scale parameters and six latent correlations. Each
endpoint's level-specific marginal SD also includes the corresponding known
diagonal multiplier. This is not a direct-SD model and it does not make the
residual correlation `rho12` known-matrix structured.

## What to inspect

After fitting, inspect the relatedness layer before interpreting it:

| Output | Use |
| --- | --- |
| `check_drm(fit)` | Confirm the `relmat()` layer was recognized and review matrix/replication diagnostics. |
| `summary(fit)$parameters` | Read the fitted latent known-matrix scale `s`; level `i` has marginal SD `s sqrt(K[i, i])`. |
| `ranef(fit, "relmat_mu")` | Inspect conditional deviations on the fitted relatedness layer. |
| `summary(fit)$covariance` | Check how the `relmat()` SD or q=2 correlation is reported beside other covariance layers. |
| `profile_targets(fit)` | See which relatedness SD or q=2 correlation targets can be profiled directly; q=4 correlations are derived-only in the first slice. |
| `corpairs(fit, level = "relmat")` | Read the fitted q=2 relatedness mean-mean correlation, or the six q=4 relatedness endpoint correlations when the model has the all-four block. |

## Rendered checks

The examples below keep the matrix, standard deviations, and correlation row in
separate visual grammars. The known matrix is raw input structure, not
uncertainty. The latent relmat scale `s` uses the validated default
location-axis bias-corrected small-sample-t Wald channel for the exact K-matrix
q1 row. The plot multiplies its point and interval endpoints by each known
`sqrt(K[i, i])` before comparing node marginal SDs with residual `sigma`, whose
interval uses its ordinary Wald route.
The intercept-only q=2 relatedness correlation is a point estimate with a dotted
zero line. Its profile mechanics are diagnostic-only; interval calibration and
coverage are still planned.

```{r relmat-guide-fit}
relmat_example <- simulate_relmat_guide_data()
relmat_dat <- relmat_example$data
K <- relmat_example$K
Q <- relmat_example$Q

fit_relmat <- drmTMB(
  bf(
    seed_mass ~ temperature + treatment + relmat(1 | line, K = K),
    sigma ~ 1
  ),
  family = gaussian(),
  data = relmat_dat
)
```

```{r relmat-known-matrix, fig.width = 5.3, fig.height = 4.8, fig.cap = "Known relatedness matrix used by the `relmat()` example; this heatmap shows the supplied latent structure, not model-estimated uncertainty.", fig.alt = "Heatmap of the known relatedness matrix for the relmat example. Values are highest on the diagonal and fade as line identifiers are farther apart in the simulated ordering."}
if (requireNamespace("ggplot2", quietly = TRUE)) {
  relmat_matrix <- relatedness_heatmap_data(K)

  ggplot2::ggplot(
    relmat_matrix,
    ggplot2::aes(column, row, fill = relatedness)
  ) +
    ggplot2::geom_tile() +
    ggplot2::coord_equal() +
    ggplot2::scale_fill_gradientn(
      colours = c("#F7FCF5", "#C7E9C0", "#74C476", "#006D2C"),
      name = "Known\nrelatedness"
    ) +
    relmat_guide_theme() +
    ggplot2::theme(
      axis.text = ggplot2::element_blank(),
      axis.ticks = ggplot2::element_blank(),
      panel.grid = ggplot2::element_blank()
    ) +
    ggplot2::labs(
      title = "Known-matrix input structure",
      subtitle = "Validated relatedness among experimental lines",
      x = "Line",
      y = "Line"
    )
}
```

```{r relmat-sd-figure, fig.width = 6.6, fig.height = 3.4, fig.cap = "Fitted residual `sigma` and level-specific known-matrix marginal SDs `s sqrt(K[i,i])` from a univariate Gaussian `relmat()` model. Points and bars are response-scale estimates with 95% Wald intervals transformed by the known diagonal multipliers.", fig.alt = "Horizontal interval display comparing residual sigma with known-matrix node marginal standard deviations calculated by multiplying the fitted latent scale and interval endpoints by the square root of each known covariance diagonal."}
if (requireNamespace("ggplot2", quietly = TRUE)) {
  relmat_ci <- confint(fit_relmat, parm = "variance_components")
  relmat_parameters <- summary(fit_relmat)$parameters
  node_multiplier <- sqrt(diag(K))
  relmat_sd <- relmat_ci[c(1, rep(2, length(node_multiplier))), , drop = FALSE]
  relmat_sd$estimate <- c(
    unname(exp(coef(fit_relmat, "sigma")["(Intercept)"])),
    relmat_parameters[
      relmat_parameters$parm == "sd:mu:relmat(1 | line)", "estimate"
    ] * node_multiplier
  )
  relmat_sd$lower[-1] <- relmat_sd$lower[-1] * node_multiplier
  relmat_sd$upper[-1] <- relmat_sd$upper[-1] * node_multiplier
  relmat_sd$label <- c(
    "Residual\nsigma",
    rep("relmat node\nmarginal SD", length(node_multiplier))
  )

  ggplot2::ggplot(relmat_sd, ggplot2::aes(y = label)) +
    ggplot2::geom_vline(
      xintercept = 0,
      linewidth = 0.45,
      linetype = "dashed",
      colour = "grey55"
    ) +
    ggplot2::geom_errorbar(
      ggplot2::aes(xmin = lower, xmax = upper),
      width = 0,
      linewidth = 1.1,
      colour = "#009E73"
    ) +
    ggplot2::geom_point(
      ggplot2::aes(x = estimate),
      shape = 21,
      size = 3.5,
      stroke = 1,
      fill = "white",
      colour = "#009E73"
    ) +
    ggplot2::scale_x_continuous(
      expand = ggplot2::expansion(mult = c(0.02, 0.05))
    ) +
    relmat_eye_theme() +
    ggplot2::labs(
      title = "Known-matrix marginal SD is separate from residual sigma",
      subtitle = "Node SD = s sqrt(Kii); bars are transformed 95% Wald intervals",
      x = "Fitted standard deviation",
      y = NULL
    )
}
```

For two response means, the known-matrix row is a fitted latent correlation
among line-level deviations. It is not residual `rho12` and it is not a raw
correlation among observations.

```{r relmat-q2-correlation-fit}
relmat_q2_example <- simulate_relmat_q2_guide_data()
relmat_q2_dat <- relmat_q2_example$data
K <- relmat_q2_example$K

fit_relmat_q2_example <- drmTMB(
  bf(
    mu1 = seed_mass ~ age + sex +
      relmat(1 | p | line, K = K),
    mu2 = plant_height ~ age + sex +
      relmat(1 | p | line, K = K),
    sigma1 = ~ 1,
    sigma2 = ~ 1,
    rho12 = ~ 1
  ),
  family = biv_gaussian(),
  data = relmat_q2_dat,
  REML = TRUE
)

relmat_q2_pairs <- corpairs(
  fit_relmat_q2_example,
  level = "relmat"
)
relmat_q2_pairs
```

```{r relmat-q2-confidence-eye, fig.width = 6.4, fig.height = 2.7, fig.cap = "`relmat()` intercept-only q=2 location-location point estimate from `corpairs()`; the dotted vertical line marks zero correlation and no interval is shown because calibration remains planned.", fig.alt = "Single-row point plot for the known-matrix relmat mean-mean correlation, with a hollow point estimate to the right of the dotted zero reference line."}
if (requireNamespace("ggplot2", quietly = TRUE)) {
  relmat_q2_display <- relmat_q2_pairs
  relmat_q2_display$display_label <- "relmat\nmu1-mu2"

  plot_corpairs(
    relmat_q2_display,
    colour = NULL,
    label = "display_label",
    facet = NULL
  ) +
    relmat_eye_theme() +
    ggplot2::labs(
      title = "Known-matrix latent mean correlation",
      subtitle = "Point estimate only; dotted line marks zero",
      x = "Correlation estimate"
    )
}
```

## Boundaries

The following `relmat()` routes remain planned:

- multiple structured slopes such as `relmat(1 + x + z | id, K = K)`;
- known-matrix intercept-slope correlations;
- multiple or labelled residual-scale structured slopes beyond the exact q1
  `relmat(1 + x | id, K = K)` sigma route;
- predictor-dependent `relmat()` `corpair()` regressions;
- generic direct-SD grammar for known-matrix standard deviations;
- bivariate relmat REML with `Q`, slopes, q4 or larger blocks, scale-side
  terms, extra random effects, incomplete pairs, non-unit weights, or
  nonconstant `sigma1`, `sigma2`, or `rho12` formulas;
- non-Gaussian known-matrix structured effects beyond the first
  recovery-grade gates noted below.

Several non-Gaussian families already fit a first q1 unlabelled `relmat()`
structured gate with route-specific evidence. The recovery-grade routes support
point estimates but are not coverage-validated: `Gamma()`, `poisson()`, and
`nbinom2()` each accept a `mu` intercept plus one independent slope;
`nbinom2()` also accepts a `sigma` intercept plus one independent slope.
Separately, `truncated_nbinom2()` (with its hurdle NB2 alias) accepts a
diagnostic-only intercept-only `hu` gate. Use that route only to confirm
fit/extractor feasibility; it does not establish point-estimate recovery.
The `beta()` family and the families with no structured layer still
reject `relmat()`.

Known sampling covariance is a different model layer, not a separate
relatedness route. If the matrix describes known observation-level
uncertainty for effect-size estimates, use `meta_V(V = V)` -- deprecated
`meta_known_V(V = V)` remains a compatibility alias -- and the
[meta-analysis documentation](meta-analysis.html) instead of `relmat()`.

Use the [structural-dependence overview](structural-dependence.html) when you
are choosing among `animal()`, `phylo()`, `spatial()`, and `relmat()`. Use the
detailed [structural-dependence tutorial](https://itchyshin.github.io/drmTMB/articles/phylogenetic-spatial.html) when you
need the current fitted examples, equations, and broader parity ladder.