The Estimation of Tree Posterior Probabilities Using Conditional Clade Probability Distributions

Principle of Conditional Independence

The approximation equation (5) is based on the assumption that clades in different parts of the tree are approximately independent. To specify this condition more explicitly, it is useful to consider unrooted trees. The same concepts apply to any rooted tree by considering the corresponding unrooted tree. Each edge of an unrooted tree T corresponds to a split which partitions the taxa into two sets. Let s be a split in T where the removal of the corresponding edge from T leaves behind rooted subtrees T₁ and T₂. I say that a split s separates splits s₁ and s₂ in tree T if s₁ corresponds to an edge in T₁ and s₂ corresponds to an edge in T₂, or vice versa. The principle of conditional independence of separated subtrees states that if s separates s₁ and s₂, then s₁ and s₂ are approximately independent so that P(s₁ ∩ s₂|s) ≈ P(s₁|s₂) P(s₂|s) and P(s₁|s, s₂) ≈ P(s₁|s). More generally, if each of the m splits a₁,...,a_m is separated from each of the n splits b₁,...,b_n by split s in a tree T, then splits a₁,...,a_m are approximately conditionally independent of splits b₁,...,b_n given s, and the following probability relationships are approximately true:

The following computational algorithm is based on this principle.

Computational Algorithm

In order to apply Equation (5), it is necessary to estimate various clade conditional probabilities from a sample of trees. An algorithm to calculate and store the necessary probabilities efficiently is described in online Appendix 1 but described briefly here. One pass through the sample of trees finds the unique trees and a count of how often each is sampled. In a pass through this summary, two maps are formed. Map 1 (m₁) counts for each clade that appears in at least one tree the number of trees that contain the clade (including clades of size 1 and the clade of all taxa, the totals of which equal the sample size). The second map (m₂) counts for each clade of size 2 or more and its two direct subclades the number of trees that contain this triple of clades. The algorithm also stores the total number of trees sampled by summing the counts.

Once these two maps are stored, the probability of a tree is estimated by taking the product over clades C of size 2 or more with subclades C₁ and C₂ the conditional probabilities P(C₁ ∩ C₂|C) ≈ m₂(C,C₁,C₂)/m₁(C). If a tree contains a pair of sister clades that are not in the map (and therefore were unsampled together), its probability is estimated to be 0.

Demonstration on Carnivora Example

The carnivora example consists of the cox1 gene for 12 taxa including six dog-like and six cat-like taxa. The sample was taken using MrBayes 3.2 (Ronquist et al. 2012) with the BEAGLE library (Ayres et al. 2012) under the HKY85 likelihood model (Hasegawa et al. 1985) with separate parameters for each codon position and using approximate gamma distributed rates (Yang, 1994) with four rate categories. The total sample size was 1.1 million with the first 100 000 trees discarded as burn in and the remaining million trees subsampled every 10 trees. Examination of trace plots of the log-likelihood and comparisons of split probabilities across independent runs are consistent with good mixing and sufficient burn in. There are 654 729 075 possible unrooted trees with 12 taxa. The sample contains only 229 of these trees. The most probable tree was sampled 33 925 times and 50 trees were sampled only once.

For each of these 229 trees, I calculate the probability of the tree first using the simple sample proportions and second using Equation (5). These 229 points are plotted in Figure 3 using both regular axes and logarithmic axes. There is strong agreement between the two calculations for both trees with relatively high and low probability. The correlation between the calculations is 0.997. The largest absolute difference between the two calculations is about 0.017 where the SRF estimate of the second most probable tree (0.1625) is slightly higher than the CCD estimate (0.1454). The sum of the 229 CCD estimated probabilities is 0.9985 which suggests that only about 0.15% of the total probability is missing from the sample. Figure 3b shows that 50 trees sampled only once have CCD probability estimates ranging from 3×10⁻⁸ to 1.4× 10⁻⁴, the larger probability more than 1000 times larger than the smaller, whereas each of these trees has an SRF probability estimate of 1×10⁻⁵.

Validation

