Phylogenetic response-column slopes

Adds slope-only predictor deviations across response columns, with a tree relating those columns. If B is the response-column by predictor coefficient matrix, the model is $$\mathrm{Cov}(\mathrm{vec}(B^\mathsf{T})) = K_\mathrm{phy} \otimes \Sigma_\mathrm{predictor}.$$ The RHS must be the response-column factor resolved from trait =, and the tree's tips must match its levels. A single | estimates a full predictor covariance; || gives separate predictor variances with zero covariance. Both bars are identical with one predictor. The term never adds a random intercept; put column intercepts in the main formula with 0 + trait.

Usage

phylo_slope(formula, tree = NULL, vcv = NULL, A = NULL, Ainv = NULL)

Arguments

formula

A slope-only x1 + ... + xP | trait or || trait formula.

tree

An ape::phylo object whose tip labels match the resolved response-column levels. Canonical.

vcv

A labelled response-column covariance/correlation matrix. Legacy alias of A =.

A

A labelled response-column relatedness matrix; alias of vcv =.

Ainv

A labelled response-column precision matrix (inverse of A). Sparse inputs are preserved.

Value

A formula marker; never evaluated.

Details

Scope

This route is covered for Gaussian long-format data with one or more bare, finite numeric predictors. Wide-format column slopes, non-Gaussian responses, transformed/factor bases, latent predictor covariance, and intervals are not supported in this release. Use extract_Sigma(fit, level = "column_slope") for \(\Sigma_\mathrm{predictor}\).

For trait-specific intercept-and-slope covariance, use the corresponding augmented phylo_indep(), phylo_latent(), or phylo_dep() syntax.

A historical non-trait RHS route remains accepted to preserve existing fits. It is compatibility behaviour, not the teaching API; new code puts the resolved response-column factor on the RHS.

See also

slope(), animal_slope(), kernel_slope(), spatial_slope(), phylo_dep(), phylo_indep(), extract_Sigma().

Examples

if (requireNamespace("ape", quietly = TRUE)) {
  set.seed(1); tree <- ape::rcoal(6); tree$tip.label <- paste0("sp", 1:6)
  dat <- expand.grid(site = factor(paste0("s", 1:12)),
                     trait = factor(tree$tip.label, levels = tree$tip.label))
  dat$elevation <- rnorm(nrow(dat)); dat$forest_cover <- rnorm(nrow(dat))
  # Slopes that covary along the tree (K_phy is the tree's correlation).
  L <- t(chol(ape::vcv(tree, corr = TRUE)))
  b_elev <- drop(L %*% rnorm(6, sd = 0.5))
  b_forest <- drop(L %*% rnorm(6, sd = 0.3))
  k <- as.integer(dat$trait)
  dat$value <- b_elev[k] * dat$elevation + b_forest[k] * dat$forest_cover +
    rnorm(nrow(dat), sd = 0.3)
  # \donttest{
  fit <- gllvmTMB(
    value ~ 0 + trait +
      phylo_slope(elevation + forest_cover | trait, tree = tree),
    data = dat, trait = "trait", unit = "site", family = gaussian()
  )
  extract_Sigma(fit, level = "column_slope")
  # }
}
#> $Sigma
#>                elevation forest_cover
#> elevation     0.27684371  -0.09301428
#> forest_cover -0.09301428   0.05665309
#> 
#> $R
#>               elevation forest_cover
#> elevation     1.0000000   -0.7427122
#> forest_cover -0.7427122    1.0000000
#> 
#> $level
#> [1] "column_slope"
#> 
#> $part
#> [1] "dep"
#> 
#> $source
#> $source$type
#> [1] "phylo"
#> 
#> $source$grouping
#> [1] "trait"
#> 
#> $source$labels
#> [1] "sp1" "sp2" "sp3" "sp4" "sp5" "sp6"
#> 
#> 
#> $predictors
#> [1] "elevation"    "forest_cover"
#> 
#> $column_labels
#> [1] "sp1" "sp2" "sp3" "sp4" "sp5" "sp6"
#> 
#> $note
#> [1] "Slope-only response-column coefficients from the phylo source with a full predictor covariance. The response-column source covariance supplies the other Kronecker factor; this is not a trait-level Sigma."
#>