Commonly used Bayesian diversification methods lead to biologically meaningful differences in branch-specific rates on empirical phylogenies

Abstract

Identifying along which lineages shifts in diversification rates occur is a central goal of comparative phylogenetics; these shifts may coincide with key evolutionary events such as the development of novel morphological characters, the acquisition of adaptive traits, polyploidization or other structural genomic changes, or dispersal to a new habitat and subsequent increase in environmental niche space. However, while multiple methods now exist to estimate diversification rates and identify shifts using phylogenetic topologies, the appropriate use and accuracy of these methods are hotly debated. Here we test whether five Bayesian methods—Bayesian Analysis of Macroevolutionary Mixtures (BAMM), two implementations of the Lineage-Specific Birth–Death–Shift model (LSBDS and PESTO), the approximate Multi-Type Birth–Death model (MTBD; implemented in BEAST2), and the Cladogenetic Diversification Rate Shift model (ClaDS2)—produce comparable results. We apply each of these methods to a set of 65 empirical time-calibrated phylogenies and compare inferences of speciation rate, extinction rate, and net diversification rate. We find that the five methods often infer different speciation, extinction, and net-diversification rates. Consequently, these different estimates may lead to different interpretations of the macroevolutionary dynamics. The different estimates can be attributed to fundamental differences among the compared models. Therefore, the inference of shifts in diversification rates is strongly method dependent. We advise biologists to apply multiple methods to test the robustness of the conclusions or to carefully select the method based on the validity of the underlying model assumptions to their particular empirical system.

Lay Summary

Understanding why some groups of organisms have more species than others is key to understanding the origin of biodiversity. Theory and empirical evidence suggest that multiple distinct historical events, such as the evolution of particular morphological features (e.g., the flower, the tetrapod limb) and competition among species, can produce this pattern of divergent species richness. Identifying when and where on the tree of life shifts in diversification rates occur is important for explaining the origin of modern-day biodiversity and understanding how disparity among species evolves. Several statistical methods have been developed to infer diversification rates and identify these shifts. While these methods each attempt to make inferences about changes in the tempo of diversification, they differ in their underlying statistical models and assumptions. Here we test if these methods draw similar conclusions using a data set of 65 time-calibrated phylogenies from across multicellular life. We find that inferences of where rate shifts occur strongly depend on the chosen method. Therefore, biologists should choose the model whose assumptions they believe to be the most valid and justify their model choice a priori, or consider using several independent methods to test an evolutionary hypothesis.

diversification-rate analyses, phylogeny, BAMM, RevBayes, BEAST, ClaDS2, macroevolution

Introduction

Understanding the patterns and processes that shape the tree of life is one of the central pursuits of biology. However, inferring the tempo of evolution among lineages—the patterns of speciation and extinction that gave rise to our extant biodiversity—remains a difficult problem both theoretically and computationally (Louca & Pennell, 2020; Moore et al., 2016; Rabosky, 2010).

Several methods estimate diversification rates (speciation and extinction rates, individually), assuming that rates are constant across the tree (Morlon, 2014). Recently developed methods have built on constant-rate models by allowing diversification parameters to vary depending on the state of a focal character (Maddison et al., 2007) or, even more recently, among branches of the phylogeny, which allows for lineage-specific diversification rate estimates (e.g., Barido-Sottani et al., 2020; Höhna et al., 2019; Maliet & Morlon, 2022; Rabosky, 2014).

Such lineage-specific methods have the potential to offer powerful insights into our understanding of evolution, such as the potential time dependency of macroevolutionary diversification (Henao Diaz et al., 2019), the macroecological and macroevolutionary causes of the latitudinal diversity gradient (Givnish et al., 2018; Rabosky et al., 2018), and macroevolutionary support of Darwinian and Simpsonian theories of microevolution within adaptive zones (Cooney et al., 2017).

The application of diversification-rate estimation methods, however, has been marred by controversy over their implementation (Meyer et al., 2018; Meyer & Wiens, 2018; Moore et al., 2016; Rabosky et al., 2017; Rabosky, 2018) and by theoretical findings that seemingly undermine the general reliability of inferring diversification parameters from phylogenies of extant species (Helmstetter et al., 2021; Kopperud et al., 2023b, ; Louca & Pennell, 2020). These issues are liable to discourage empiricists, who may wonder if the disagreements among model developers and theorists correspond with biologically relevant inference differences in empirical studies.

To address this question, we assess how inferences under five leading contemporary Bayesian methods—Bayesian Analysis of Macroevolutionary Mixtures (BAMM; Rabosky, 2014); the Lineage-Specific Birth–Death–Shift model (LSBDS; Höhna et al., 2019) and its MCMC-free implementation: Phylogenetic Estimation of Shifts in the Tempo of Origination (PESTO; Kopperud & Höhna, 2023a, ); the approximate Multi-Type Birth–Death model (MTBD; Barido-Sottani et al., 2020); and the Cladogenetic Diversification Rate Shift model (ClaDS2; Maliet et al., 2019)—compare to each other.

While all five methods aim to estimate lineage-specific diversification rates, they differ in how and where rate shifts are allowed to occur.

