Multi-response phylogenetic mixed models: concepts and application

Abstract

The scale and resolution of trait databases and molecular phylogenies is increasing rapidly. These resources permit many open questions in comparative biology to be addressed with the right statistical tools. Multi-response (MR) phylogenetic mixed models (PMMs) offer great potential for multivariate analyses of trait evolution. While flexible and powerful, these methods are not often employed by researchers in ecology and evolution, reflecting a specialised and technical literature that creates barriers to usage for many biologists. Here we present a practical and accessible guide to MR-PMMs. We begin with a review of single-response (SR) PMMs to introduce key concepts and outline the limitations of this approach for characterising patterns of trait coevolution. We emphasise MR-PMMs as a preferable approach for analyses involving multiple species traits, due to the explicit decomposition of trait covariances. We discuss multilevel models, multivariate models of evolution, and extensions to non-Gaussian response traits. We highlight techniques for causal inference using graphical models, as well as advanced topics including prior specification and latent factor models. Using simulated data and visual examples, we discuss interpretation, prediction, and model validation. We implement many of the techniques discussed in example analyses of plant functional traits to demonstrate the general utility of MR-PMMs in handling complex real-world data sets. Finally, we discuss the emerging synthesis of comparative techniques made possible by MR-PMMs, highlight strengths and weaknesses, and offer practical recommendations to analysts. To complement this material, we provide online tutorials including side-by-side model implementations in two popular R packages, MCMCglmm and brms.

Keywords: evolutionary ecology; generalised linear mixed models; multivariate statistics; phylogenetic comparative methods; trait evolution; variance partitioning.

Fig. 1

Bivariate trait data (y ₁, y ₂) simulated from a basic multi‐response phylogenetic mixed model (MR‐PMM) (see Equation 10) containing different levels of phylogenetic (ρ ^phy) and residual (ρ ^res) correlation (A–D). Simulation conditions for each panel are inset at the top left of each scatterplot, with ellipses providing a visual representation of the strength of correlation within each variance component. Data are plotted in scatterplots, as well as heatmaps arranged against the generating phylogeny. In A, y ₁ and y ₂ have no phylogenetic signal (σ ^phy = 0), but a strong positive residual correlation (ρ ^phy = 0, ρ ^res = 0.9). Clades overlap completely in the scatter plot and the heatmap shows bands of colour across y ₁ and y ₂ that appear random with respect to phylogeny. B shows the opposing situation, y ₁ and y ₂ are positively correlated but entirely with respect to phylogeny (ρ ^phy = 0.9, ρ ^res = 0), with no residual variation in either trait (σ ^res = 0). The scatter plot shows clearly distinguishable clades and a tendency for both within‐ and between‐clade correlation. The extent of between‐clade correlation (the tendency for clades to arrange along a positive slope) depends on the topology of the tree, with deep splits promoting separation of clades along the major axis of co‐variation. The heatmap shows phylogenetic structure weakening across clades from green, to orange to blue as the topology becomes more deeply nested, that is as subclades become less clearly separated in evolutionary time. In C, y ₁ and y ₂ have equal phylogenetic and residual variances, and a strong positive correlation operating on both levels. This scenario, where phylogenetic and residual correlations are similar in sign and magnitude, is likely to be common for many biological traits. In D, both traits contain phylogenetic and residual variance, but correlation is only present on the phylogenetic level. This shows how easily conserved correlations are obscured when residual sources of variation contribute substantially to trait variance. Notably, D represents a set of conditions for which a single‐response PMM (SR‐PMM), such as PGLS, will typically fail to detect a significant association between y ₁ and y ₂ (Westoby et al., 2023).

Fig. 2

Correlation $R=\left(\frac{{\mathrm{\sum }}_{\mathit{ij}}}{\sqrt{{\mathrm{\sum }}_{\mathit{ii}}{\mathrm{\sum }}_{\mathit{jj}}}}\right)=\left({\mathrm{\rho }}_{\mathit{ij}}\right)$ and partial correlation $P=\left(\frac{-{\mathrm{\Omega }}_{\mathit{ij}}}{\sqrt{{\mathrm{\Omega }}_{\mathit{ii}}{\mathrm{\Omega }}_{\mathit{jj}}}}\right)=\left({\mathrm{\rho }}_{\mathit{ij}|\mathit{kl}\dots }\right)$ matrices between species traits x, y and z, where Ω = ∑⁻¹ and kl… indexes all variables other than i and j. The undirected network graphs below provide a qualitative representation of trait relationships, where the absence of an edge between two traits corresponds to a zero off‐diagonal element in the matrix above. To evaluate whether the relationship between x and y can be explained by a covariate z, we compute partial correlations from the precision matrix of trait covariances, which quantifies the linear relationship between x and y when controlling for z. In this example, x and y are strongly correlated, however this is fully explained by their relationships to the covariate z, that is x and y are independent, conditional on z.

Fig. 3

Posterior predictive checks from a fitted multi‐response phylogenetic mixed model (MR‐PMM) of leaf traits across 457 species of Eucalyptus. For each trait, black points represent observed values, while five‐point summaries show the median (dark blue points), 0.5 CI (blue bars), and 0.95 CI (light blue bars) of the posterior predictive distribution for each observation. Predictions are made by conditioning on all random effects, with observations ordered by the predictive mean. δ¹³C, carbon isotope ratio in leaf tissue; N_mass, nitrogen content per dry mass of leaf tissue; LMA, leaf mass per unit area.

Fig. 4

Scatterplots (left) and heat‐maps (right) of three leaf traits across 361 species of Eucalyptus (data filtered to complete cases for plotting). Trait values have been log‐transformed and scaled. For heat maps (right), trait values are aligned with the corresponding species in the phylogeny (centre). For δ¹³C versus N_mass (top left), opposing phylogenetic and non‐phylogenetic correlations reported by the model (Fig. 5) are obscured at the level of species phenotypes. δ¹³C, carbon isotope ratio in leaf tissue; N_mass, nitrogen content per dry mass of leaf tissue; LMA, leaf mass per unit area.

Fig. 5

Results from a multi‐response phylogenetic mixed model (MR‐PMM) of leaf traits across 457 species of Eucalyptus. (A) Variance decomposition reveals large differences between response traits in phylogenetic signal, as well as the relative contributions of different sources of variance. (B) Phylogenetic and non‐phylogenetic between‐species correlation coefficient estimates. Points represent the posterior mean for each estimate, with 50% and 95% confidence intervals (CIs) represented by heavy and light wicks respectively. (C) Correlations and partial correlations between traits represented as network graphs. Edge widths are proportional to the posterior mean of each coefficient estimate, with line type indicating significance at the 95% CI. δ¹³C, carbon isotope ratio in leaf tissue; N_mass, nitrogen content per dry mass of leaf tissue; LMA, leaf mass per unit area.