The relationship between prior sensitivity and extinction probabilities in BAMM

By Brian R. Moore, Sebastian Höhna, Michael R. May, Bruce Rannala, and John P. Huelsenbeck

Synopsis

Our recent study (Moore et al., 2016) explored the statistical behavior of BAMM (Rabosky, 2014). One of the pathologies demonstrated in our paper is the extreme prior sensitivity of BAMM. That is, the estimated number of diversification-rate shifts inferred using BAMM is highly sensitive to the assumed prior on the expected number of diversification-rate shifts. In a recent post, Rabosky attempts to refute this conclusion by claiming that: (1) we used modified and/or undocumented features of BAMM normally “unavailable” to users, and; (2) the prior sensitivity we observe is an artifact of how we forced BAMM to compute extinction probabilities.

In this post, we address these claims by clarifying: (1) the version of the BAMM we used was unmodified, and all the settings we used are available to all BAMM users; (2) the provenance of the options for computing extinction probabilities, and; (3) the relationship between the prior sensitivity of BAMM and its options for computing extinction probabilities. Specifically, the options for computing extinction probabilities were implemented by the BAMM developers (claims that we modified the source code are simply untrue), and these options are freely available to any BAMM user. Critically, the default option for computing estimation probabilities—the “if_different” opion—does not solve the prior sensitivity issue with BAMM, although it uniquely conceals the prior-sensitivity problem. We emphatically agree with the criticism that the option used in our study (the “random” option) for computing extinction probabilities is undocumented. More troubling is the fact that these major changes have not been published or subjected to peer review, and the default setting—which substantially impacts the statistical behavior of BAMM—is clearly invalid (see our explanation below).

1. Introduction

BAMM (Rabosky, 2014) attempts to estimate: (1) the number and location of diversification-rate shifts, and; (2) the diversification-rate parameters (speciation and extinction rates) of every branch in the study tree. Our recent paper (Moore et al., 2016) demonstrates several problems with BAMM: (1) the likelihood function is incorrect (it cannot correctly compute the probability of the data when rates of diversification vary across lineages); (2) the compound Poisson process (CPP) prior model is incoherent (it does not correctly specify the prior distribution for the number and location of diversification-rate shifts across the tree), and; (3) the estimated number of diversification-rate shifts is highly sensitive to the assumed number of diversification-rate shifts a priori (i.e., BAMM exhibits extreme prior sensitivity). Although each of these problems constitutes a serious flaw with BAMM, the prior-sensitivity issue is perhaps the least pernicious. This is because prior sensitivity is relatively straightforward to detect in empirical analyses; that is, simply by repeatedly inferring the posterior number of diversification-rate shifts under a range of prior values on the expected number of diversification-rate shifts. In fact, several careful empirical studies using BAMM have noted that the method exhibits extreme prior sensitivity using this strategy (e.g., Callaghan and McPeek, 2016).

This same strategy was used by Rabosky (2014) and Moore et al., (2016) to explore the prior sensitivity of BAMM under simulation (see post 1). These results are summarized in Figure 1 of our previous post 1; the estimated posterior distribution closely resembles the corresponding prior distribution for the expected number of diversification-rate shifts. Nevertheless, Rabosky (2014) used a procedure for summarizing these results that concealed the prior sensitivity of BAMM (see Figure 2 in our previous post 1). Clearly, if the intention is to assess the prior sensitivity of BAMM, Rabosky’s (2014) procedure for summarizing these results is invalid because it cannot detect prior sensitivity: for any value of the prior, the histogram of posterior modes will always have a value of one, simply because the prior distribution will always have a mode of 1, and the posterior always closely resembles the prior (Figure 1).

Rabosky has claimed that our demonstration of the prior sensitivity of BAMM is “scientifically unjustified” because we did not use the “if_different” option for computing extinction probabilities (introduced in BAMM v.2.5). This is a curious criticism of our results, as it applies equally to the results in his original study (Rabosky, 2014, which similarly did not use the “if_different” option). In any case, let’s explore Rabosky’s various claims to see if they have merit.

2. Our results are based on an unmodified version of BAMM v.2.5 using options available to all BAMM users

