# BAMM provides unreliable parameter estimates

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

## Synopsis

In our recent study (Moore et al., 2016), we explored the ability of BAMM to provide reliable estimates of branch-specific diversification-rate parameters (speciation and extinction rates). To make our results comparable to those presented in the original study (Rabosky, 2014), we summarized the results of our simulation using the identical procedure as that used in the original study. The results of our study—when summarized in a comparable and appropriate manner—reveal that variation in branch-specific estimates of diversification parameters (speciation and extinction rates) are uncorrelated with the true (simulated) values. In a recent post on the BAMM website, Rabosky attempts to dismiss this conclusion by claiming that: (1) he is unable to replicate our results; (2) our demonstration that BAMM cannot reliably estimate diversification-rate parameters is an artifact of our use of the “random” option for computing extinction probabilities, and; (3) finally, he claims that BAMM provides reliable parameter estimates when using the default “if_different” option for computing extinction probabilities. In this post, we discuss each of these claims, demonstrating that: (1) the results regarding parameter estimates presented in Moore et al. (2016) are summarized in the same way as those in the original study (Rabosky, 2014); (2) parameter estimates using BAMM are biased regardless of whether we use the random or if_different options for computing extinction probabilities, and; (3) Rabosky’s claim that BAMM provides reliable parameter estimates is an artifact of his adoption of a new (and entirely invalid) way of summarizing these results that happens to conceal the pathological behavior of BAMM.

## 1. We summarized our results using the same procedure as Rabosky (2014)

BAMM 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. Although we have previously demonstrated that BAMM cannot estimate the number and location of diversification-rate shifts (see post 1 and post 2), we might hope that BAMM may still be able to estimate the branch-specific diversification-rate parameters. If this were the case, it would be possible to informally infer the pattern of diversification-rate shifts by inspecting how the posterior estimates of diversification rates vary across branches. In our paper, we demonstrated that BAMM cannot reliably estimate diversification-rate parameters; in fact, we demonstrated under simulation that estimated variation in branch-specific speciation and extinction rates is generally uncorrelated with the true parameter values when we summarize the results following the procedure advocated by Rabosky (2014). Here, we show that our results are summarized identically to Rabosky (2014). Rabosky has responded to our paper claiming that he is unable to replicate the results of our study. We are at a loss to explain this claim, as we simply followed the summary procedure that he used in his study (Rabosky, 2014).

### 1.1. How did Rabosky (2014) summarize diversification-rate parameter estimates?

To assess the accuracy of parameter estimates obtained using BAMM, Rabosky (2014) simulated 500 trees under each of five alternative scenarios where diversification rates varied across branches. He then estimated branch-specific speciation- and extinction-rate parameters using BAMM, and computed the mean of the posterior distribution of the speciation and extinction rates for each branch of each tree. Next, for each tree, he computed the “correlation” between the true rates and estimated rates by fitting a linear regression, using the true rates as the predictor variable and the estimated rates as the response variable. What do the slopes of the correlation coefficients indicate? A positive regression slope indicates that the relative parameter estimates are reliable; the method is correctly inferring the difference between branches with relatively high rates and branches with relatively low rates in the corresponding simulated tree. Conversely, a negative regression slope indicates that the relative parameter estimates are negatively correlated; the method is incorrectly inferring that relatively high-rate branches have relatively low rates (and vice versa). Finally, a regression slope of zero indicates that the relative parameter estimates are uncorrelated; the method fails to correctly distinguish between branches with relatively high rates and branches with relatively low rates. In other words, positive slopes indicate that BAMM provides reliable parameter estimates (the accuracy of the relative-rate estimates increase as the slopes approach one), whereas negative and zero slopes indicate that BAMM provides unreliable parameter estimates. For each of the five rate-variable simulation scenarios, Rabosky (2014) generated 500 slopes (1 per tree).

Rabosky also assessed the performance of parameter estimates by computing the proportional error of the branch-specific parameter estimates, defined as:

where refers to a particular parameter (e.g., speciation or extinction rate) and refers to a particular branch; a proportional error near 1 indicates that the estimated value is close to the true value. Rabosky then computed the mean proportional error for each tree.