1. BAMM models diversification rates as varying across lineages by testing among models that include different numbers of diversification-rate regimes (sets of speciation and extinction parameters) and different placements of those regimes in the tree; however BAMM does not model rate shifts on extinct (thus unobserved) branches (Rabosky, 2014).

2. The LSBDS model, as implemented in RevBayes (Höhna et al., 2016), samples rate regimes from a prior distribution discretized into a fixed number of rate categories; this discretization facilitates computation and allows the method to model shifts on extinct branches (Höhna et al., 2019).

3. PESTO is a new implementation of the LSBDS model that analytically computes the posterior mean speciation and extinction rates conditional on a set of hyperparameters without the need for Monte Carlo sampling (Kopperud & Höhna, 2023a, ).

4. The MTBD method is based on a multitype birth–death process that infers the number of rate regimes as well as the transition rate $\gamma$ between rate regimes (Barido-Sottani et al., 2020). This approach allows for the same rate regime to be present in different parts of the tree. The approximate MTBD, tested here, assumes that no rate changes occur in the extinct parts of the tree; this approximation, when applied with a high transition rate prior, has been found to not substantially differ from the exact MTBD method, which allows rates changes along extinct lineages (Barido-Sottani et al., 2020).

5. Finally, in the ClaDS2 model, diversification rates only change at speciation events. Descendant lineages inherit the speciation rate via a stochastic process that is influenced by the $\alpha$ parameter, which represents the long-term trend (i.e., increase or decreases) of the speciation rate (Maliet et al., 2019). This model results in many small and frequent shifts in diversification rate regimes, unlike the other methods, which tend to infer a few large shifts in rate regimes (Maliet et al., 2019; Maliet & Morlon, 2022). Another aspect of ClaDS2 is that extinction rates are not inferred per branch. Instead, the model estimates a global turnover parameter (⁠$ϵ={\mu }_i/{\lambda }_i$⁠). However, shifts are allowed to occur along extinct branches.

Other methods, not tested here, leverage hidden states using a maximum likelihood framework (e.g., Vasconcelos et al., 2022).

To assess whether the theoretical and computational differences among these methods result in biologically meaningful differences, we reanalyze 65 empirical data sets, compiled from Henao Diaz et al., (2019), using each of BAMM, LSBDS, PESTO, MTBD, and ClaDS2. We address the question: do different analytical methods for estimating branch-specific diversification rates produce significantly different results across an array of empirical data sets?

Methods

Empirical data

Our empirical data are derived from the set of 104 chronograms compiled and analyzed with BAMM by Henao Diaz et al., (2019). From the set by Henao Diaz et al., we excluded trees with fewer than 30 extant taxa in order to concentrate on more informative data sets, resulting in our final set of 76 chronograms.

Model settings

Our goal was to apply each method as a typical diligent user might. For each chronogram, we used the incomplete-sampling fraction collected from the original study by Henao Diaz et al., (2019) and applied that sampling fraction when we ran each of the five inference methods. While the methods differ in their specific parameterizations of the birth–death process, we attempted to use comparable settings and priors across methods.

For BAMM analyses, we modified the control files from Henao Diaz et al., (2019). We set lambda to be time constant rather than time variable in order to more closely match the inferences of other methods. We set BAMM priors for each phylogeny using the setBAMMpriors() function in the BAMMtoolsR package (Rabosky et al., 2014b). This function computes data set-specific priors by estimating metrics from the data set such as the root age of the chronogram and then estimating reasonable and broad expectations for shifts and rates. We ran BAMM v2.5.0 using the BAMMtools priors and control files, which determined the phylogeny-specific number of generations for a single MCMC chain. We removed the first 10% of the MCMC samples as burnin and assessed convergence by computing estimates of effective sample size (ESS) using the R package coda (Plummer et al., 2006). We specifically looked for convergence of the log-likelihood parameter and the “number of distinct regimes” parameter, as is recommended (Rabosky et al., 2014a). Analyses that did not reach convergence were run for additional generations until they converged.

For LSBDS analyses, we used the same set of priors for all phylogenies (except for sampling fraction) with eight categories for speciation and for extinction (64 total rate categories). The number of rate categories was chosen after performing a test on one representative phylogeny, which found that increasing the rate categories above 64 did not result in a significant change in model fit. For each chronogram, we ran four MCMC chains for 5,000 generations. Convergence was assessed for each chain by checking that the ESS values for all model parameters in the log files were greater than 200 using the R package coda (Plummer et al., 2006). Chains that did not reach convergence were restarted and run for an additional 5,000 generations. We merged the posteriors, retaining the last 4,000 generations from the MCMC (10% burnin for nonrestarts and 60% burnin forrestarts).

We applied PESTO in a three-step fashion. First, we estimated the parameters of a constant-rate birth–death process and treated these as hyperparameters: the speciation rate (⁠$\lambda$⁠) and the extinction rate (⁠$\mu$⁠). Second, we set up a state-dependent speciation-extinction (SSE) model. In this model, we used rate values that correspond to 10 quantiles of 2 lognormal distributions with medians $\lambda$ and $\mu$⁠, and standard deviation $0.587$⁠. In the SSE model, we used all pairwise combinations of these (i.e., 100 rate categories). Furthermore, we fixed the shift rate parameter (⁠$\eta$⁠) such that we would expect 10 rate shifts across the phylogenies (i.e., $\eta =10/$tree-length). Third, we calculated the posterior state probabilities along each branch. Finally, we plotted the posterior mean rates averaged over the time span for each individual branch.