Rabosky has claimed that we have used a version of the BAMM source code that implements options not available to non-developers. However, this claim is entirely specious. As noted by careful readers, the version of BAMM we used differs only trivially from BAMM v.2.5 in that it writes extinction probabilities to a file; this modification was necessary to assess the degree of the error made by BAMM in computing the extinction probability, and does not in any way alter how the quantity itself is computed. For easy comparison, our (trivially) modified version of BAMM v.2.5. is available on a GitHub repository provided as supplementary material for Moore et al. (2016).

Further, we note that the various options for the combineExtinctionAtNodes setting (which we describe below) are available to all users of BAMM without modifying any source code—it is simply an option specified by the user in the BAMM input file. That’s how you use the program: if you don’t specify an option for a setting in the input file, the program uses the default option. As we explained in our previous post, we wished to evaluate the performance of BAMM in the most favorable possible circumstances, and so we chose to avoid using obviously incorrect features (i.e., “bugs”), such as the “if_different” option, because doing so would constitute evaluating implementation errors (a “bug report”) rather than identifying critical theoretical and statistical errors inherent to this popular method.

3. The “if_different” option is pathological

3.1. The “if_different” option conceals the prior-sensitivity problem

Recall (post 1) that birth-death processes result in lineages that may go extinction before the present. Consequently, computing the probability of an observed tree requires that we allow for the possibility of unobserved speciation events (i.e., where one of the daughter lineages has gone extinct before the present). The probability that a lineage goes extinct before the present is called the extinction probability, and depends only on the speciation- and extinction-rates under which the lineage arose and its age. BAMM computes extinction probabilities using a recursive algorithm that begins at the tips of the tree and moves down the branches in small time steps. At each time step, it computes the change in the extinction probability using a set of differential equations (see Moore et al., 2016; SI1.2 Likelihood function). When the algorithm reaches an internal node (in versions prior to BAMM v.2.5), the extinction probability for the node is set to the extinction probability of its left descendant lineage. More recent versions of BAMM (> v.2.5) implement five options for computing extinction probabilities at an internal node:

  • if_different: if the extinction probabilities of the left and right descendants branches of an internal node are different (i.e., they exceed some arbitrary numerical tolerance), set the extinction probability of the internal node as the product of the extinction probabilities of the left and right descendant branches; otherwise, set the extinction probability of the internal node to that of either descendant branch (since they are essentially the same).
  • left: set the extinction probability of an internal node to the extinction probability of its left descendant branch (e.g., the hard-coded behavior prior to BAMM v.2.5).
  • right: set the extinction probability of an internal node to the extinction probability of its right descendant branch.
  • random: set the extinction probability of an internal node to the extinction probability of its left or right descendant branch (chosen uniformly at random).
  • favor_shift: set the extinction probability of an internal node to the extinction probability of the descendant branch with the most recent rate shift.

In Moore et al. (2016), we analyzed our simulated data using what seems to be the most sensible option available, “random”. However, we note that all of these options are wrong, and none of them have been subjected to peer review. How do these alternative options for computing extinction probabilities impact the primary symptom of prior sensitivity? We re-analyzed our constant-rate simulations (as described here) using each of these options, and compared the prior to each of the resulting posterior distributions (Figure 1, below).

Figure 1. Prior sensitivity for every combineExtinctionAtNodes option currently implemented in BAMM. We simulated 100 constant-rate birth-death trees similar to the cetacean phylogeny (Steeman et al., 2009). We analyzed each tree under a range of prior values for the expected number of diversification-rate shifts, (columns; gamma = 10, 2, 1, 0.5, 0.1), and using each combineExtinctionAtNodes option available in BAMM (rows; combineExtinctionAtNodes = "if_different", "favor_shift", "random", "left", "right"). Each plot summarizes the prior distribution (grey lines) and the corresponding posterior distribution (blue lines) for the number of diversification processes. With the exception of the "if_different" option, the estimated posterior distribution closely follows the assumed prior distribution over the entire range of parameter space explored. Notably, combineExtinctionAtNodes = "random", "left" and "right" generate effectively identical posterior distributions.