Finally, Rabosky generated histograms of the regression slopes (Figure 1, left column) and mean proportional errors (Figure 1, right column) for each of the five variable-rate scenarios (one for each of the rows in Figure 1). The summaries presented by Rabosky (2014) are tightly peaked around 1, which he argues indicates the ability of BAMM to provide accurate and unbiased estimates of diversification parameters.

Figure 1. Precision and bias of branch-specific speciation-rate estimates using BAMM (this is Figure 6 reproduced from Rabosky, 2014). Rabosky (2014) simulated 500 trees for each of five alternative scenarios where rates vary across branches (rows). For each scenario, he estimated the relationship between the estimated branch-specific rates and their corresponding true values by fitting a linear model to each simulated tree, using the true branch-specific speciation rate as predictor variable and the estimated branch-specific speciation rate as the response variable, and then summarized the set of 500 slopes for each simulation scenario as a histogram (left column). He also estimated the mean proportional error in the branch-specific speciation-rate estimates for each tree, and generated histograms of the proportional error across all branches of all 500 trees (right column). Both of these histograms have a mean near 1, suggesting that the parameter estimates are relatively accurate and unbiased.

### 1.2. How did Moore et al. (2016) summarize diversification-rate parameter estimates?

We followed the above procedure to make our results comparable to those presented in the original study (Rabosky, 2014), with one small exception. Rather than computing the mean proportional error of all the branches per tree, we simply computed (and plotted) the proportional error across all branches (i.e., without taking the mean of the error per tree); this decision was motivated by the concern that taking the mean of the errors across a tree could mask systematic errors within individual trees, but it turns out to make very little difference (see Figure 3 below). Our results demonstrate that branch-specific variation in diversification parameters (speciation and extinction rates) estimated by BAMM is uncorrelated with the true (simulated) values (see Figure 2, which is Figure 6 from Moore et al., 2016) when data are simulated under the BAMM model. To increase the ease of comparison, we have also summarized our results as histograms, following Rabosky (2014) (see Figure 3, below).

Figure 2. Precision and bias of branch-specific speciation-rate estimates (top row) and extinction-rate estimates (bottom row) using BAMM (this is Figure 6 reproduced from Moore et al. 2016). In our study, (Moore et al., 2016), we simulated 100 trees under the BAMM model. We then estimated the relationship between the estimated branch-specific rates and their corresponding true values by fitting a linear model to each tree, using the true branch-specific rates as predictor variables and the estimated branch-specific rates as the response variables, and plotted these linear models across the 100 trees (left column). We also estimated the proportional error in the branch-specific estimates, and plotted the proportional error as a function of the true parameter value across all branches of all 100 trees (right column).

Figure 3. Precision and bias of branch-specific speciation-rate estimates (top row) and extinction-rate estimates (bottom row) using BAMM. In our study, (Moore et al., 2016), we simulated 100 trees under the BAMM model. We then estimated the relationship between the estimated branch-specific rates and their corresponding true values by fitting a linear model to each tree, using the true branch-specific rates as predictor variables and the estimated branch-specific rates as the response variables, and generated histograms of these slopes across the 100 trees (left column). We estimated the proportional error in the branch-specific rate estimates, and generated histograms of the proportional error across all branches of all 100 trees (middle column). We also estimated the mean proportional error in the branch-specific rate estimates for each tree, and generated histograms of the proportional error across all branches of all 100 trees (right column). Our results demonstrate that, when data are simulated under the birth-death model used by BAMM, variation in diversification-rate estimates is, on average, uncorrelated with the true variation in rates, and, on average, diversification-rate estimates are biased upward by ~30% for the speciation rate and ~5000% for the extinction rate.

## 2. The inability of BAMM to estimate diversification-rate parameters is not an artifact of our analysis options

