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 (FALSE) { # \dontrun{
library(ape)
tree   <- rcoal(20); tree$tip.label <- paste0("sp", 1:20)
Cphy   <- vcv(tree, corr = TRUE)
fit <- gllvmTMB(value ~ 0 + trait + phylo_rr(species, d = 2),
                data      = df,
                trait     = "trait",
                unit      = "species",
                cluster   = "species",
                phylo_vcv = Cphy)
} # }