We ran the MTBD model under default priors (implemented in BEAST2; Bouckaert et al., 2014; Barido-Sottani et al., 2020). We ran three MCMC chains for 100,000,000 generations per phylogeny. We removed the first 25% as burnin and assessed MCMC convergence by checking that ESS values were higher than 200 for all rates.

We ran ClaDS2 using the default priors (as described in Maliet & Morlon, 2022). We ran three MCMC chains for each data set and took a 25% burnin, as is the default setting for ClaDS2. We assessed convergence by calculating the Gelman statistic (Gelman et al., 2014) every 1000^th generation and stopping the analysis once it achieved a Gelman statistic of 1.05, following the standard guidelines for ClaDS2. While we set up our BAMM analyses to be time constant, ClaDS2 does not allow users to fix the inherited speciation-rate parameter (⁠$\alpha$⁠). Thus our analyses estimate $\alpha$ and infer trends of increasing or decreasing speciation rate throughtime.

All methods are conditioned on “survival”: the occurrence of the most recent common ancestor and survival of the two descendent lineages (Nee et al., 1994). However, there are method-specific differences in the details of how methods condition for different rate-shift scenarios that could potentially contribute to different rate estimates across methods.

Processing model output

In cases where MCMC convergence was difficult, we aimed to determine the potential underlying cause. To assess whether the subset of trees where one or more methods failed to converge was substantially different from the subset that did converge, we compared descriptive metrics including phylogeny size, phylogeny age, incomplete sampling fraction, branch length variance, and multidimensional scaling (MDS) via Robinson-Foulds (RF; Robinson & Foulds, 1981) and Kuhner–Felsenstein (KF; Kuhner & Felsenstein, 1994) distances.

Processing model output

We obtained estimates of the relevant diversification parameters (e.g., speciation rate, extinction rate, etc.) from each model. BAMM posterior estimates of speciation rate and extinction rate were extracted using the getMarginalBranchRateMatrix() function in the BAMMtoolsR package (Rabosky et al., 2014b).

We extracted LSBDS posterior distributions from the stochastic branch rate log file produced by the mnStochasticBranchRate() function in RevBayes. In the PESTO analyses, we computed the branch rates averaged across the branch. If ${\lambda }_k$ and ${\mu }_k$ are the rate values in state $k$⁠, and $P_k\left(t\right)$ is the posterior probability of being in state $k$ at time $t$⁠, then the average net-diversification rate along a branch is

$$\frac{1}{t_1-t_0}{\int }_{t_0}^{t_1}\left[\sum _k\left({\lambda }_k-{\mu }_k\right)P_k\left(t\right)\right]𝑑t,$$

(1)

where $t_0$ is the youngest and $t_1$ is the oldest end point of the branch.

The posterior distributions of speciation and extinction rates of the MTBD model were obtained from the extended Newick file produced by BEAST2 using a modified read.beast() function from the treeio package (Wang et al., 2020). As ClaDS2 does not directly infer extinction rates, we calculated extinction rates per branch by multiplying the inferred global turnover value (⁠$ϵ$⁠) by the branch-specific speciation rates (⁠${\mu }_i={\lambda }_iϵ$⁠). For all branches and models, we calculated net diversification by subtracting the extinction rate from speciation rate (⁠${\lambda }_i-{\mu }_i$⁠) per MCMC generation.

Comparing model inferences

To compare inferences among the five models, we (a) visualized rate estimates on individual chronograms, (b) summarized inferences across all chronograms in the data set to reveal systematic differences, (c) identified differences in the location and magnitude of inferred shifts among methods, and (d) tested for overlap in the 95% HPD interval of the posterior distributions (described in Supplementary Section S4).

Visualizing rates on trees

The canonical way of presenting the results of branch-specific diversification-rate analyses is by coloring the branches of the tree by the estimated rates. For each tree, we colored each branch by the posterior median estimate of speciation, extinction, and net diversification to visualize if the methods inferred similar shifts in similar locations on the tree.

Comparing rate estimates by method

To understand whether the methods displayed any consistent differences across the chronograms, we calculated six summary statistics for each tree. For each diversification rate (i.e., speciation rate, extinction rate, and net-diversification rate), we calculated the posterior medians for each branch, and from those posterior medians, we calculated the tree-wide mean and variance in branch rates for each phylogeny. For each of the six summary statistics (mean and variance for each of the three rates), we set up a linear mixed-effect model:

$$\log \left(\text{summary statistic}\right)=X\beta +Zu+r,$$

(2)

with inference method as a fixed-effect categorical predictor (effect sizes $\beta$⁠), phylogeny as a random effect categorical predictor (⁠$u$⁠), and an error term $r$⁠. $X$ and $Z$ are design matrices for the fixed and random effects. We visually checked that the residuals (⁠$r$⁠) were normally distributed and did not suffer from heteroscedasticity; phylogenies that violated these assumptions were excluded from this analysis. For each linear model, we tested if the least-square means of each pair of methods were statistically different using Tukey’s corrected p-value for multiple comparisons.

