Skip to contents

Coevolution accessors for q = 4 structured bivariate location-scale fits

Source: R/coevolution-accessors.R
coevolution-accessors.Rd

coevolution_cor(), coevolution_vc(), and coevolution_summary() read the among-axis structure of a q = 4 "coevolution" fit: a bivariate Gaussian location-scale model whose four axes mu1, mu2, sigma1, sigma2 share one structured random effect (phylo(), animal(), relmat(), or spatial()) with a dense 4 x 4 among-axis covariance Sigma_a. They are ported term-for-term from DRM.jl's src/coevo_accessors.jl (#1118) and return the same quantities:

Usage

coevolution_cor(object)

coevolution_vc(object)

coevolution_summary(object)

Arguments

object

A drmTMB fit (native engine) or a drmTMB_julia fit (engine = "julia") of a q = 4 structured bivariate location-scale model.

Value

All three return a plain list with the field names of the DRM.jl originals:

  • coevolution_cor(): cor (the 4 x 4 correlation matrix, symmetric, unit diagonal, dimnames = axes) and axes (the axis labels c("mu1", "mu2", "sigma1", "sigma2"), the row/column order).

  • coevolution_vc(): axes, variance (named numeric, diag(Sigma_a)), sd (named numeric, its square root), and cov (the 4 x 4 Sigma_a).

  • coevolution_summary(): axes, variance, sd, pair (a 6 x 2 character matrix with columns from and to), correlation and covariance (numeric length 6, named from:to, in pair order), cor, and cov.

Details

  • coevolution_cor(): the 4 x 4 among-axis correlation matrix R = D^{-1/2} Sigma_a D^{-1/2}. Its off-diagonals are the coevolutionary correlations – R["mu1", "mu2"] is the among-species correlation of the two trait means, R["sigma1", "sigma2"] the correlation of the two log-scales, and the mean-scale entries the lability couplings.

  • coevolution_vc(): the per-axis variance components diag(Sigma_a), their square roots, and the full covariance.

  • coevolution_summary(): both of the above in a tidy long form – one entry per unordered axis pair (upper triangle, in combn(4, 2) order: mu1:mu2, mu1:sigma1, mu1:sigma2, mu2:sigma1, mu2:sigma2, sigma1:sigma2).

Every value is a deterministic map of the covariance the fit already stores; nothing is re-optimised and no uncertainty is reported (this matches DRM.jl, which reports point estimates here and leaves interval work to its bootstrap). For the residual between-response correlation see rho12() and corpairs(), which also lists these six correlations one row each under level = "phylogenetic" (or the matching structured level).

Which fits qualify

The fit must be family = biv_gaussian() with one structured marker shared by all four axes under a single covariance label, e.g. phylo(1 | p | species, tree = tree) on mu1, mu2, sigma1, and sigma2. Any other fit – a residual-only bivariate model, a univariate model, a q = 2 block on the two means only, an intercept-and-slope q = 4 block on one trait, or a block whose SD is itself modelled with sd_phylo(...) ~ – stores no single 4 x 4 Sigma_a over the four location-scale axes and is refused with an error, as in DRM.jl.

Block-diagonal q = 4 blocks (two labelled 2 x 2 blocks) are accepted: the cross-block correlations are exactly zero by construction and are reported as such.

Scale convention for engine = "julia" fits

DRM.jl reports Sigma_a on the raw branch-length scale of the tree, whereas native drmTMB standardises the phylogenetic covariance to unit height (ape::vcv(tree, corr = TRUE)). For an engine = "julia" fit these accessors rescale the stored covariance by the tree height (variances by height, SDs by sqrt(height)) so that the same model reports the same numbers under both engines – the conversion profile_targets() and confint() already apply to that fit's axis SDs. Correlations are scale-free and are unaffected. On a unit-height tree the two conventions coincide.

Examples

# \donttest{
if (requireNamespace("ape", quietly = TRUE)) {
  set.seed(1)
  n_tip <- 16L
  tree <- ape::compute.brlen(ape::stree(n_tip, type = "balanced"), 1)
  tree$tip.label <- paste0("t", seq_len(n_tip))
  C <- ape::vcv(tree, corr = TRUE)
  Sigma_a <- diag(c(0.8, 0.7, 0.4, 0.4)) %*%
    matrix(c(1, 0.6, 0, 0, 0.6, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1), 4) %*%
    diag(c(0.8, 0.7, 0.4, 0.4))
  A <- t(chol(C)) %*% matrix(rnorm(n_tip * 4), n_tip, 4) %*% chol(Sigma_a)
  rows <- rep(seq_len(n_tip), each = 4)
  x <- rnorm(length(rows))
  dat <- data.frame(
    species = tree$tip.label[rows],
    x = x,
    y1 = rnorm(length(rows), 2 + 0.5 * x + A[rows, 1], exp(-0.7 + A[rows, 3])),
    y2 = rnorm(length(rows), -1 + 0.3 * x + A[rows, 2], exp(-0.5 + A[rows, 4]))
  )
  fit <- drmTMB(
    bf(
      mu1 = y1 ~ x + phylo(1 | p | species, tree = tree),
      mu2 = y2 ~ x + phylo(1 | p | species, tree = tree),
      sigma1 = ~ 1 + phylo(1 | p | species, tree = tree),
      sigma2 = ~ 1 + phylo(1 | p | species, tree = tree),
      rho12 = ~1
    ),
    family = biv_gaussian(),
    data = dat,
    control = drm_control(se = FALSE)
  )
  coevolution_cor(fit)$cor["mu1", "mu2"] # coevolution of the two trait means
  coevolution_vc(fit)$sd
  coevolution_summary(fit)$correlation
}
#> Warning: `drmTMB()`: optimizer reported non-convergence (code 1: singular convergence
#> (7)).
#> ℹ Treat the estimates and standard errors with caution. The optimizer escalates
#>   its preset ladder automatically, so inspect `fit$optimizer_attempts` and run
#>   `check_drm()` to diagnose; consider rescaling the data, within-group
#>   replication, or a penalized/MAP fit.
#>       mu1:mu2    mu1:sigma1    mu1:sigma2    mu2:sigma1    mu2:sigma2 
#>     0.3713559     0.6356354    -0.7922738    -0.4807374    -0.8607459 
#> sigma1:sigma2 
#>    -0.0325600 
# }