I have shown that the CCD method based on the theoretical assumption of approximate conditional independence of separated subtrees is accurate on a single sample from a single data set. To provide additional evidence that this method is valid, I examine its behavior with repeated samples on the Carnivora data set, on a sample from a uniform distribution over tree space where exact calculations are possible, and for multiple other real data sets.

Carnivora example with repeated samples.—Using the Carnivora data set from the previous section, I added 9 more runs with different streams of random variables using the identical model and MCMC parameters and summarize the 10 total runs. The number of unique trees was 249.7 ±16.2 (mean ± standard deviation) and the estimated proportion of the total tree space covered by each individual sample was 0.9985±0.0002. Over all 10 runs, a total of 470 trees were sampled, so many trees appear in only one or a few samples. The correlation between the SRF and CCD probability estimates for trees sampled in each data set is 0.9970±0.0007, indicating strong agreement between these estimates for each sample. The SRF probabilities of the trees in Figures 1 and 2 are 0.3367 ±0.0092 and 3.00×10⁻⁶ ±4.83×10⁻⁶, respectively, when compared with the CCD probability estimates of 0.3411±0.0094 and 1.23×10⁻⁶ ±4.49×10⁻⁷. The CCD probability estimate of Figure 1 tree is slightly larger than the SRF estimate in each of the 10 DSs (0.0045±0.0029), which confirms that there is some bias in the CCD probability estimation of the true distribution. However, the sample to sample variability is about the same for both estimation methods and at this MCMC sample size, the size of the approximation bias is smaller than the standard deviation due to sampling variation for the most probable tree. Note that the SRF estimates of the tree in Figure 2 are based on three estimates of 1/100 000 and 7 of 0. For the low probability tree in Figure 2, the CCD probability estimates are substantially less variable. This suggests that using the CCD method introduces only a small amount of bias at little penalty in variance when estimating probabilities of highly probable trees, but greatly reduces the variance in estimating probabilities of more unlikely trees.

Figure 4a displays the relative standard deviations (standard deviation divided by the mean of the SRF estimates) of the CCD estimates versus the SRF estimates across the 10 samples for each of the 470 trees found in at least one sample. The line is drawn where these values are equal. Most of the points fall below the line, which means there is less sample to sample variation for the CCD method than the SRF method for almost all trees for this data set. This is consistent with the idea that clade probability estimates are estimated with less variation from samples than are tree probabilities. Thus, the approximations based on CCDs are much less variable from sample to sample for most of the trees with nonnegligible probability. The line of points to the right of 3 on the horizontal axis are from the trees sampled once in one sample and not at all in other samples.

Figure 4b plots the ratio of the absolute mean difference in the two estimates to the standard deviation of the SRF estimates against means of the SRF estimates for all trees sampled at least once with both variables plotted on a log scale to spread out the points near the origin. This display shows that the size of the absolute bias is smaller than the standard deviation of the SRF estimates for the majority of the trees. A large majority of the points fall below the dotted line at zero where the absolute bias due to the approximation is less than the size of the uncertainty due to sampling variability. The positive trend in this plot indicates that the size of the potential estimation bias relative to the size of the error due to sample to sample uncertainty tends to be smaller for very low probability trees than for more probable trees.

To summarize, replacing simple estimates with conditional clade probability estimates introduces a small amount of bias for each tree. However, the size of this bias is fairly small and often may be smaller than the MCMC uncertainty of the simple estimates from a single MCMC sample. For less probable trees in particular, the conditional clade probability estimates from a single sample are likely both more accurate and less variable than simple estimates from a single sample, at least for this single data set.