Location and magnitude of rate shifts

We additionally tested whether the methods inferred consistent locations and magnitudes of rate shifts, using the rooted version of the Kuhner–Felsenstein distance (KF; Kuhner & Felsenstein, 1994). To do this, we first replaced branch lengths of each timetree with the posterior median rate estimate, from a given method, then scaled each branch by the total tree height. This produces a method-dependent tree with branch lengths that provide information regarding the magnitude and location of rate shifts but with identical topology. We calculated KF distances between the rescaled trees from each pair of methods; this distance is equivalent to the mean square error (MSE) given that the two trees being compared have the same topology, as they do in our analyses. For each tree and for each diversification parameter, we computed the mean square error among the different methods:

$$\text{MSE}=\frac{1}{N}\sum _i^N{\left({\lambda }_i-{\lambda }_i^{\prime }\right)}^2,$$

(3)

where ${\lambda }_i$ (or similarly ${\mu }_i$⁠, or $\left({\lambda }_i-{\mu }_i\right)$⁠) is the scaled diversification rate parameter for branch $i$⁠.

A large MSE tells us that the two methods being compared infer different rate magnitudes and/or rate shifts in different locations. A small MSE, however, indicates that the two methods give us similar results.

Diversification rates may vary according to clade age, potentially due to sampling (Henao Diaz et al., 2019; Louca et al., 2022); in particular, younger clades may exhibit faster diversification rates. To pull apart the potentially compounding effects of method and clade age on our results, we also tested whether the method-specific patterns we observe still hold for phylogenies of different ages by splitting our data set into older and younger phylogenies (where older and younger trees are defined as having a root age greater or less than the median root age) and repeating the above-described procedure for both groups.

Computation

We ran all diversification analyses either locally, on the Savio HPC at UC Berkeley, or using the CIPRES Science Gateway v 3.3 (Miller et al., 2010).

We performed all comparison analyses in R v 3.6.0 (R Core Team, 2013). We performed data manipulation with the R packages phytools, (Revell, 2012), tidyverse (Wickham, 2017), reshape2 (Wickham, 2012), readr (Wickham & Hester, 2020), plyr (Wickham, 2011b), and coda (Plummer et al., 2006). We generated plots with R packages see (Lüdecke et al., 2021), ggplot2 (Wickham, 2011a), ggpubr (Kassambara, 2018), ggtree (Yu et al., 2018), ggsignif (Ahlmann-Eltze, 2017), ggExtra (Attali & Baker, 2016), cowplot (Wilke, 2016), RevGadgets (Tribble et al., 2022), and pdftools (Ooms, 2020). We obtained linear mixed models using the R package lmer (Bates et al., 2015) and obtain emeans estimates using the R package eemeans (Lenth, 2020). We additionally used smacof (Mair et al., 2022) and phangorn (Schliep, 2011) to perform MDS and to calculate RF and KF distances. Citations for R packages were generated with RefManageR (McLean, 2014).

Results and discussion

Convergence

Our full data set contains 76 chronograms from multicellular organisms, with 31–4,161 extant tips, root ages of 4.9–349.8 million years ago (MYA), and 0.014%–0.100% of extant species sampled (Supplementary Figure S1A and Supplementary Table S1). All methods converged for 43 trees (the “complete subset”; Supplementary Figure S1B). Of the methods tested, LSBDS had the most difficulty achieving MCMC convergence (it converged for 46 trees). All methods except LSBDS converged for 65 trees (the “partial subset”; Supplementary Figure S1C); PESTO directly computes the posterior mean and thus “convergence” does not apply. Trees that did not converge have poorer taxon sampling (i.e., the ratio of sampled species to total species richness; p-value = .039), older root ages (p-value = .0001), and greater branch-length variance (p-value = .00006) than the converged trees, but sample size (number of tips) was not an important factor (p-value = .076; Supplementary Figure S2C–F). The branch-length variance is consistent with the degree of spread between the KF and RF MDS analyses (Supplementary Figure S2A and B); the KF MDS—which accounts for branch lengths as well as topology—has a larger spread than the RF MDS. Overall, these results fit with our intuitive understanding of the challenges in inferring shifts in diversification rates. We expect that older trees and trees with greater variation in branch lengths should undergo more rate shifts than younger trees and those with less variation in branch lengths. Thus inferring the diversification-rate parameters of these trees should be generally more challenging for the MCMC. These results suggest that users should be particularly attentive to MCMC convergence if their chronogram(s) are poorly sampled, old, or have a lot of branch-length variation, and especially so if they are using LSBDS. In these cases, even more so than usual, it is important to run each MCMC multiple times independently, to assess both stationarity and convergence.

Comparison of methods

Visualizing rates on trees

None of the 43 phylogenies in our “complete subset” had concordant estimates among all methods given our evaluation criteria (Figure 1). For some phylogenies, the methods inferred similar shifts in net diversification (e.g., Figure 1A–E), whereas for others the inferred shifts differed slightly (e.g., Figure 1F–J) or strongly (e.g., Figure 1K–O). We would expect some differences when comparing different modeling approaches, as there are patterns to the differences in our results that can be attributed to the fundamental differences between the models. We illustrate these patterns using a few example phylogenies, which are representative of the patterns one will find when perusing the full set of trees in the supplemental materials (Supplementary Section S5).