Rabosky also attempts to dismiss the conclusion that BAMM does not provide reliable parameter estimates by claiming that our conclusions are an artifact of using the (more sensible) “random” option for computing extinction probabilities, and that BAMM provides reliable estimates of diversification-rate parameters when using the default “if_different” option for computing extinction probabilities. Fortunately, it is straightforward to evaluate this claim: to this end, we simply reanalyzed all of the rate-variable simulated trees using the “if_different” option for computing extinction probabilities, and summarized the resulting parameter estimates following the protocol used in Rabosky (2014) and replicated in Moore et al., 2016. These results reveal Rabosky’s claim to be patently false: our conclusions are demonstrably not an artifact of using the (more sensible) “random” option for computing extinction probabilities. Rather, these results confirm that variation in diversification-rate parameters estimated using BAMM is uncorrelated with the true variation in rates regardless of the option used to compute extinction probabilities (Figure 5). The estimated error in diversification rates is virtually identical using either the “random” or “if_different” option.

Figure 4. Precision and bias of branch-specific speciation-rate estimates (top row) and extinction-rate estimates (bottom row) using BAMM with the “combineExtinctionAtNodes = if_different” option for computing extinction probabilities. We repeated the analysis performed in Moore et al. (2016), this time using the “combineExtinctionAtNodes = if_different” option, and generated summaries as in Figure 2 above. Results for both options are nearly identical, indicating that the results presented in Moore et al. (2016) are not an artifact of the “combineExtinctionAtNodes = random” option for computing extinction probabilities.

Figure 5. Precision and bias of branch-specific speciation-rate estimates (top row) and extinction-rate estimates (bottom row) using BAMM with the either option for computing extinction probabilities. We repeated the analysis performed in Moore et al. (2016), this time using the “combineExtinctionAtNodes = if_different” option for computing extinction probabilities (red histograms), and replicated summaries using the “random” (depicted from Figure 3, above; grey histograms) to facilitate comparison. As with Figure 4, results for both options are virtually identical, indicating that the results presented in Moore et al. (2016) are not artifacts of the “combineExtinctionAtNodes = random” option for computing extinction probabilities.

## 3. Lies, damned lies, and summary statistics

### 3.1. Pooling estimates of branch-specific diversification rates across trees is misleading

BAMM is intended to estimate diversification-rate parameters (speciation and extinction rates) of every branch in the study tree. Therefore, evaluating the ability of BAMM to accurately estimate rate parameters must focus on individual trees. Rabosky (2014) recognized this, and accordingly summarized the relationship between the true branch-specific rates and estimated branch-specific rates for each simulated tree. As we discussed above, we simply followed this procedure in our study (Moore et al., 2016). However, in attempting to invalidate the results presented in our PNAS paper, Rabosky chose to assess the relationship between true and estimated diversification rate parameters in a new—and misleading—way. Specifically, rather than assessing correlations between true and estimated rates for individual trees, Rabosky instead infers correlations by pooling the true and estimated branch-specific rates across all trees in the simulation. This new procedure has the effect of conflating between-tree rate variation and within-tree rate variation, thus confounding our ability to assess whether BAMM is correctly estimating rate variation within individual trees.

Consider two trees with among-lineage speciation-rate variation, such that some lineages within each tree are speciating twice as fast as the other lineages in the tree. Further, imagine that the “average” speciation rate of the first tree is twice that of the second tree. Now imagine we use BAMM to estimate branch-specific speciation rates for each of those trees, but that it is unable to accurately characterize how speciation rates vary across the branches of each tree (the quantity we actually care about!), but BAMM can nonetheless estimate the “average” speciation rate of each tree. Estimating the correlation between the true and estimated speciation rates for each tree separately (correctly) reveals that the correlation is zero. By contrast, estimating the correlation between the true and estimated speciation rates after pooling all of the branches from the two trees conceals the poor performance of BAMM; we observe a positive relationship in the pooled sample simply because BAMM estimates higher average rates for one tree compared to the other! Clearly, summarizing the relationship between true and estimated branch-specific diversification rates without controlling for among-tree differences in diversification rates induces a spurious and positively misleading impression of the BAMM’s performance. We demonstrate this phenomenon with real trees from our simulation study (Figure 6); pooling results among trees conceals the poor performance of rate-parameter estimates within individual trees (Figure 7).

