Multivariate Phylogenetic Comparative Methods

Introduction

In this tutorial, we will investigate the relationship between two traits in a multivariate setting, specifically the co-evolution between ovipositor length in wasps and style length in the flowers they use to oviposit. We will treat the flower style length as part of the wasp’s niche, and therefore as a trait of the wasp’s biology (which is debatable but useful for this example). This simple two-trait system will help us explore multivariate phylogenetic comparative methods (PCMs).

The mvSLOUCH package will be our primary tool, as it allows for a flexible array of multiple evolutionary models. This package uses the PCMbase optimization algorithm that is further enhanced by PCMbaseRcpp, which significantly increases computation speed. While mvSLOUCH will run without PCMbaseRcpp, it will be considerably slower. Therefore, ensure that all mentioned packages are installed, even if they are not being explicitly called.

Key Concepts in Multivariate PCMs

Before diving into the analysis, it’s important to understand some key concepts:

1. Brownian Motion (BM): A simple model where traits evolve by random drift with a constant rate
2. Ornstein-Uhlenbeck (OU): A model that incorporates both random drift and selection toward an optimal value
3. Selection Matrix (A): In multivariate models, this matrix describes the strength and direction of selection
4. Rate Matrix (Σ): Describes the rate and correlation of random changes in traits

Loading Required Packages

First, let’s load all the required packages for our analysis:

# Load core phylogenetic packages
library(ape)        # For phylogenetic tree manipulation
library(slouch)     # For univariate OU models
library(mvSLOUCH)   # For multivariate OU models

# Load data manipulation packages
library(dplyr)      # For data manipulation
library(plyr)       # Additional data manipulation tools

# Additional packages for visualization
library(ggplot2)    # For enhanced plotting
library(phytools)   # For phylogenetic simulations and tools
library(ape)        # For basic tree manipulation

Loading and Exploring the Phylogenetic Tree

Next, we’ll load our phylogenetic tree data and examine its basic structure:

# Load phylogenetic tree
tree <- read.tree("waspTree.tre")
print(tree)

Phylogenetic tree with 39 tips and 38 internal nodes.

Tip labels:
  Blastophaga_intermedia, Blastophaga_malayana, Eupristina_verticillata, Waterstoniella_brevigena, Platyscapa_corneri, Platyscapa_fischeri, ...

Rooted; includes branch length(s).

The phylogeny contains 39 species with a total of 38 internal nodes. Let’s visualize this tree to better understand its structure:

# Plot the tree with smaller labels for better visibility
plot(tree, cex = 0.5, main = "Wasp Phylogenetic Tree")

Loading and Preparing Trait Data

Now we’ll import our trait data and prepare it for analysis:

# Import trait data and convert to matrix format
Traits <- read.csv("Traits.csv", row.names = 1) %>% as.matrix()
head(Traits)

                           PollOvip  PlantStyle
Blastophaga_intermedia   -2.9957323 -0.91629073
Blastophaga_malayana     -2.1202635 -0.91629073
Eupristina_verticillata  -0.1863296 -0.01005034
Waterstoniella_brevigena  0.1133287  0.44468582
Platyscapa_corneri       -0.3011051  0.12221763
Platyscapa_fischeri      -0.5798185 -0.10536052

The data is stored as a matrix with species names as row names. The first column represents the ovipositor length and the second column represents the style length. Let’s perform a quick summary of these traits:

# Basic statistics for our traits
data.frame(
  Trait = c("Ovipositor_Length", "Style_Length"),
  Mean = apply(Traits, 2, mean),
  SD = apply(Traits, 2, sd),
  Min = apply(Traits, 2, min),
  Max = apply(Traits, 2, max)
)

It’s crucial to verify that all species in our tree are represented in our trait data, and vice versa:

# Check if species names match between tree and trait data
all_match <- all(rownames(Traits) %in% tree$tip.label, 
                tree$tip.label %in% rownames(Traits))
all_match

[1] TRUE

In this case, all data have been pre-processed to match.