Figure 1

Three representative phylogenies with Z-transformed (mean centering and scaling to unit variance) posterior median estimate of net diversification painted on the branches. Columns show estimates from BAMM (A, F, K), LSBDS (B, G, L), PESTO (C, H, M), MTBD (D, I, N), and ClaDS2 (E, J, O). (A–E) Phylogeny of Lindsaeaceae (necklace ferns; Testo & Sundue, 2016), (F–J) Phylogeny of Ruminants (tetrapod; Toljagić et al., 2018), and (K–O) Phylogeny of Odonates (dragonflies and damselflies; Waller & Svensson, 2017). The rate values are in units of events per lineage per million years.

Occasionally, two methods generally identified similar patterns. For example, BAMM and LSBDS identified a similar shift of about the same magnitude in speciation rates for some clades, for example, the Lindsaeaceae (necklace ferns, clade i; Figure 1A and B). Nonetheless, there are still differences between the two methods, for example, a second nested rate shift in the LSBDS and PESTO estimates (clade ii).

In the ruminants (tetrapods) phylogeny (Figure 1F–J), we find that even for results that overall appear similar between methods, there are meaningful differences between their estimates. For example, BAMM, LSBDS, and PESTO inferred a shift around the ancestor of clade i, but LSBDS and PESTO also find approximately two more shifts (Figure 1G, clades ii and iii). Likewise, MTBD differs from the latter two as it infers several shifts in the largest clade and low net-diversification rates on the backbone of that lineage (Figure 1I). Similarly, ClaDS2 infers a slightly different history from all of them, including multiple slowdowns as well as an increase in net diversification within clade i. BAMM, LSBDS, and PESTO identify a shift in approximately the same node (indicated by i), while MTBD infers many replicated increases in rate within clade i. LSBDS and PESTO infer the same shift, which is expected as they are based on the same underlying model andassumptions.

Multiple diversification shifts across a phylogeny are common to many of the MTBD trees (Figure 1I and N; Barido-Sottani et al., 2020). This pattern is caused by the bimodal posterior distribution commonly inferred by this method (Barido-Sottani et al., 2020). Point-estimate summary statistics (e.g., posterior median) of these types of distributions are susceptible to small variation between the ancestor-descendant branches, which causes point estimates to switch between the two optima producing the rapid switching pattern (Figure 1I and N).

Likewise, ClaDS2 is the only model in our analysis that directly models trends in rates through time. The inherited speciation rate (⁠$\alpha$⁠) determines if daughter lineages inherit generally faster or slower rates than the parental lineage. $\alpha$ only changes at cladogenic (speciation) events, which results in many small changes at cladogenetic events, rather than the few large changes that characterize the other methods (Maliet et al., 2019; Maliet & Morlon, 2022). When $\alpha <1$ evolutionary slowdowns occur where the ancestral lineages have higher net-diversification rates than the descendant lineages, a pattern observed in our data (Figure 1O; Moen & Morlon, 2014).

On the other hand, BAMM, LSBDS, and PESTO are similar models and therefore we may expect them to infer similar diversification rates and shifts (Ronquist et al., 2021). While this is sometimes true (e.g., BAMM and LSBDS agree: Figure 1K and L), other times there are pronounced differences (e.g., BAMM is different from both LSBDS and PESTO: Figure 1F and G). Differences between BAMM and the other two models may be due to the well-known differences between the models, namely assuming either rate shifts can (LSBDS and PESTO) or cannot (BAMM) occur on extinct lineages or unsampled lineages (Moore et al., 2016).

Likewise, LSBDS and PESTO are not always in agreement (e.g., Figure 1L and M). While they are mathematically equivalent models, their inferences may differ for two reasons: (a) The overall method of parameter inference and prior distributions. LSBDS infers the mean speciation rate, mean extinction rate, and the shift rate jointly using full Bayesian inference via MCMC using stochastic character mapping. In contrast, PESTO uses empirical Bayes hyperparameters (similar to BAMM), where the mean speciation and extinction rates are first estimated using a constant-rate birth–death process. This means that branch rates estimated using PESTO are conditional on the parameters instead of integrating over the parameter values in the MCMC as in LSBDS. (b) We used a fixed shift rate in PESTO, (⁠$\eta =10/$tree-length) instead of estimating it from the data. While the estimates of the shift rate in LSBDS are sensitive to the shift rate prior, estimates of branch-specific speciation and extinction rates are more robust to their priors (Höhna et al., 2019). We argue that the inference method (empirical Bayes vs. full Bayesian hierarchical model) between PESTO and LSBDS as well as using 64 versus  100 rate categories largely explains the differences in the branch-specific rate estimates.

Comparisons of rate estimates by method

To gain a global perspective of the differences between these models, we calculated two tree-wide summary statistics and distance metrics in order to compare these methods across the entire data set.

We ran all comparisons on the complete subset (the 43 trees that converged for all methods; Supplementary Figure S4) and on the partial subset (the 65 trees that converged for all methods except LSBDS; Figure 2). Comparisons between these two subsets reveal only one small difference (e.g., compare BAMM vs. MTBD—Figure 2C vs. Supplementary Figure S4C), and the most significant differences did not change. Given the large number of data sets that did not converge for LSBDS but converged for all other methods (unconverged data sets = 22) as well as the theoretical similarities between PESTO and LSBDS, we report the following results for the partial data set (but see Figure S4A–F for summaries from the complete data set).