Figure 1. Prior sensitivity for every combineExtinctionAtNodes option currently implemented in BAMM. We simulated 100 constant-rate birth-death trees similar to the cetacean phylogeny (Steeman et al., 2009). We analyzed each tree under a range of prior values for the expected number of diversification-rate shifts, (columns; γ = 10, 2, 1, 0.5, 0.1), and using each combineExtinctionAtNodes option available in BAMM (rows; combineExtinctionAtNodes = “if_different”, “favor_shift”, “random”, “left”, “right”). Each plot summarizes the prior distribution (grey lines) and the corresponding posterior distribution (blue lines) for the number of diversification processes. With the exception of the “if_different” option, the estimated posterior distribution closely follows the assumed prior distribution over the entire range of parameter space explored. Notably, combineExtinctionAtNodes = “random”, “left” and “right” generate effectively identical posterior distributions.

Each of these settings demonstrate some degree for prior sensitivity; however, the “if_different” option is the only option that causes the estimated posterior to depart substantially from the corresponding prior distribution. Notably, the “left”, “right”, and “random” options generate effectively identical results (i.e., they exhibit extreme prior sensitivity). Of all available options, the one that exhibits the least prior sensitivity—the “if_different” option—is the default setting for computing extinction probabilities in BAMM. Of course, one might argue that, if extreme prior sensitivity is a problem, anything that reduces the apparent prior sensitivity is a good thing. However, this argument fails because the “if_different” option is an arbitrary and invalid change to the code that causes errors in the computed extinction probabilities, and these errors cause the estimated posterior to artifactually depart from the assumed prior. In other words, this is a bug that only conceals the symptom of prior sensitivity (i.e., by causing the posterior to depart from the prior).

3.2. The “if_different” option results in nonsensical extinction probabilities

Consider a tree with a single diversification-rate shift from a “blue” rate category to an “orange” rate category (Figure 2, below). For simplicity, we assume that—within each of these rate categories—speciation- and extinction-rates are constant (i.e., they do not depend on time). If we follow a path from a particular tip (in the orange rate category) to the root of the tree, we can examine how the computed extinction probabilities change through time for each of the five combineExtinctionAtNodes options (Figure 3, below; under a constant-rate birth-death process, we can compute the extinction probabilities analytically, which we show as the blue and orange lines in each panel). The correct behavior for this extinction probability (which is not available as an option in BAMM) is that, as soon as there is a rate shift, the extinction probability immediately becomes the extinction probability of the new rate category (Figure 5, left panel). (Rabosky calls this procedure the “recompute” method because it recomputes extinction probabilities for every diversification-rate regime; he further claims that the “recompute” equations are invalid, which we will refute in a forthcoming post.) Now let us consider how the “if_different” option computes the extinction probability.

Figure 2

Figure 2. An example tree with a diversification-rate shift. Imagine a tree that started in the “blue” diversification-rate category, with one rate shift to the “orange” diversification rate category. We can follow the extinction probabilities computed by BAMM along a path (indicated by the dashed line) from a tip in this tree to the root under each of the combineExtinctionAtNodes options to understand how they behave (Figure 3).

Under the “if_different” option, the extinction probabilities at nodes 1 and 2 (Figure 2, above) are unaffected; this is because both descendant lineages belong to the same rate category, and therefore the extinction probabilities will be identical when they are evaluated at each node. However, when we reach node 3, the extinction probabilities from the left-descendant lineage (≈0.21, Figure 3, below) and the right-descendant lineage (≈0.19) are different, so the initial extinction probability for the branch subtending node 3 is set to the product of those numbers (≈0.04). These quantities imply that extinct lineages arising immediately after node 3 have approximately the same extinction probability (≈0.2), while extinct lineages arising immediately before node 3 have a very different (and much smaller) extinction probability (≈0.04), even though lineages arising immediately before and after the node are arising in the same (blue) rate category at essentially the same moment in time!

