Update, 01 August 2016, 2:50PM. This post has been updated to include information contained in the supplemental material of Rabosky et al. 2017, and clarify the difference between branch-specific and tree-wide rate variation.
Back in August, I summarized the main points of a debate over the
reliability of the popular macroevolutionary modeling program BAMM. At the time,
critics Moore et al. (hereafter “MEA”) had published a high profile paper in PNAS arguing that several
crucial aspects of BAMM’s implementation and theoretical underpinnings
hindered its ability to accurately estimate diversification rates. Though the paper stimulated a vigorous online debate, a formal rebuttal from the program’s developers Rabosky et al. (hereafter “REA”) did not appear until this week, published early access in Systematic Biology.
Before getting to the major takeaways of REA’s response, a quick cautionary note: like many scientific discussions that blossom over the internet, the debate over BAMM has involved both peer reviewed and non-peer reviewed dimensions. While blog posts from both MEA and REA are useful to understand different perspectives on the issue at hand, it’s important to understand that these dialogues are on essentially parallel tracks — and so the failure to discuss a particular online critique in a published paper should not be taken to mean that the authors have nothing to say on the topic if the argument has not yet been made in the peer reviewed literature. (REA have themselves stated they will only respond to peer reviewed critiques going forward.) For the purposes of this update, I’ve focused only what has been presented in the PNAS and Systematic Biology papers, but I would encourage users of the program to read both REA’s initial rebuttal on the BAMM website and MEA’s in-depth posts on the Treethinkers blog.
With that out of the way, here’s where things stand.
On likelihood functions and unobserved rate shifts…
Background: To describe locations of diversification rate shifts, BAMM samples locations from a prior distribution, mapping them on to a user-inputted phylogeny. Because this study tree unavoidably excludes extinct or otherwise unobserved lineages, MEA contended the methodology biased extinction probabilities. They thus provided their own likelihood function in turn, allowing for MCMC integration over a term describing unobserved rate shifts.
**Figure 2** from Rabosky et al. 2017. Likelihood calculation in BAMM on a hypothetical phylogeny (A), showing working assumptions that likelihood is conditioned on a subtree leading to the shift (B) and that extinct lineages branching before a rate shift but after the most recent surviving node are also assumed to have a rate shift (C)
Response: REA claim that 1) the alternative likelihood calculation method presented by MEA is flawed and 2) that rate shifts on unobserved lineages do not significantly affect BAMM’s performance when given empirically parameterized values. Their first argument revolves around an analysis of the behavior of the MEA likelihood function transformed to a formal probability. REA claim that because MEA calculate extinction probability as a strict function of speciation and extinction rates on the root segment of the tree, a high chance of extinction near the root can therefore cause the probability of the data to exceed a standard 0 to 1 interval, increasing arbitrarily to infinity. In contrast, the authors note, BAMM conditions extinction probabilities by recursively passing down previous computed values from the tips of the tree to the root.
**Figure 4c** from Rabosky et al. 2017, showing differences between MEA’s paramaterizations for rate shift frequencies used in assessing the influence of extinction on unobserved lineages, and empirically-derived parameters for Cetaceans.
Their second point relies on a reanalysis of data provided along with MEA’s PNAS manuscript. While acknowledging that rate shifts on unobserved lineages are potentially an important issue in development of macroevolutionary models, REA contest that MEA implement an analysis with parameterizations of rate shift frequencies that are much higher than values obtained for their own empirical analyses. After running the same Monte Carlo simulator with adjusted, empirically-parameterized values — rate shifts 100 to 1000 less
frequent than those presented in MEA’s critique — REA found no significant effect of shifts on unobserved lineages affecting extinction probabilities.
On posterior rate shift distributions with extreme prior sensitivity
Background: The prior distribution BAMM uses to sample the number and location of diversification rate shifts is known as a “compound Poisson process (CPP) relaxed molecular clock model,” where the waiting time between rate shifts is exponentially distributed and the locations of rate shifts are uniformly distributed across the tree. MEA claim that 1) this model is statistically incoherent, with a theoretically infinite number of parameterizations in which the data that have equal likelihood, and 2) the results of BAMM show extreme sensitivity to the shape of the prior.
Response: REA dismiss concerns about the statistical incoherency of the CPP model by claiming that while investigating prior sensitivity is important, MEA fail to demonstrate a link between this property and flaws in BAMM’s performance — “demonstrable pathologies in the shape of the posterior.” They emphasize that all Bayesian methods show some degree of prior sensitivity, which is not inherently problematic, as sufficiently informative data should shift the posterior towards its “true” value. Finally, they highlight the use of Bayes factors to assess the evidence in favor of rate variation (and therefore filter out models heavily influenced by liberal prior choice).
On the question of combineExtinctionAtNodes = "random"
Background: Overshadowing much of discussion between REA and MEA is a conflict over differing implementations of a setting to calculate extinction probabilities. Following the initial PNAS critique, REA claimed that MEA had implemented a hidden (or “developer only”) setting (combineExtinctionAtNodes = "random", hereafter “random”), producing pathological results that were otherwise impossible to replicate under default parameters.
Response: REA maintain their stance that the “random” setting was unjustified, and that many of the reported performance issues with BAMM cannot be reproduced under default parameters. However, to test the extent of the influence of this nonstandard implementation on BAMM’s performance, REA repeated all MEA’s analyses with both options for combineExtinctionAtNodes=, and found that only in cases where liberal priors are selected does the “random” option significantly affect posterior distributions. Additionally, REA offer a formal proof and simulation results for why they believe this setting is mathematically unreliable, which they present supplemental material associated with the paper.
On unreliable simulation results
Background: After simulating trees under both constant and variable birth-death processes, MEA claimed that BAMM accurately characterized diversification rate parameters when when diversification rates were constant, but failed to accurately estimate parameters when rates of speciation and extinction vary.
**Figure 8a** from Rabosky et al. 2017, showing the distribution of MEA’s simulated trees by information content, with minimum information content required to detect a rate shift marked.
Response: REA contend that the MEA’s simulated rate trees were characterized by small total size and small tip-to-rate-shift ratios, therefore containing low information content with respect to rate heterogeneity. To support their argument, REA used MEA’s simulation code and identical input parameters, retaining all trees with 2000 or fewer extant tips. While MEA reported a data set with a mean tree size of 89 extant tips, REA found a mean size of 342 tips, a difference they claim indicates significant ascertainment bias toward trees with low information content in MEA’s analysis. To quantify this bias, REA calculated information content (delta log likelihood) for each tree in MEA’s variable rates data set, arguing that a minimum value of 3 is required to detect a rate shift. They found that of 435 rate shifts in 100 phylogenies, only 24 meet this standard. Furthermore, REA argue, MEA’s choice of summary statistic fails to distinguish between low information trees and poor performance: in these cases, branch-specific rate variation estimates will have low accuracy, even as the tree-wide accuracy rate is high. Analyzing only the set of trees with sufficient information content, REA conclude that with a “reasonable” tip-to-shift ratio (10:1 or greater), BAMM’s estimates will narrow in on true rates across all branches.
Conclusions
I am going to go out on a limb and say that we probably haven’t heard the last word on BAMM. But while most of us probably lack sufficient expertise to critically evaluate alternate likelihood functions and Bayesian models, we should all appreciate the high standard of debate MEA and REA have set — a rigorous analysis of the theoretical underpinnings of a widely used tool, complete with publicly accessible data and code. For evolutionary biologists grappling with the difficult task of describing diversification regimes across the tree of life, it’s a welcome discussion.
References
Rabosky, D.L. 2014. Automatic Detection of Key Innovations, Rate Shifts, and Diversity-Dependence on Phylogenetic Trees. PLoS ONE. DOI: 10.1371/journal.pone.0089543
Rabosky, D.L., Mitchell, J.S., Chang, J. 2017. Is BAMM flawed? Theoretical and practical concerns in the analysis of multi-rate diversification models. Systematic Biology. DOI: 10.1093/sysbio/syx037
Moore, B.R., Höhna, S., May, M.R., Rannala, B., Huelsenbeck, J.P. 2016. Critically evaluating the theory and performance of Bayesian analysis of macroevolutionary mixtures. PNAS. DOI: 10.1073/pnas.1518659113