Figure 2

Comparison of tree-wide summary statistics across methods for the partial subset (n = 65). (A–C) Tree-wide mean of posterior median estimates of the branch-specific rate parameters, plotted on a log scale. Asterisks correspond to the p-value of linear mixed model, calculated on the natural log of the rates (*.05 > p-value >.01; **.01 > p-value >.001; ***.001 > p-value). (D–F) Pairwise mean squared error (MSE) between inference methods of phylogenies with branch lengths scaled by rates (speciation, extinction, and net diversification), plotted on a log scale. Split colors correspond to inference method color in (A–C). Distributions closer to zero indicate that the inference methods produced more similar rate estimates, whereas higher values indicate greater dissimilarity. (D) MSE of speciation-scaled phylogenies; (E) MSE of extinction-scaled phylogenies; (F) MSE of net-diversification-scaled phylogenies.

We recover consistent differences in rate estimates among methods, particularly between PESTO and all other models; ClaDS2 also was an outlier, albeit to a lesser extent (Table 1, Supplementary Table S2). In contrast, BAMM and MTBD tended to infer similar speciation and extinction rates. We find that tree-wide average speciation and extinction estimates of PESTO are statistically different from all other methods (Figure 2A and B).

Table 1

Post hoc pairwise comparisons of inference methods using the tree-wide average of summary statistics: speciation, extinction, and net-diversification rates. Columns contain the summary statistics, contrasts of inference methods, the ratios of geometric means, standard errors, degrees of freedom, t-ratios, Tukey-adjusted p-values, significances, and the percent variances explained by the random effect.

Summary statistic	Contrasts	Means ratio	SE	DF	T-ratio	Adj. p-Value	Sig.	% Var.
Speciation	BAMM/ClaDS2	1.0114	0.0336	189	0.3402	.986	N.S.	94.57
	BAMM/PESTO	0.8833	0.0294	189	-3.7326	.001	**
	BAMM/MTBD	1.0352	0.0344	189	1.041	.726	N.S.
	ClaDS2/PESTO	0.8733	0.029	189	-4.0728	.000	***
	ClaDS2/MTBD	1.0236	0.034	189	0.7009	.897	N.S.
	PESTO/MTBD	1.172	0.039	189	4.7736	.000	***
Extinction	BAMM/ClaDS2	1.2616	0.3006	189	0.9752	.764	N.S.	58.65
	BAMM/PESTO	3.0505	0.7268	189	4.6808	.000	***
	BAMM/MTBD	1.119	0.2666	189	0.4718	.965	N.S.
	ClaDS2/PESTO	2.418	0.5761	189	3.7056	.002	**
	ClaDS2/MTBD	0.887	0.2113	189	-0.5034	.958	N.S.
	PESTO/MTBD	0.3668	0.0874	189	-4.209	.000	***
Net diversification	BAMM/ClaDS2	0.6467	0.0406	188.067	-6.9392	.000	***	81.83
	BAMM/PESTO	0.6666	0.0419	188.067	-6.4555	.000	***
	BAMM/MTBD	0.8278	0.052	188.067	-3.0073	.016	*
	ClaDS2/PESTO	1.0309	0.0644	188.0003	.4863	.962	N.S.
	ClaDS2/MTBD	1.2802	0.08	188.0003	3.9523	.001	***
	PESTO/MTBD	1.2419	0.0776	188.0003	3.466	.004	**

While PESTO inferred higher tree-wide average speciation values, the magnitude of the differences is small (ratio of means < 1.2 for all significant contrasts; Table 1). Conversely, PESTO inferred lower tree-wide averages of extinction rates with larger magnitude changes (ratio of geometic means > 1.2; Table 1). The significant difference between PESTO and other methods holds for tree-wide average net-diversification as well, except for the comparison between PESTO and ClaDS2 (Figure 2C), which is not significant.

Additionally, ClaDS2 tree-wide average net-diversification estimates are significantly different from BAMM and MTBD (Figure 2C). A significant difference in net diversification could be driven by the ClaDS2 parameterization of extinction: extinction is not directly estimated in ClaDS2. Therefore, the net diversification rates of ClaDS2 are scaled speciation rates. Alternatively, the differences between methods could be due to the wider variance of net diversification estimates that both BAMM and MTBD have compared to ClaDS2 (Supplementary Figure S3A–C). However, similar to tree-wide average speciation, the magnitude of the difference between the contrast is not large (Table 1). There is also a weakly significant difference between speciation rates of BAMM and MTBD in our partial subset that was not found in the smaller complete subset. Despite these differences, the methods tend to agree on the top three data sets with the fastest net-diversification rate (Supplementary Figure S5), generally inferring that the same set of lineages is evolving remarkably rapidly.