Figure 6. Pooling parameter estimates among trees can be misleading. We estimated branch-specific parameter estimates for two simulated trees using BAMM, and compared the true speciation rate for each branch (x-axis) against the estimated speciation rate for each branch (y-axis). Right panel. BAMM is intended to estimate the branch-specific parameters (speciation and extinction rates) of individual trees. When we appropriately estimate the correlation between the estimated and true branch-specific speciation rates for each tree separately, we observe that the slope of the regression for each tree is ~zero (i.e., the true variation and estimated variation in branch-specific speciation rates is uncorrelated). Left panel. When we assess the relationship between the true speciation rate and the estimated speciation rate (either by fitting a linear regression or computing Spearman’s correlation) across both trees, we observe a very strong (and completely spurious) positive relationship. This positive correlation is spurious because it is an artifact of the difference in the average (overall) speciation rate in these two trees.

Figure 7. The relationship between true and simulated speciation rates for a random sample of our trees. We plotted the correlation between the true and estimated speciation rates for 25 trees from our simulation study. The slopes indicate the correlation between estimated (y-axis) and true (x-axis) values of the branch-specific speciation rates. For many trees, the true variation and estimated variation in speciation rates is uncorrelated (the slope is zero) or even negatively correlated (c.f., Figure 3).

### 3.2. Summarizing absolute errors in diversification-rate parameter estimates is misleading

The relevant aspect of the branch-specific diversification-rate parameter estimates is their relative values. This is because we are interested in changes in the relative magnitude of diversification rates across branches of the study tree. For example, we are interested to discover that rates of speciation in the dolphin clade is three times higher than the speciation rate in their cetacean relatives (Steeman et al.,2009); by contrast, we are generally uninterested in differences in the absolute values of the dolphin and non-dolphin speciation rates. These considerations inform how we should assess error in estimates of the diversification-rate parameters. Imagine, for example, that the estimated speciation rate for a branch differes from the true speciation rate by an absolute value of 0.03. If the true speciation rate was 0.01, the absolute error (i.e., 0.03 – 0.01 = 0.02.) may seem small. However, this degree of estimation error implies a three-fold difference in the speciation rate! Presumably, such considerations informed Rabosky’s decision to use proportional error to characterize the ability of BAMM to estimate diversification-rate parameters in his study (Rabosky, 2014). Proportional error measures the relative error of parameter estimates, and is simply computed as the estimated value divided by the true value:

where refers to a particular parameter (e.g., the speciation or extinction rate) and refers to a particular branch; a proportional error near 1 indicates that the estimated value is close to the true value.

In our study, we followed Rabosky (2014), and used proportional error to assess the reliability of diversification-rate parameters using BAMM. The results of our simulation study—using the same summary to quantify estimation error as that used by Rabosky (2014)—demonstrate that BAMM does not provide reliable parameter estimates (Figure 3). For reasons that are unknown to us, the recent post by Rabosky attempts to discredit this conclusion by using a different (and misleading) procedure to summarize the results of our simulation. Specifically, Rabosky uses the mean-squared error to summarize the error in diversification-rate parameters. Mean-squared error measures the squared absolute difference between the true and estimated diversification rates, which is computed as:

As discussed above, assessing error in speciation and extinction rates (which are small absolute values) is both deeply misleading and completely irrelevant. Moreover, because the speciation and extintion rates are small, taking the square of the difference has the effect of further shrinking the apparent error (Figures 8 and 9). Accordingly, Rabosky appears to prefer the (correct) summary measure—the proportional error—unless it reveals that BAMM does not provide reliable parameter estimates. When this relevant summary exposes the fact that BAMM does not provide reliable diversification-rate parameter estimates, Rabosky suddenly searches for alternative (and inappropriate) summaries that conceal the flaws in BAMM.