Samples from a uniform distribution.—To explore the behavior of the CCD method in a case where the probability distribution is diffuse and the true probabilities are known exactly, I took an independent sample of 100 000 trees from a uniform prior distribution on the space of unrooted trees with 10 taxa, of which there are 2 027 025. This sample size was large enough to estimate the clade distributions well, but too small for the SRF method to be accurate. Each tree has true probability about 4.9×10⁻⁷, but a tree sampled once will have SRF probability of 1.0×10⁻⁵ and some multiple of this if it is sampled multiple times. There were 97 722 unique trees included in the sample, one of which was sampled four times. The SRF probability estimates for each of these trees is far larger than the true probability. In contrast, the CCD probability estimates for these sampled trees vary from 6×10⁻⁸ to 1.6×10⁻⁶ and the median estimate is 5.2×10⁻⁷ which is very close to the true value. The estimate of the total coverage of the sample (sum of all approximate probabilities) is 0.0524 which is close to the correct value of 97 722/2.027.025 ≐ 0.04821.

To further explore the accuracy of the CCD approximation of this true uniform distribution, I took a second random sample of 100 000 trees and used the estimated conditional clade distributions from the first sample to estimate probabilities of trees in the second sample. The second sample had 97 516 unique trees, only about 5% of which were also in the first sample. The estimated probabilities of these trees ranged from 0 to 1.6×10⁻⁶ and the median estimated probability was 4.8×10⁻⁷, very close to the true value of 4.9×10⁻⁷. Only 36 of these trees had estimated probabilities of zero based on the conditional clade distributions of the first sample; the lowest estimated probability for those for which it was positive was 5×10⁻⁸. The mean and standard deviation of the CCD estimated probabilities was 4.9×10⁻⁷ ±1.5×10⁻⁷ and was centered over the true value. The sum of estimated probabilities of trees in this second sample is 0.0480, very close to the true value of 97 516/2 027 025 ≐ 0.0481.

This example shows that the CCD estimates can be very accurate, even when applied to trees that are not in the original sample. In this instance where the true probability of the most probable tree is substantially less than the reciprocal of the sample size, the CCD estimates are much more accurate and less variable from sample to sample than are the SRF estimates.

A much larger data set.—Poux et al. (2006) constructed and analyzed a data set using three unrelated genes from 62 mammalian taxa that included two marsupial outgroups, representatives of all extant placental mammalian orders, and heavy representation of South American rodents and primates, the primary foci of the study. I used a much simpler model than the authors and treated all 3768 sites in the alignment identically using the general time reversible (GTR) model (Tavaré, 1986) with approximate gamma distributed rates using four rate categories.

Using MrBayes 3.2 with the BEAGLE library I ran MCMC for 5 500 000 generations, discarded the first 500 000 as burn in, and subsampled every 100th tree, resulting in a sample of 50 000 trees. Comparisons split probability estimates with a second run of the same size support the conclusions that the size of this sample was sufficiently large and that the chain mixed well (data not shown). Figure 5 shows a graph of the simple and conditional clade estimates of the probabilities of each of the 4297 sampled trees. The correlation between these two measures is 0.990. The sum of the conditional clade probability estimates is 0.9150, so this sample likely represents over 90% of the total probability, but about 8.5% of the probability resides on unsampled trees. The largest absolute difference between the two estimates is 0.008. The most probable tree has SRF estimate 0.0655 and CCD estimate 0.0604. This example demonstrates that the CCD method can be accurate even when tree space is enormous; there are about 6.97 ×10⁹⁸ unrooted trees with 62 taxa, much greater than the estimated 10⁸⁰ atoms in the observable universe. A sample of 50 000 trees cannot estimate all of these probabilities accurately, but the CCD method provides a means to estimate accurately the probabilities of the sampled trees and for many relatively probable trees not included in the sample.

Additional data sets.—Lakner et al. (2008) and Höhna and Drummond (2012) examine the performance of many MCMC sampling procedures across 11 DSs with various numbers of taxa and nucleotide sites. Further details about these data sets are included in these papers and in Table 2. For each of these data sets, I used MrBayes 3.2 with the BEAGLE library treating all sites as identically distributed with the GTR model with approximate gamma distributed rates using four rate categories to sample 5 500 000 trees, discarding the first 500 000 as burn in and subsampling each 50th tree of the remainder for a sample of 100 000 trees from each data set. The results are summarized in Table 3. A second set of runs under the same conditions, but with different random numbers, shows very similar results, indicating that these MCMC samples are likely not to suffer from poor convergence (data not shown).