All methods generally had comparable tree-wide average extinction-rate estimates with the exception of PESTO, which sometimes infers much lower extinction rates for some trees than the other methods (though, on average, it infers higher extinction rates). The inference of extinction rate has been the subject of substantial debate, particularly in how failures to account for diversification shifts along extinct branches can impact the likelihood function (Moore et al., 2016; Rabosky et al., 2017). Regardless of the theoretical importance of correctly inferring extinction rates, we demonstrate that differences between extinction and speciation rates manifest in statistically different estimates of net diversification in empirical studies. Therefore, our results indicate that method-dependent tree-wide bias in diversification parameter inference may influence the interpretation of evolutionary shifts in diversification rates.

We find discrepancies between results derived from tree-wide summary statistics and our visual inspection of trees (see section “Visualizing rates on trees”). For example, we find that ClaDS2 and PESTO show no statistical difference in average net diversification (Figure 2C). However, visual inspection of many trees suggests that ClaDS2 and PESTO often differ greatly in the number and position of inferred rate shifts (e.g., Figure 1). Conversely, BAMM and LSBDS often look very similar when we assess individual phylogenies and yet significantly differ when we compare speciation, extinction, and diversification averages (Supplementary Figure S3A–C). This discrepancy reveals the difficulty of summarizing diversification rate estimates across phylogenies to reveal general patterns and motivates the topology-informed rate comparisons, discussed in the following section.

Location and magnitude of rate shifts

We also test whether the models recover similar locations and magnitudes of rate shifts by comparing the mean squared error (MSE) of branch rates; this metric bridges the discrepancies between the global metrics and the observed patterns across the trees (both described above; Figure 2).

When quantifying differences in the location and magnitude of shifts in speciation and net diversification rates, ClaDS2 differs the most (larger MSE), compared with the other methods (Figure 2D and F) and it is by far the biggest outlier across the models when visually inspecting the trees (Figure 1E, J, and O). This result indicates that ClaDS2 estimates differ strongly from those of the other methods in the degree of the shifts it infers, and in their locations. This result is in contrast to the tree-wide averages presented above (see also Figure 2A–C), where ClaDS2 is unexceptional. These results are also corroborated by analyses that take into account uncertainty in rate estimates (see Supplementary Section S4, especially Supplementary Figure S6 and Supplementary Table S3). We recover the same pattern when we subset our data by older and younger trees (Supplementary Figure S7), indicating that our comparisons are robust to clade age.

Tools for assessing methods

Inferred rates generally differ depending on the analysis method; how then should an empirical biologist choose which method to use? A reasonable conclusion from our analyses is that these methods are inherently unreliable. However, given the continued interest (both via use and development) in these methods, we advise empirical users who wish to conduct diversification-rate analyses to take one of two paths.

The first path is to carefully select a method based on the model assumptions. The methods presented in this analysis have theoretical differences in their approach, which appear to produce corresponding differences in results. For example, methods differ in whether shifts in diversification rates are allowed on extinct or unsampled lineages (LSBDS, PESTO, and the “exact MTBD” not tested here), whether diversification rates of each regime are drawn from a continuous distribution (BAMM, MTBD, and ClaDS2) or from a set of discrete rate categories (LSBDS and PESTO), and if shifts occur at cladogenetic events (ClaDS2) or along lineages (BAMM, MTBD, LSBDS, and PESTO). The models make additional assumptions, such as whether shifts in diversification rates affect the process-intrinsic parameters (the speciation and extinction rates) or transformations thereof (e.g., the net-diversification or turnover rate) and whether shifts affect single parameters or combinations of parameters. These assumptions lead to notably different interpretations of how values change through time. Choice of method can be supported by taxon-specific data such as species distribution, fossil record, or phenotypic data (Morlon, 2014). Thus, users should also familiarize themselves with how these models parameterize and estimate diversification rates and ensure that these modeling choices reflect the user’s assumptions about biological processes.

The second path is to critically compare multiple methods when performing diversification analyses. We have shown that—despite the difference in models—in some cases multiple methods produce results with similar biological interpretations. To facilitate the adoption of this practice, we provide R code to easily visualize the results of multiple diversification-rate models across the same phylogeny: https://github.com/Jesusthebotanist/CompDiv_processing_and_plotting.

Regardless of the path taken, these methods are known to be sensitive to prior choice (e.g., LSBDS: prior on number of rate shifts $\eta$⁠, Höhna et al., 2019, MTBD: prior on state change rate $\gamma$⁠, Barido-Sottani et al., 2020, and BAMM: prior on the rate of shifts, Mitchell & Rabosky, 2017; Moore et al., 2016). Method developers generally recommend making inferences under a range of priors to characterize the impact of prior choice on the analysis. While assessing the impact of the prior choice is outside of the scope of this article, we strongly recommend researchers investigate prior sensitivity when using these methods.

The future of diversification analyses

The rise of methods aiming to identify shifts in diversification speaks to the importance of these analyses for understanding the drivers and impacts of important evolutionary events. However, we advocate for caution, for two reasons described below.

First, taking a cautious approach is especially important in light of the many potential problems with these methods, including the controversy surrounding the identifiability of birth–death models (Louca & Pennell, 2020, but see also Helmstetter et al. 2021; Kopperud et al., 2023b ; Legried & Terhorst, 2022; Morlon et al., 2022, among others).