If you get a “FALSE” in the above logical operation, you should trim your tree and/or dataset to ensure they match each other. There are multiple ways to do this, with the packages geiger and treeplyr offering helpful functions for this purpose.

Visualizing Trait Data

Let’s visualize the relationship between our two traits:

# Create a more informative scatter plot
ggplot(data.frame(Traits), aes(x = PlantStyle, y = PollOvip)) +
  geom_point() +
  labs(
    x = "Ovipositor Length",
    y = "Style Length",
    title = "Relationship Between Wasp Ovipositor Length and Flower Style Length"
  ) +
  theme_bw()

The plot shows the raw correlation between traits without accounting for phylogenetic relationships. Our subsequent analyses will properly account for the shared evolutionary history.

Fitting Multivariate Evolutionary Models

With our traits and tree properly aligned, we can now fit various multivariate models to our data.

Multivariate Brownian Motion Model

Let’s start with the simplest model: a multivariate Brownian Motion (BM), which assumes traits evolve by random drift with a constant rate:

# Fit multivariate Brownian Motion model
bm_out <- BrownianMotionModel(tree, Traits)

You may notice that the Sxx parameter estimate, which corresponds to the rate matrix of stochastic evolution, is not symmetric. To get the correct result, we need to take the transpose cross-product of that entry:

# Calculate the correct rate matrix by transpose cross-product
tcrossprod(bm_out$ParamsInModel$Sxx)

           PollOvip PlantStyle
PollOvip   3.715832   1.675935
PlantStyle 1.675935   1.455431

This result is also stored in bm_out$ParamSummary. The rate matrix tells us about the variance in evolutionary rates for each trait (diagonal elements) and the covariance between traits (off-diagonal elements).

Multivariate Ornstein-Uhlenbeck Models

Now, let’s fit a more complex model, an Ornstein-Uhlenbeck (OU) model. The OU model incorporates both random drift and selection toward an optimal value. In mvSLOUCH, we can specify the shape of the multivariate rate matrix (Σ) and the selection strength matrix (A).

Let’s start with a full OU model:

# Fit full OU model with diagonal stochastic matrix
ou_out <- ouchModel(tree, Traits, Syytype = "Diagonal")

Let’s examine the key parameters:

# Selection strength matrix
ou_out$FinalFound$ParamSummary$expmtA

                PollOvip  PlantStyle
PollOvip    0.3636448237 0.623861612
PlantStyle -0.0008654139 0.003923432

# Rate matrix
ou_out$FinalFound$ParamSummary$StS

           PollOvip PlantStyle
PollOvip   3.315415   0.000000
PlantStyle 0.000000   2.477881

# Optimal trait values
ou_out$FinalFound$ParamSummary$mPsi.rotated

                reg.1
PollOvip   -0.2772541
PlantStyle -0.3379831

The selection strength matrix (A) indicates how strongly traits are pulled toward their optimal values. The rate matrix (Syy) describes the rate and pattern of random changes. The optimal values represent the trait values that selection favors.

Now, let’s fit additional models that specify different shapes for the A matrix:

# Fit OU model with diagonal A matrix (independent trait selection)
ou_out_d <- ouchModel(tree, Traits, 
                    Atype = "Diagonal",
                    Syytype = "Diagonal")

# Fit OU model with lower triangular A matrix
ou_out_l <- ouchModel(tree, Traits, 
                    Atype = "LowerTri",
                    Syytype = "Diagonal")

# Fit OU model with upper triangular A matrix
ou_out_u <- ouchModel(tree, Traits, 
                    Atype = "UpperTri",
                    Syytype = "Diagonal")

These models represent different assumptions about trait evolution:

• Diagonal A matrix: Assumes that selection acts independently on each trait
• Lower triangular A matrix: Assumes that trait 2 (style length) influences the trait 1 (ovipositor length)
• Upper triangular A matrix: Assumes that evolution of trait 1 (ovipositor length) influences the evolution of trait 2 (style length)

Style optimum affects ovipositor optimum

Multivariate SLOUCH Model

Finally, let’s explore a regression model in which ovipositor length is modeled as a function of style length using the mvSLOUCH approach:

# Fit mvSLOUCH model (regression with OU process)
slouch_out <- mvslouchModel(tree, Traits, 
                          kY = 1, 
                          Atype = "DecomposablePositive",
                          Syytype = "Diagonal",
                          diagA = NULL)

Let’s visualize the relationship between traits with regression lines:

# Enhanced plot of the relationship with regression lines
trait_data <- data.frame(
  OvipositorLength = Traits[, 1],
  StyleLength = Traits[, 2]
)

# Calculate regression lines
opt_intercept <- slouch_out$FinalFound$ParamSummary$mPsi.rotated
opt_slope <- slouch_out$FinalFound$ParamSummary$optimal.regression
evol_slope <- slouch_out$FinalFound$ParamSummary$evolutionary.regression

# Create plot
ggplot(trait_data, aes(x = StyleLength, y = OvipositorLength)) +
  geom_point(size = 3, alpha = 0.7) +
  geom_abline(intercept = opt_intercept, slope = opt_slope, 
              color = "blue", size = 1.2) +
  geom_abline(intercept = opt_intercept, slope = evol_slope, 
              color = "red", linetype = "dashed", size = 1.2) +
  labs(
    x = "Style Length",
    y = "Ovipositor Length",
    title = "Relationship Between Traits with Evolutionary Regressions",
    subtitle = paste("Optimal regression (blue):", round(opt_slope, 3), 
                    "| Evolutionary regression (red):", round(evol_slope, 3))
  ) +
  theme_bw() +
  annotate("text", x = min(trait_data$StyleLength), 
          y = max(trait_data$OvipositorLength), 
          label = "Optimal regression", 
          color = "blue", hjust = 0, vjust = 1) +
  annotate("text", x = min(trait_data$StyleLength), 
          y = max(trait_data$OvipositorLength) - 0.5, 
          label = "Evolutionary regression", 
          color = "red", hjust = 0, vjust = 1)

Notice that the optimal regression (solid blue line) and the evolutionary regression (dashed red line) are very close to each other. This similarity suggests that the evolutionary process has had sufficient time to reach the optimal state, or that the strength of selection is strong relative to the phylogenetic signal.

Adding Phylogenetic Signal to Simulated Data

Let’s explore what happens when we add phylogenetic signal to our data:

# Set seed for reproducibility
set.seed(123)

# Create a copy of the original traits
Traits2 <- Traits

# Add phylogenetic signal to the first trait (ovipositor length)
Traits2[, 1] <- Traits2[, 1] + phytools::fastBM(tree, sig2 = 0.5)

# Fit mvSLOUCH model with data containing phylogenetic signal
slouch_out2 <- mvslouchModel(tree, Traits2, 
                           kY = 1, 
                           Atype = "DecomposablePositive",
                           Syytype = "Diagonal",
                           diagA = NULL)

trait_data2 <- data.frame(
  OvipositorLength = Traits2[, 1],
  StyleLength = Traits2[, 2]
)

opt_intercept <- slouch_out2$FinalFound$ParamSummary$mPsi.rotated
opt_slope <- slouch_out2$FinalFound$ParamSummary$optimal.regression
evol_slope <- slouch_out2$FinalFound$ParamSummary$evolutionary.regression


ggplot(trait_data2, aes(x = StyleLength, y = OvipositorLength)) +
  geom_point(size = 3, alpha = 0.7) +
  geom_abline(intercept = opt_intercept, slope = opt_slope, 
              color = "blue", size = 1.2) +
  geom_abline(intercept = opt_intercept, slope = evol_slope, 
              color = "red", linetype = "dashed", size = 1.2) +
  labs(
    x = "Style Length",
    y = "Ovipositor Length",
    title = "Relationship Between Traits with Evolutionary Regressions",
    subtitle = paste("Optimal regression (blue):", round(opt_slope, 3), 
                    "| Evolutionary regression (red):", round(evol_slope, 3))
  ) +
  theme_bw() +
  annotate("text", x = min(trait_data$StyleLength), 
          y = max(trait_data$OvipositorLength), 
          label = "Optimal regression", 
          color = "blue", hjust = 0, vjust = 1) +
  annotate("text", x = min(trait_data$StyleLength), 
          y = max(trait_data$OvipositorLength) - 0.5, 
          label = "Evolutionary regression", 
          color = "red", hjust = 0, vjust = 1)

