Known-matrix relatedness with relmat

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:

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:

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 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:

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:

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	This one-slope model fits independent structured intercept and slope variation. Other structured-slope designs and intercept–slope correlations are not yet supported.
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. For the documented one-slope K/Q route, standard-error intervals use an uncorrected Wald approximation on log standard deviation, so read them cautiously; with eight groups, profile intervals are for diagnosis rather than a calibrated interval estimate. Matching univariate mu and sigma relmat() intercepts estimate one known-matrix mean–scale correlation. Models with several structured scale slopes, labelled scale slopes beyond the documented bivariate models, or correlations between structured slopes are not yet supported.
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. For native REML, we have checked point estimates only in one narrow model: a supplied K, intercept-only location formulas, constant residual formulas, complete pairs, unit weights, and no other model layer. Do not assume interval coverage has been established for other REML models. 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 four-parameter location–scale model. corpairs(level = "relmat") reports six constant latent relatedness correlations; profile intervals for those derived correlations are not yet available.
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	Not yet supported. Wait for a documented guide before using these models.

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:

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

For the tested univariate Gaussian REML models, keep sigma ~ 1, use an unlabelled intercept or independent intercept-plus-one-numeric-slope shape, and set REML = TRUE:

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

Simulation checks 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; its interval coverage has not been tested across simulated datasets. These models 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:

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. For this REML model, we have checked point estimates but not interval coverage.

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

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, treat this intercept-only bivariate known-matrix correlation as a point estimate. The profile calculation can run, but interval calibration and coverage have not yet been established for this model:

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

That plot shows a latent known-matrix correlation, not residual rho12. Its interval has not been established. The all-four-endpoint models below also provide point estimates only: their relmat() correlations are derived values.

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:

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")

This all-four-endpoint model is constant across levels and responses: it estimates four latent known-matrix 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 bivariate mean–mean correlation is reported beside other covariance layers.
profile_targets(fit)	See which relatedness SD or mean–mean correlation targets can be profiled directly; all-four-endpoint correlations are derived point estimates in the current model.
corpairs(fit, level = "relmat")	Read the fitted bivariate mean–mean correlation, or the six 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. For the documented one-response K-matrix model, the latent relatedness scale s uses the default small-sample Wald interval on the location scale. 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 bivariate relatedness correlation is a point estimate with a dotted zero line; its interval calibration and coverage have not yet been established.

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
)

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"
    )
}

Known relatedness matrix used by the relmat() example; this heatmap shows the supplied latent structure, not model-estimated uncertainty.

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
    )
}

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.

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.

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
#>    level group block from_dpar to_dpar   from_coef     to_coef from_response
#> 1 relmat  line     p       mu1     mu2 (Intercept) (Intercept)     seed_mass
#>    to_response     class                                       parameter
#> 1 plant_height mean-mean cor(mu1:(Intercept),mu2:(Intercept) | p | line)
#>    estimate       min       max n_values link_estimate  link_min  link_max
#> 1 0.3286778 0.3286778 0.3286778        1     0.3413456 0.3413456 0.3413456
#>   modelled   conf.status interval_source
#> 1    FALSE not_requested   not_available

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"
    )
}

relmat() intercept-only bivariate mean–mean point estimate from corpairs(); the dotted vertical line marks zero correlation, and no interval is shown because calibration has not yet been established.

Boundaries

The following relmat() routes are not yet supported:

• 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 documented one-slope 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, all-four-endpoint or larger models, 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 models described below.

For several non-Gaussian response distributions, relmat() currently fits one unlabelled structured intercept and, where stated, one independent slope. Simulation studies have checked whether the point estimates recover known values, but have not yet assessed how well the corresponding intervals cover the truth. 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 model. Use that route only to check that a fit and its reported values can be obtained; it does not establish point-estimate recovery. The beta_family() 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 instead of relmat().

Use the structural-dependence overview when you are choosing among animal(), phylo(), spatial(), and relmat(). Use the detailed structural-dependence tutorial when you need the current fitted examples, equations, and broader parity ladder.