Extinction probabilities are small numbers (i.e., < 1), so invocation of the “if_different” option will always result in increasingly smaller extinction probabilities (i.e., taking the product of two small numbers gives a much smaller number). Accordingly, when “if_different” option is invoked at a given internal node (i.e., where the extinction probability is computed as the product of the extinction probabilities of its descendant lineages), this will cause a sharp decrease in the extinction probability. This sharp decrease in extinction probability will, in turn, trigger the “if_different” option at all more inclusive nodes along the path the root of the tree. For example, when we reach node 4, the extinction probability from the right lineage (≈0.04) is guaranteed to be different than from the left lineage (≈0.21, since the extinction probability of the right lineage is based on the product taken at node 3), resulting in a “cascade” of products of ever smaller extinction probabilities (second panel in Figure 3, below). The consequence of this procedure is that the extinction probabilities (and thus the likelihood) computed by the “if_different” option depart in a substantial and unpredictable way from the true extinction probabilities! By contrast, the remaining options never take products of extinction probabilities, so that, even they they are still not technically correct, they are at least less incorrect than the “if_different” procedure (Figure 3, remaining panels).

Figure 3

Figure 3. Extinction probability calculations are severely compromised by the combineExtinctionAtNodes = if_different option. Consider tree with a diversification-rate shift depicted in Figure 2, above. For the path indicated by the dashed line in Figure 2, we can compare the extinction probabilities as they are computed by each of the combineExtinctionAtNodes options; for reference, we show the extinction probabilities of a constant-rate birth-death process in either the blue or orange diversification-rate categories (blue and orange lines). We show the correct behavior in the left panel for reference (under a model with no diversification-rate shifts on extinct lineages); here, the extinction probability changes immediately when there is a diversification-rate shift (note: this option is not available in BAMM). The “if_different” option behaves very inappropriately, with the extinction probability drastically decreasing at the node immediately subtending the diversification-rate shift, followed by a “cascade” of decreases at each ancestral node. [NB: The “random” option has multiple possible extinction-probability trajectories; in this case, there are actually four possibilities (two per node), resulting in four overlain dashed lines.]

Discussion

In this post, we have explored the relationship between prior sensitivity and the computation of extinction probabilities in BAMM. Prior sensitivity is generally diagnosed by comparing the estimated posterior distribution to the corresponding prior distribution; i.e., the similarity between the inferred posterior and the assumed prior distributions is the primary diagnostic symptom of prior sensitivity. This diagnostic strategy was used by Rabosky (2014) to assess the prior sensitivity of BAMM; the results of his study demonstrate that estimates of the number of diversification-rate shifts inferred using BAMM are extremely sensitive to the assumed prior on the expected number of diversification-rate shifts. Even a cursory inspection of the prior and posterior distributions makes it clear that BAMM exhibits extreme prior sensitivity (e.g., see Figure 1 from our previous post). However, as we discussed in our previous post, Rabosky (2014) summarized these results in a manner that concealed the extreme prior sensitivity of BAMM (e.g., see Figure 2 from our previous post).

On September 15, 2015, we submitted our initial manuscript to PNAS, which exposed the error in Rabosky’s (2014) summary procedure, and thereby revealed the extreme prior sensitivity of BAMM v.2.3. At our request, Rabosky reviewed the initial submission of our manuscript. On November 4, 2015, Rabosky released a new version of BAMM (v.2.5). This new version implemented four new options for computing extinction probabilities. Specifically, these options provide alternative ways to specify the extinction probability of an internal node based on the extinction probabilities of its two immediate descendant lineages. We have demonstrated in this post that the inferred posterior follows the assumed prior for four of these five options—including the previously hard-coded “left” option—for computing extinction probabilities. However, the “if_different” option uniquely causes the inferred posterior to depart from the assumed prior. Accordingly, the “if_different” option causes the symptom of prior sensitivity to disappear. The “if_different” option was made the default setting in BAMM v.2.5.

The results of our study (like those of his original study, Rabosky 2014) are not based on the new default “if_different” option for computing extinction probabilities. Rabosky has seized on this point in an attempt to dismiss our conclusion that BAMM is prior sensitive by claiming: (1) that we modified the BAMM source code [this statement is untrue]; (2) that we used a “developers-only” (e.g., a “hidden” and “undocumented”) option that is not available to other users [this statement is untrue]; (3) that prior sensitivity is not exhibited by BAMM when using the new, default “if_different” option for computing extinction probabilities. Does this last claim mean that the “if_different” option for computing extinction probabilities solves the underlying prior sensitivity afflicting BAMM?

