An Empirical Bayes approach for the study of phenotypic evolution from high-dimensional data

This paper introduces a fast and memory-efficient Empirical Bayes maximum likelihood approach implemented in the R package mvMORPH that enables the analysis of high-dimensional phenotypic traits and complex evolutionary models, such as multi-optima Ornstein-Uhlenbeck processes, overcoming the computational limitations of existing phylogenetic comparative methods.

The Gist

Imagine you are trying to understand how a group of animals evolved over millions of years. You have a family tree (a phylogeny) and you want to study their physical traits.

In the past, scientists could only study a few traits at a time, like the length of a leg or the size of a tooth. But today, technology allows us to scan entire bodies in 3D, capturing thousands of tiny details at once (like the curve of a jawbone or the shape of a skull). This is like trying to describe a person not just by their height and weight, but by every single freckle, wrinkle, and contour on their face.

The problem? The math breaks.

When you have thousands of traits (data points) but only a few hundred species (lineages), the standard mathematical tools used to analyze evolution get stuck. They try to calculate a massive "relationship map" (a covariance matrix) between all those traits, but the map becomes too messy to solve—it's like trying to find a specific grain of sand in a desert while the wind is blowing.

The Solution: The "Empirical Bayes" Shortcut

The authors of this paper, Montoya and colleagues, invented a new mathematical "shortcut" called an Empirical Bayes approach.

Here is how it works, using a simple analogy:

The Old Way (The "Brute Force" Method):
Imagine you are a detective trying to solve a crime with 1,000 suspects and only 50 clues. The old method tries to write down a detailed profile for every single possible combination of suspects and clues to find the perfect match. It requires a supercomputer, takes days to run, and often crashes because there are too many possibilities.

The New Way (The "Smart Guess" Method):
The new method is like a seasoned detective who knows the general rules of the game. Instead of checking every single possibility, they use a "smart guess" (a statistical prior) based on what they already know about how these things usually behave.

• They don't try to map every single grain of sand.
• They assume the sand is mostly uniform, but they allow for a few special spots.
• They use a mathematical trick to "integrate out" the messy parts, focusing only on the signal that matters.

This allows them to solve the puzzle 10 times faster and use 20 times less computer memory than previous methods. It's like switching from a steam engine to a high-speed electric train.

What Did They Discover? (The Mammal Jaw Study)

To prove their new tool works, they applied it to a real-world mystery: How did mammal jaws evolve to fit different diets?

They looked at the 3D shapes of jaws from 95 different mammals (both living and extinct), including carnivores (meat-eaters) and herbivores (plant-eaters).

The Findings:

1. Convergence is Real: Despite being on different branches of the family tree (like marsupials vs. placentals), meat-eaters and plant-eaters evolved remarkably similar jaw shapes.
2. The "Meat" vs. "Plant" Design:
   • Carnivores evolved jaws with a specific shape to handle tearing meat and delivering a strong bite.
   • Herbivores evolved jaws that are deeper and taller, designed to grind tough plant fibers.
3. The "Suckling" Constraint: They also checked if the way babies are raised (marsupials suckle for a long time vs. placental mammals) changed the jaw shape. Surprisingly, diet was the main driver, and the way they were raised as babies didn't override the need to eat specific foods as adults.

Why This Matters

This paper isn't just about math; it's about unlocking the future of biology.

• Speed: Because the new method is so fast, scientists can now analyze massive datasets that were previously impossible to study.
• Complexity: It allows researchers to test much more complex theories about how animals adapt to their environments.
• Accuracy: By looking at all the traits together (rather than just picking a few), they get a much clearer picture of evolution.

In a nutshell: The authors built a faster, smarter engine for evolutionary biology. This engine can now handle the "Big Data" of modern biology, helping us understand how nature shapes life with incredible detail and speed. They proved that even when the data is huge and messy, we can still find the clear patterns of evolution hiding inside.

Technical Summary

1. The Problem

The paper addresses a critical bottleneck in phylogenetic comparative methods (PCMs): the analysis of high-dimensional phenotypic datasets (where the number of traits $p$ exceeds the number of species $n$).

• The Singularity Issue: In high-dimensional settings ($p \geq n$), the trait variance-covariance matrix becomes singular (non-invertible). Since standard multivariate PCMs rely on calculating the likelihood of a multivariate normal distribution, which requires inverting this matrix and computing its determinant, these methods fail or become computationally prohibitive.
• Limitations of Existing Solutions:
  • Dimension Reduction (PCA): Discards covariance information and can lead to inaccurate model selection.
  • Pairwise Composite Likelihood (PCL): Ignores full covariance structures and is not rotation-invariant (problematic for geometric morphometrics).
  • Penalized Likelihood (PL): Uses regularization (shrinkage) to make matrices invertible but relies on computationally intensive cross-validation (CV) to tune the penalty parameter. This makes PL impractical for datasets with thousands of traits or complex models (e.g., Ornstein-Uhlenbeck with multiple optima).

2. Methodology: The Empirical Bayes Approach

The authors propose a new Empirical Bayes (EB) framework to estimate parameters and compare models for high-dimensional data without explicitly computing or storing the full covariance matrix during optimization.