Louca & Pennell (2020) presented a class of birth–death models that are unidentifiable if the rate functions are time-varying (but homogeneous across lineages) and allowed to take any continuous shape. Nonetheless, hypothesis-driven approaches are not allowed to take any arbitrary shape. Since the rate shapes are designed to test diversification scenarios, defined a priori, it has been shown that this approach is not prone to the identifiability issue (Morlon et al., 2022). Even time-varying models that are more agnostic about prior hypotheses are not typically allowed to take any continuous rate shape. Among the “agnostic” models, the piecewise-constant model is the most eminent (Magee et al., 2020; Stadler, 2011), and this model has been proven to be asymptotically identifiable provided there are not too many pieces (Legried & Terhorst, 2022).

However, in spite of the nonidentifiability, inferences of rapidly changing speciation and extinction rates are still typically robust (Kopperud et al., 2023b). The issue of nonidentifiability remains to be investigated thoroughly in lineage-heterogeneous models. These models are more parameter-rich than their homogeneous cousins, and so we do not expect the issue of nonidentifiability to be any simpler here.

Second, we caution against relying too heavily on the estimates from a single method without justifying the assumptions encoded into the model’s choices regarding parameterization and estimation, as we describe in detail in “Tools for assessing methods.”

Some users may wish to use model testing or model adequacy analyses to choose among methods, but we currently do not recommend these approaches. Model testing—especially across models implemented in different software packages—is likely to pick up signals from differences in prior specifications that may mislead or otherwise confuse the results. In addition, Bayes factors for these models, if calculated from marginal likelihoods (e.g., by stepping-stone or path sampling), are unreliable due to the conflation of prior and likelihood in the tree model: such Bayes factors may offer decisive support for the wrong model (May & Rothfels, 2023). Model adequacy approaches have not yet been implemented (to our knowledge) in any of the software for these types of analyses. For example, posterior predictive simulations in PESTO or LSBDS would not condition on the observed rate shifts and instead would simulate from the prior distribution of rates. The development of appropriate methods for model adequacy under these models would be a highly useful tool for future research. However—for now—we recommend that users look toward other sources of data on diversification rates (e.g., the fossil record, expected time to speciation, etc.) to verify whether the models’ estimated rates fall within reasonable expectations given the particularities of the system.

The methods investigated in this article vary in their underlying model and assumptions, but are theoretically related (Ronquist et al., 2021). Due to these model differences, we expect differences in inferences which, in turn, could translate into different biological interpretations. Using a set of empirically derived phylogenies, we show that this is true (Figure 1): no two methods inferred the same shifts for any phylogeny. In some cases, methods generally agreed on the location and timing of inferred shifts, but in other cases, methods strongly disagreed. Method-dependent differences of individual trees were corroborated by tree-wide summary statistics, which indicated small but significant differences between methods (Figure 2; Table 1). While these results hold up when we take into account the uncertainty in rate estimates, we also urge caution in relying too heavily on summary statistics and encourage users to carefully examine their posterior distributions, as 95% HPD intervals vary among methods and distributions may be bimodal, which may mislead common Supplementary statistics (Supplememtary Figure S8, see also Barido-Sottani et al., 2020).

Regardless, it is clear there will be a continued interest in using diversification analyses with a renewed appreciation for the complexities of these methods and the details of how rates are parameterized and estimated.

Data and code availability

All scripts, data, and outputs can be found on Dryad at (doi:10.6078/D18Q68) upon publication. A set of R functions to help users analyze outputs of studied Bayesian methods can be found at https://github.com/Jesusthebotanist/CompDiv_processing_and_plotting.

Author contributions

J.M.-G., M.J.S., and C.M.T. designed and performed the analyses. J.M.-G., M.J.S., C.M.T., S.H., B.T.K., W.A.F, C.D.S., and C.J.R. wrote the manuscript.

Funding

This research used the Savio computational cluster resource provided by the Berkeley Research Computing program at the University of California, Berkeley (supported by the UC Berkeley Chancellor, Vice Chancellor for Research, and Chief Information Officer). This material is based on work supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE-1650441 for J.M.-G. and No. DGE-1752814 for C.M.T. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation. This material is based on work supported by the NSF Postdoctoral Research Fellowships in Biology Program under Grant No. 2209159 to J.M.-G. and No. 2109835 to C.M.T. J.M.-G. was additionally supported by a Cornell Provost Diversity Fellowship. This work was supported by the Deutsche Forschungsgemeinschaft (DFG) Emmy Noether-Program (Award HO 6201/1-1 to S.H.) and by the European Union (ERC, MacDrive, GA 101043187). Views and opinions expressed are, however, those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. The authors declare no conflicts of interest.

Conflict of interest The authors declare no conflict of interest.

Acknowledgments

We thank Dr. Francisco Henao Diaz and Dr. Matthew Pennell for providing the materials necessary to replicate their analyses. We thank Dr. Matthew Pennell and two anonymous reviewers for extensive feedback on the manuscript. We thank Dr. Michael R. May, members of the Rothfels Lab sensu lato, and the Specht lab for comments on the manuscript, we thank the Bonnie Lane Brain Trust for rapid brainstorming and the Cornell University Statistical consulting unit for providing valuable advice on linear mixed models.

Author notes

J.M.-G., M.J.S., and C.M.T. contributed equally to this study.