The comparisons between the SRF and CCD methods vary across these data sets. I will discuss the results separately for groups of data sets for which the behavior is somewhat similar.

DS 2, DS 3, and DS 7 are somewhat similar to the 62-taxon mammal DS from the previous section in that the sample of trees contains a large fraction of the total posterior distribution (about 75% or higher) and the correlation between the CCD and SRF estimates of tree probabilities within the sample is greater than 0.97. There is generally very strong agreement between the estimates, as seen in Figure 6b for DS 2.

DS 1 and DS 4 have estimated coverage between 60% and 70% and lower correlations between the estimates (0.777 and 0.606). In these two DSs, there is a fair amount of variation in the CCD estimates among the many trees with relatively low SRF estimates, which lowers the correlation between these measures when compared with the data sets with larger coverage. In addition, DS 1 has some quirks where there are several trees that are sampled only once, but have CCD estimates much larger than 1 in 100 000. Many of the points are below the line in Figures 6a and c due to the fact that the coverage is substantially below 1.

DS 6 and DS 8 have coverage between 2% and 20%, maximal SRF estimates between 0.0001 and 0.0002 as no trees are sampled more than 20 times out of 100 000, and maximal CCD estimates about the same. There is significant positive correlation between the two methods for computing tree probabilities, but the association is weaker than the previous examples because the SRF method simply does not discriminate much among the sampled trees. In contrast, the CCD method varies over several orders of magnitude among trees sampled the same number of times and is measuring the posterior probabilities of these trees much more precisely.

The total coverage of the sample according to CCD for DS 5, DS 9, DS 10, and DS 11 is estimated to be much less than 1%. In these cases, the posterior distribution is very diffuse and even the most probable trees individually have very low probabilities. In some cases, the coverage is less by many orders of magnitude. For example, in DS 9 and DS 11, the coverage of the entire samples is estimated by the CCD method to be on the order of 10⁻¹² and 10⁻¹⁵, far smaller than the 10⁻⁵ estimated probability for almost all single trees using the SRF method. Figure 6f shows the comparison between SRF and CCD estimates for DS 10. The CCD estimates are close to zero for nearly all trees and are not substantially higher for trees sampled more than once compared with those sampled only once. Consequently, the correlation between the measures is nearly zero. This reflects that the SRF estimates are quite inaccurate and highly variable. In fact, for DS 11, a second random sample contains 99 790 distinct trees, none of which are in the first sample. For this data set, samples of more than 10¹⁶ trees would be required for even the most probable trees to be sampled with high probability more than a few times. Clearly, when the posterior distribution is so diffuse that even the most probable trees have smaller probability than one over the sample size, the SRF estimates of tree probabilities are far too large for trees in the sample. However, clade and conditional clade probabilities can be estimated more accurately and the CCD method provides a means to approximate tree probabilities from these distributions.

To further illustrate this point, DS 11 contains 71 taxa, and the true unrooted tree will have 68 nontrivial splits (those that do not separate just a single taxon from the rest). In the single MCMC sample of 100 000 trees, there were 21 splits with probability greater than 0.99, 29 with probability greater than 0.9, and 40 with probability greater than 0.5. The majority rule consensus tree has a fair amount of resolution, even though a lot of uncertainty remains in the tree. Many additional splits are probable enough that trees that contain them appear many times in the sample. The probabilities of these splits can be estimated with relatively small sampling error and the CCD method uses these accurate clade probability estimates to produce accurate tree probabilities when the SRF method fails.

To summarize, in cases where the MCMC sample represents a substantial fraction of the total posterior distribution, there is strong agreement between SRF and CCD estimates of tree probabilities. In cases where the coverage is low, SRF estimates for individual trees are unreliable, and so correlations between SRF and CCD estimates are weak.

Some Mathematical Results