• Core Concept: Instead of maximizing the likelihood of the data given a fixed covariance matrix, the authors integrate the covariance matrix ($\mathbf{R}$) out of the likelihood function using a conjugate prior.
• Mathematical Formulation:
  • They assume the traits follow a matrix-variate normal distribution.
  • They place an Inverse Wishart prior on the trait covariance matrix $\mathbf{R}$.
  • By integrating $\mathbf{R}$ analytically, the marginal likelihood of the data follows a Matrix-Variate T-distribution.
  • This allows the estimation of evolutionary parameters (e.g., selection strength $\alpha$, decay rate $r$) and the regularization parameter ($\mu$) directly via Maximum Likelihood (or Restricted Maximum Likelihood - REML) without iterative cross-validation.
• Regularization Targets: The approach allows for different "shrinkage targets" for the prior scale matrix ($\mathbf{\Psi}$):
  • Scaled Identity ($\mathbf{\Psi}_{iden}$): Rotation-invariant, suitable for geometric morphometric data.
  • Diagonal ($\mathbf{\Psi}_{diag}$): Preserves individual trait variances, suitable for heteroscedastic data.
• Implementation: The method is implemented in the R package mvMORPH within the function mvgls() (method = 'EmpBayes'). It supports Brownian Motion (BM), Early Burst (EB), Pagel's $\lambda$, and Ornstein-Uhlenbeck (OU) models, including the complex OU with multiple optima (OUM).
• Post-Estimation: While the matrix is integrated out for fitting, a regularized estimator of $\mathbf{R}$ can be derived a posteriori as the mean of the posterior Inverse Wishart distribution (Minimum Mean Square Error estimator).

3. Key Contributions

1. Computational Efficiency: The EB approach eliminates the need for cross-validation and explicit matrix inversion/storage during optimization. It is at least 10 times faster and requires up to 50 times less memory than state-of-the-art Penalized Likelihood (PL) methods.
2. Scalability: The method scales effectively to datasets with thousands of traits (demonstrated up to $p=4000$), whereas PL methods often fail or become intractable at these dimensions.
3. Model Complexity: The efficiency gain enables the fitting of complex models previously impossible for high-dimensional data, specifically the Ornstein-Uhlenbeck process with multiple optima (OUM).
4. Robust Model Selection: The computational speed makes feasible the use of computationally expensive model selection tools like the Extended Information Criterion (EIC) and bootstrap-based Likelihood Ratio Tests (LRT), which are more accurate than AIC/BIC in high-dimensional, small-sample scenarios.

4. Results

The authors validated the method through extensive simulations and an empirical application.

• Simulation Performance:

  • Parameter Accuracy: The EB approach accurately estimated evolutionary parameters (e.g., $\alpha$, $\lambda$, $r$) across various models and correlation structures, performing comparably to PL methods and significantly better than naive or pairwise approaches.
  • Covariance Estimation: The regularized covariance estimates showed low Kullback-Leibler (KL) divergence from the true matrices. While PL with quadratic ridge performed slightly better on highly correlated data using Quadratic Loss, EB offered a superior balance of accuracy and computational feasibility.
  • Model Selection: The EIC (Extended Information Criterion) was identified as the most reliable criterion for selecting the correct generating model, outperforming AIC and BIC, particularly in avoiding the misclassification of Brownian Motion as OU.
  • Efficiency: For $n=50$ and $p=4000$, EB was ~10x faster than the fastest PL implementation and ~1000x faster than standard PL with leave-one-out cross-validation. Memory usage was reduced from ~80 GB (PL) to ~3.3 GB (EB).

• Empirical Application (Mammalian Jaw Evolution):

  • Dataset: 3D geometric morphometric data ($p=342$ landmarks) from 95 mammal species (extant and extinct).
  • Hypothesis: Testing adaptive convergence of jaw shape to diet (carnivorous vs. herbivorous) and developmental constraints (marsupials vs. placentals).
  • Findings:
    • The OUM model with two optima (diet-based) was the best-supported model by AIC, EIC, and LRT, indicating strong morphological convergence between carnivores and herbivores across both marsupials and placentals.
    • Models separating optima by developmental mode (marsupial vs. placental) were not supported, suggesting diet is the primary driver of convergence.
    • Morphological Details: Herbivores evolved deeper mandibular corpora and taller rami with specific adaptations for mastication (e.g., masseteric fossa size), while carnivores showed adaptations for bite force and caudal condylar orientation.
    • Selection Strength: The selection strength ($\alpha$) was weak, implying a slow evolutionary response (phylogenetic half-life ~1.1 times the tree age), likely due to the averaging of selection across the entire jaw.

5. Significance

This paper represents a major advancement in evolutionary biology and statistical phylogenetics.

• Unlocking High-Throughput Data: It enables the analysis of massive, high-resolution datasets (e.g., 3D scans, gene expression profiles) that were previously unusable with multivariate PCMs due to computational limits.
• Preserving Covariance: Unlike PCA or distance-based methods, it retains the full covariance structure of traits, which is essential for understanding complex evolutionary constraints and correlations.
• Enabling Complex Hypothesis Testing: By making bootstrapping and complex model fitting (like multi-optima OU) computationally feasible, it allows researchers to rigorously test hypotheses about adaptive radiation, convergence, and developmental constraints in high-dimensional trait spaces.
• Software Accessibility: The implementation in the widely used mvMORPH R package ensures immediate accessibility for the scientific community.