Reduced-rank phylogenetic random effect (PGLLVM)

Adds a reduced-rank phylogenetic random effect on the species dimension. For each of d phylogenetic latent factors, the species-level scores are drawn from a multivariate normal with the user-supplied phylogenetic correlation matrix; the factors are then linearly combined into trait-specific contributions via a lower-triangular n_traits x d loading matrix. Mathematically:

Usage

phylo_rr(species, d = 1)

Arguments

species

An unquoted column name giving the species factor. Levels must match the rownames of phylo_vcv passed to gllvmTMB().

d

Integer; the number of phylogenetic latent factors.

Value

A formula marker; never evaluated.

Details

$$p_{it} = \sum_{k=1}^{d} \Lambda_{\mathrm{phy},tk}\, g_{ik}, \qquad g_{\cdot k} \sim \mathcal{N}(\mathbf{0}, \mathbf{A}_{\mathrm{phy}}).$$

When the user supplies phylo_tree (an ape::phylo object), the implementation builds the sparse inverse \(\mathbf{A}_{\mathrm{phy}}^{-1}\) over tips + internal nodes natively (only ape + Matrix, no MCMCglmm dependency) and evaluates the prior through the quadratic form \(g^\top \mathbf{A}_{\mathrm{phy}}^{-1} g\), exploiting tree-topology sparsity. The construction is the deterministic Hadfield & Nakagawa (2010) algorithm (Journal of Evolutionary Biology 23:494-508), adopted from that method; see Hadfield (2010) Journal of Statistical Software 33(2):1-22 for the MCMCglmm reference implementation used for numerical cross-checks. When the user supplies only phylo_vcv (the dense tip-only covariance matrix), the package falls back to inverting it densely via Matrix::solve() – correct, but does not exploit sparsity.

Used inside a gllvmTMB() formula:

# recommended: pass tree
gllvmTMB(value ~ 0 + trait + phylo_rr(species, d = 2),
         data       = df,
         trait      = "trait",
         unit       = "species",
         cluster    = "species",
         phylo_tree = tree)

# legacy: pass dense covariance only
gllvmTMB(value ~ 0 + trait + phylo_rr(species, d = 2),
         data      = df,
         trait     = "trait",
         unit      = "species",
         cluster   = "species",
         phylo_vcv = Cphy)

Examples

if (requireNamespace("ape", quietly = TRUE)) {
  set.seed(1)
  tree <- ape::rcoal(60); tree$tip.label <- paste0("sp", seq_len(60))
  sim <- simulate_site_trait(
    n_sites = 1, n_species = 60, n_traits = 3,
    mean_species_per_site = 60,
    Cphy = ape::vcv(tree, corr = TRUE),
    sigma2_phy = rep(1, 3), sigma2_eps = 0.2, seed = 1
  )
  df <- sim$data
  df$species <- factor(tree$tip.label[as.integer(df$species)],
                       levels = tree$tip.label)
  Cphy <- ape::vcv(tree, corr = TRUE)
  # \donttest{
  # phylo_rr() is a deprecated alias of phylo_latent() and prints a message.
  fit <- gllvmTMB(value ~ 0 + trait + phylo_rr(species, d = 2),
                  data      = df,
                  trait     = "trait",
                  unit      = "species",
                  cluster   = "species",
                  phylo_vcv = Cphy)
  # }
}
#> ! Formula keyword `phylo_rr()` is a deprecated alias; use `phylo_latent()` for
#>   new code.
#> ℹ `phylo_latent(species, d = K)` =/=> `phylo_rr(species, d = K)`
#> → Aliases will be dropped at the next minor release. See `?phylo_latent`.