I have presented the CCD method for estimating tree probabilities from conditional clade distributions on the basis of an assumption of approximate conditional independence among clades in separated subtrees and have shown close agreement between these estimates and SRF estimates when MCMC samples are large enough for SRF estimates to be accurate for many data sets. I justified the method with an appeal to a property on unrooted trees, but described the computational algorithm with regard to rooted trees. In addition, I have made other statements about mathematical properties of the CCD method without justification. The following list summarizes mathematical results about the CCD method for estimating tree probabilities. The rigorous mathematical arguments to support these results are described in online Appendix 2.

1. The CCD method for estimating tree probabilities of unrooted trees is invariant to different rootings of the trees.

2. The computational method in Equation (5) for rooted trees is equivalent to the following expression for unrooted trees:

   (6)

   where the node n is incident to 3 edges that correspond to the splits s₁(n), s₂(n), and s₃(n).

3. For any full specification of CCP distributions, the sum of the corresponding probabilities in Equation (5) taken over all possible rooted trees for the given set of taxa equals 1; hence, the CCD method defines a valid probability distribution on the set of all possible trees.

A derivation of the result that approximate conditional independence of separated subtrees follows from the form of the calculation of the likelihood in phylogenetic models is shown in online Appendix 3. A theoretical argument that smaller samples are needed to estimate probability distributions on trees that follow the conditional independence of separated subtree principle than to estimate general tree distributions is described in online Appendix 4.

Software

I have written C++ code that implements the concepts of this article. The code is available currently on my web site for download under the name ccd-probs (Larget, 2012) and is available at sourceforge.net/projects/ccdprobs. The software runs from a command line interface and includes separate components to: (i) summarize raw MrBayes output files with the mbsum program from Larget et al. (2010) and (ii) apply the CCD method to estimate the probabilities of all sample trees. The program creates an additional output file with the single split and joint triple split maps that can in principle be read by other programs that wish to apply the CCD probability distribution from a given set of trees without reprocessing the sample summaries. In addition, the program can take as input a second file of trees that are not part of the data for forming the maps, but whose probability can also be evaluated. The software uses the dynamic bitset functionality of the BOOST C++ library (Siek and Allison, 2012) to efficiently store clades.

Ramifications

Interpretations of consensus trees.—The results in this article will not drastically change the basic interpretations found in typical Bayesian phylogenetic analyses, but do offer modest improvements. Authors usually summarize a Bayesian posterior distribution on trees by reporting a consensus tree of the sample and annotating edges in this tree with posterior probabilities. This tree is simply a graphical way to display the collection of most probable clades, and one point of this article is to state that these clade probabilities are generally estimated accurately in MCMC samples, even when the probabilities of specific trees may not be. The CCD method in the article, however, will provide users of Bayesian phylogenetic methods a new tool for assessing how thoroughly their MCMC sample covers tree space by calculating the sum of CCD probabilities over all sampled trees. Particularly when there are many taxa, there can be large differences from data set to data set in how diffuse the posterior distribution is. I am not aware of any previous methods that allow one to estimate the fraction of the posterior distribution represented by the posterior sample in Bayesian phylogenetics. The ability to calculate this coverage may improve interpretations of results of Bayesian analyses. Some authors summarize posterior distributions by reporting the number of trees in the sample needed to obtain a high percentage, such as 90, 95, or 99, of the posterior sample. The work in this article shows that in many instances, the number of trees needed to include a high percentage of the sample may be a gross underestimate of the number of trees needed to cover the same percentage in the true posterior distribution.

Note that the previous statement does not imply that MCMC samples with low coverage have mixed poorly or lead to poor inferences. Small samples can lead to accurate estimates. Furthermore, most interesting phylogenetic hypotheses relate to clades and not to the complete specification of the tree. Even when there is significant uncertainty in the complete specification of the tree, there is often strong information about the posterior probability of many clades.

In addition, note that the notion of coverage introduced in this article (the total posterior probability of all sampled trees) is distinct from the frequentist notion of “coverage probability” (the probability that a confidence interval contains the true value of a parameter).