Figure 8. Error in speciation-rate estimates using the “random” option implemented in BAMM. We analyzed the 100 rate-variable trees published in Moore et al. (2016) using the “combineExtinctionAtNodes = random” option with an expected number of diversification-rate shifts of 1. We computed the posterior mean of the speciation rate for each branch in each tree. We then plotted the error in the speciation-rate estimates using: (1) mean-squared error (left panel), and; (2) proportional error (right panel). The only relevant result is that depicted in the right column; the proportional error. A proportional error of 1 would imply that BAMM provides unbiased parameter estimates. By contrast, our results reveal that parameter values estimated using BAMM are, on average, quite biased.

Figure 9. Error in extinction-rate estimates using the “random” option implemented in BAMM. These results are summarized as described in the legend for Figure 8, but here pertain to error in extinction-rate estimates. The extinction-rate estimates, as expected, are more severely biased than speciation-rate estimates.

Finally, we also calculated the error in estimated diversification parameters using the “if_different” option for computing extinction probabilities (Figure 10). These results are qualitatively identical to those presented in our study (i.e., the results are the same for both the “random” and “if_different” options for computing extinction probabilities). Again, these results reveal that Rabosky’s claims are manifestly untrue: our main conclusion—BAMM does not provide reliable parameter estimates—is not an artifact option used to compute extinction probabilities.

Figure 10. Error in speciation- and extinction-rate estimates using the “if_different” option for computing extinction probabilities implemented in BAMM. The summaries follow those described in Figures 7 and 8, but are based on analyses using the default “if_different” option for computing extinction probabilities in BAMM. Contrary to recent claims by Rabosky, BAMM provides unreliable parameter estimates using either option for computing extinction probabilities.

## Summary

In this post, we have demonstrated that BAMM does not provide reliable estimates of branch-specific diversification-rate parameters (speciation and extinction rates). The results of our simulation—when summarized correctly—confirm that the estimated variation in rates is uncorrelated with the true variation in rates. We have also demonstrated that Rabosky’s claim—that he is “unable” to reproduce our results—is specious; the procedures that we used to summarize the results of our simulation exploring the performance BAMM are identical to the procedures used in his original study (Rabosky, 2014). Moreover, we have demonstrated that Rabosky’s claim—that our conclusions are “scientifically unjustified” because of the option we used to compute extinction probabilities in BAMM—is demonstrably untrue: the estimated variation in rates inferred using BAMM is uncorrelated with the true variation in rates using both the default “if_different” option and “random” option for computing extinction probabilities in BAMM. Finally, we have exposed Rabosky’s attempt to conceal the pathological statistical behavior of BAMM; he has attempted to discredit the conclusions of our study by manipulating our results using new—and completely invalid and misleading—summary procedures that conceal the poor performance of BAMM.

The demonstration that BAMM does not provide reliable estimates of branch-specific diversification-rate parameters should come as no surprise. The inability of BAMM to provide reliable estimates of the posterior distribution of branch-specific rate parameters is unsurprising, given our formal proofs that: (1) the likelihood function used by BAMM is incorrect (it cannot correctly compute the probability of the data when rates of diversification vary), and; (2) the CPP prior model used by BAMM is incoherent (it cannot correctly specify the probability distribution of rates across branches of the tree). The posterior distribution is an updated version of the prior distribution, where the prior is updated by the information in the data via the likelihood function. Given our demonstration that the likelihood function in BAMM is invalid, and our demonstration that the prior model in BAMM is invalid, it is unreasonable to expect that, somehow, the posterior estimates of the parameters will be reliable.

We wrote our paper because we believe that it is in the best interest of our community to expose the serious problems with BAMM. Our significance statment makes this clear: “Exposing the problems with BAMM is important both to empiricists (to avoid making unreliable inferences using this method) and to theoreticians (to focus their efforts on solving the problems that we identify).” We are dismayed by Rabosky’s response to our paper. Rather than trying to correct the problems afflicting BAMM (i.e., to develop a correct likelihood function and a coherent prior model), he has resorted to making demonstrably false statements about our study and attempting to dismiss our conclusions by manipulating our results. These efforts will not resolve the problems with BAMM; they can only serve to cause confusion regarding the problems with BAMM. In this regard, we are troubled that Rabosky is putting his own interests above those of our community.

## Data availability

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

## References

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. The Moore et al (2016) use of developer switch “random” was not scientifically justified. link

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.