The difference between optimal and evolutionary regression is more pronounced in the modified data. This illustrates how phylogenetic signal can affect our interpretation of trait relationships.

Model Comparison

Now, let’s compare all models using the penalized AICc (Akaike Information Criterion corrected for small sample sizes) to determine which model best fits our data:

# Create dataframe for model comparison
model_comparison <- data.frame(
  Model = c("BM", "OU_full", "OU_diag", "OU_lowertri", "OU_uppertri", "slOUCH"),
  AICc = c(bm_out$ParamSummary$aic.c,
          ou_out$FinalFound$ParamSummary$aic.c,
          ou_out_d$FinalFound$ParamSummary$aic.c,
          ou_out_l$FinalFound$ParamSummary$aic.c,
          ou_out_u$FinalFound$ParamSummary$aic.c,
          slouch_out$FinalFound$ParamSummary$aic.c)
)

# Calculate delta AICc and Akaike weights
model_comparison <- model_comparison %>%
  arrange(AICc) %>%
  mutate(
    Delta_AICc = AICc - min(AICc),
    Akaike_Weight = exp(-0.5 * Delta_AICc) / sum(exp(-0.5 * Delta_AICc))
  )

# Format table for better display
model_comparison

The model with the lowest AICc value provides the best fit to our data, taking into account both goodness-of-fit and model complexity. Akaike weights can be interpreted as the probability that a given model is the best among the set of models considered.

Based on the model comparison results, we can interpret the evolutionary process that shaped the relationship between wasp ovipositor length and flower style length. The best-fitting model suggests that these traits evolve primarily through random drift without strong directional selection.

Accounting for Measurement Error

In reality, trait measurements often contain some error. Let’s redo our analyses taking measurement error into consideration. We’ll assume a measurement error of about 10% of the total variation:

# Calculate approximate error values (10% of trait variance)
error_values <- c(
  var(Traits[, 1]) * 0.1,
  var(Traits[, 2]) * 0.1
)

# Fit Brownian Motion model with measurement error
bm_out_e <- BrownianMotionModel(tree, Traits, M.error = error_values)

# Compare BM models with and without error
bm_comparison <- data.frame(
  Model = c("BM without error", "BM with error"),
  AICc = c(bm_out$ParamSummary$aic.c, bm_out_e$ParamSummary$aic.c),
  LogLik = c(bm_out$ParamSummary$LogLik, bm_out_e$ParamSummary$LogLik)
)

bm_comparison

Let’s also fit OU models with measurement error:

# Fit OU model with measurement error
ou_out_e <- ouchModel(tree, Traits, 
                    Syytype = "Diagonal", 
                    M.error = error_values)

Adding measurement error improves our models. This highlights the importance of accounting for measurement error in phylogenetic comparative analyses. Now try adding error to all other models and observe if your conclusions change.

Exercise Questions

1. How would you interpret the biological meaning of the selection strength matrix (A) in the OU models?

2. What does it mean when the optimal regression and evolutionary regression are similar or different?

3. Try fitting models where you switch which trait is dependent and which is independent in the SLOUCH model. How does this change your interpretation?

4. How might you incorporate additional traits or environmental variables into these analyses?

5. What are the limitations of these multivariate PCM approaches?

References

1. Bartoszek, K., Pienaar, J., Mostad, P., Andersson, S., & Hansen, T. F. (2012). A phylogenetic comparative method for studying multivariate adaptation. Journal of Theoretical Biology, 314, 204-215.

2. Hansen, T. F., Pienaar, J., & Orzack, S. H. (2008). A comparative method for studying adaptation to a randomly evolving environment. Evolution, 62(8), 1965-1977.

3. Mitov, V., & Bartoszek, K. (2021). PCMBase: A framework for phylogenetic comparative methods. Methods in Ecology and Evolution, 12(11), 2314-2329.