To help understand the answer to this question, consider the following analogy. Imagine that your vehicle has a serious engine problem (e.g., it is out of engine oil because of a serious leak). You will be alerted to this fact by the “check-engine” light in your instrument panel. The check-engine light is the symptom of a serious underlying problem. A really, really bad mechanic might suggest “solving” this mechanical problem by cutting the power to your instrument panel. After all, this constitutes a “solution” to the extent that the check-engine light is no longer illuminated. But this clearly does not solve the underlying problem (the engine is still out of oil). Moreover, this “solution” is likely to cause other serious problems that will adversely impact the operation of the vehicle (e.g., you will not be able to monitor the vehicle speed, fuel level, engine temperature, etc.).

As we have demonstrated above (section 3.2), all of the options for computing extinction probabilities at internal nodes in BAMM are incorrect. However, the error in the computed extinction probabilities is most severe using the default “if_different” option (Figure 3). The errors in the extinction probabilities computed under the “if_different” option will, in turn, cause errors in the computed likelihoods. The errors in the computed likelihoods will, in turn, cause errors in the computed posterior probabilities. The errors in the computed posterior probabilities will, in turn, cause the estimated posterior probability distribution to depart from the assumed prior distribution. The “if_different” option, therefore, causes the diagnostic symptom of prior sensitivity—the similarity of the estimated posterior to the assumed prior—to disappear. However, like cutting power to the warning light, this “solution” only masks the underlying problem; i.e., the “if_different” option only conceals the prior sensitivity in BAMM, it does not solve the underlying problem. As in our automotive analogy, the “if_different” option also causes other adverse impacts on the behavior of BAMM; the increased error in the computed likelihoods will adversely impact parameter estimates, etc.

Accordingly, Rabosky’s claim that BAMM does not exhibit prior sensitivity using the default “if_different” option is technically true, but only to the extent that it conceals the symptom of prior sensitivity. In fact, the “if_different” option is a bug that causes extreme errors in the extinction probabilities (and the likelihoods and parameter estimates). The only “positive” aspect of this bug is that it conceals a pathology of BAMM. However, if we are interested in assessing the statistical behavior of BAMM, then the fact that it conceals problems with the method is clearly undesirable. For this reason, we explored the statistical behavior of BAMM (including the prior sensitivity) using the most correct option for computing extinction probabilities (the “random” option). It would be unfair and confusing to use the default “if_different” option because: (1) we know this is a serious bug, and; (2) this bug is likely to cause pathologies that complicate our attempts to characterize inherent problems with the method.

Finally, we concede that Rabosky has made one valid criticism: the “random” option for computing extinction probabilities is “hidden” and has not been “documented”. That is, this option for computing extinction probabilities in BAMM (like all of the other options for computing extinction probabilities) has not been published, has not been subjected to peer review, and is clearly incorrect. In fact, this is the case for much of the current code in BAMM: the method currently implemented in BAMM is not the method that is described in Rabosky (2014) or in any peer-reviewed paper. We encourage the BAMM developers to expose the major changes to the theory implemented in BAMM to a formal, external review process. More careful and transparent method development would hopefully help identify and avoid the most obvious errors.

Data availability

All of the phylogenetic data and BAMM config files used for this blog post are publicly available as a BitBucket repository.

References

Callahan, MS, and McPeek, MA. 2016. Multi-locus phylogeny and divergence time estimates of Enallagma damselflies (Odonata: Coenagrionidae). Molecular Phylogenetics and Evolution 94: 182–195.

Moore, BR, S Höhna, MR May, B Rannala, and JP Huelsenbeck. 2016. Critically evaluating the theory and performance of Bayesian analysis of macroevolutionary mixtures. Proceedings of the National Academy of Sciences, USA, (in press).

Rabosky DL. 2014. Automatic detection of key innovations, rate shifts, and diversity-dependence on phylogenetic trees. PLoS One, 9:e89543.

Rabosky DL. 2016. Is the posterior on the number of shifts overly sensitive to the prior? http://bamm-project.org/prior.html

Steeman, ME, MB Hebsgaard, RE Fordyce, SY Ho, DL Rabosky, R Nielsen, C Rahbek, H Glenner, MV Sørensen, and E Willerslev. 2009. Radiation of extant cetaceans driven by restructuring of the oceans. Systematic Biology, 58:573–585.