Potential improvements in BEAST.—Höhna and Drummond (2012) created their CCP method for estimating tree probabilities from conditional clade distributions in order to speed up the mixing rate of MCMC proposal distributions and implemented the method in BEAST. The idea is that an initial standard MCMC run can be performed in order to estimate conditional clade probabilities and these probabilities can then be used by subsequent proposal methods to propose new trees from a distribution close to the true posterior distribution rather than uniformly at random like many currently implemented MCMC proposal methods. The CCD method in this article does not require renormalization, and thus ought to be substantially faster. In addition, the assumption of conditional independence among separated subtrees appears to be reasonably benign in real data sets which suggests that the CCD method may come closer than the CCP method of Höhna and Drummond (2012) to accurately estimate true posterior distributions. I have not tested either of these claims rigorously; the speed claim appears self-evident because the calculations require nearly identical work up to renormalization, and the renormalization step of the CCP method is the most time-consuming step by far. The only evidence I have that the CCD method is more accurate than the CCP method for estimating tree probabilities is that the CCD method will be exactly correct under the reasonable approximating assumption of conditional independence in separated subtrees and the CCP method is based on no such theory. In addition, after renormalization, the CCP method will place a total probability of one on the set of all trees it considers. For many data sets, the true sum of these probabilities can be much smaller than one when the set of trees is small relative to the size of tree space. It is possible that replacing algorithms in BEAST that use their CCP method with an algorithm based on the CCD method in this article will improve the efficiency of MCMC sampling in this program. Similar ideas could be implemented in MrBayes.

A new class of importance sampling algorithms.—In conjunction with a method to propose branch lengths given a tree topology, the CCD method of this article could be used to generate trees at random from a distribution that mimics the true posterior distribution after an initial calculation to estimate conditional clade distributions. Independent samples from this proposal distribution could be used as the basis for importance sampling from the posterior distribution. If the CCD-based proposal distribution is an accurate enough representation of the true posterior distribution, importance sampling could replace MCMC and allow for independent samples from the posterior to be taken efficiently. Note that MCMC would still be required to set up the conditional clade probability distributions. Much work remains to develop these rough ideas into practical algorithms and tested software, but there is the distinct possibility that CCD-based importance sampling could greatly speed up Bayesian sampling from posterior distributions on phylogenetic trees, allowing for significant computational time savings for large data sets.

Improved accuracy in BCA.—BCA is sensitive to the estimation of small tree probabilities from data in some genes when those trees are highly probable in other genes. The CCD method in this article can be used to improve the accuracy of BCA in measuring the probability of discordance among gene trees by providing a means to estimate the probabilities of unsampled trees.

Alternative sampling methods for large trees.—If the assumption of conditional independence among separated subtrees is not too inaccurate, this suggests that Bayesian analysis of large trees could begin by partitioning taxa into groups that are likely to be separated on the true tree, carrying out separate analyses for these groups, and combining the results under the assumption of conditional independence. The computational effort of analyzing multiple data sets with smaller taxon sets lends itself to parallel computing and speedups due simply to working on smaller subtrees. For example, for the Carnivora data set from the beginning of the aricle (or a larger similar example), it could be possible to run separate analyses for the cat-like and dog-like taxa, perhaps with a few cross representatives each and then combine the results to obtain a sample from the larger tree space. A great deal of work is required to explore the potential gains in computational efficiency and to protect against bias from the conditional independence approximation. However, the good fit between true distributions and approximations that assume conditional independence suggests that some sort of the divide and conquer approach for large trees may work well without introducing too much bias.

Summary

The CCD method for estimating tree probabilities from CCDs solves a long-standing problem among developers of phylogenetic methods on how to use clade information to better estimate probabilities of trees. The method provides a new tool to augment standard Bayesian phylogenetic analyses by allowing the scientist to estimate the proportion of tree space covered by the MCMC sample. The method will patch a hole in BCA when implemented in the program BUCKy. More importantly, the method has the potential, when further developed and combined with other ideas, to lead to new computational approaches for Bayesian phylogenetic inference that could be substantially more efficient than the current state of the art.