Exact and approximate limit behaviour of the Yule tree's cophenetic index
Krzysztof Bartoszek
TL;DR
This paper investigates the limit distribution of the cophenetic index in Yule trees, providing theoretical results, approximations, and practical methods for hypothesis testing of phylogenetic trees.
Contribution
It establishes the almost sure convergence of the normalized cophenetic index and introduces heuristic simulation methods for the limit distribution.
Findings
Normalized index converges almost surely and in L2.
Limit distribution is similar to Quicksort distribution in certain cases.
Heuristic simulation methods produce reasonable tail approximations.
Abstract
In this work we study the limit distribution of an appropriately normalized cophenetic index of the pure-birth tree conditioned on contemporary tips. We show that this normalized phylogenetic balance index is a submartingale that converges almost surely and in L2. We link our work with studies on trees without branch lengths and show that in this case the limit distribution is a contraction-type distribution, similar to the Quicksort limit distribution. In the continuous branch case we suggest approximations to the limit distribution. We propose heuristic methods of simulating from these distributions and it may be observed that these algorithms result in reasonable tails. Therefore, we propose a way based on the quantiles of the derived distributions for hypothesis testing, whether an observed phylogenetic tree is consistent with the pure-birth process. Simulating a sample by the…
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8| Centring () | |||
| Scaling () | |||
| Centring () | |||
| Scaling () |
| Model | Sackin’s | Colless’ | ||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| mean() | variance() | |||||||||||
| Yule | ||||||||||||
| Coalescent | ||||||||||||
| birth–death | ||||||||||||
| birth–death | ||||||||||||
| biased speciation | ||||||||||||
| biased speciation | ||||||||||||
| biased speciation | ||||||||||||
| biased speciation | ||||||||||||
| biased speciation | ||||||||||||
| biased speciation | ||||||||||||
| biased speciation | ||||||||||||
| biased speciation | ||||||||||||
| Sackin’s | Colless’ | ||||||||
|---|---|---|---|---|---|---|---|---|---|
| Sackin’s | |||||||||
| Sackin’s | Colless’ | ||||
|---|---|---|---|---|---|
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Exact and approximate limit behaviour of the Yule tree’s cophenetic index
Krzysztof Bartoszek
Abstract
In this work we study the limit distribution of an appropriately normalized cophenetic index of the pure–birth tree conditioned on contemporary tips. We show that this normalized phylogenetic balance index is a submartingale that converges almost surely and in . We link our work with studies on trees without branch lengths and show that in this case the limit distribution is a contraction–type distribution, similar to the Quicksort limit distribution. In the continuous branch case we suggest approximations to the limit distribution. We propose heuristic methods of simulating from these distributions and it may be observed that these algorithms result in reasonable tails. Therefore, we propose a way based on the quantiles of the derived distributions for hypothesis testing, whether an observed phylogenetic tree is consistent with the pure–birth process. Simulating a sample by the proposed heuristics is rapid, while exact simulation (simulating the tree and then calculating the index) is a time–consuming procedure. We conduct a power study to investigate how well the cophenetic indices detect deviations from the Yule tree and apply the methodology to empirical phylogenies.
Keywords : Contraction type distribution; Cophenetic index; Martingales; Phylogenetics; Significance testing
1 Introduction
Phylogenetic trees are now a standard when analyzing groups of species. They are inferred from molecular sequences by algorithms that often assume a Markov chain for mutations at the individual positions of the genetic sequence (e.g. Ewens and Grant, 2005; Felsenstein, 2004; Yang, 2006). Given a phylogenetic tree it is often of interest to quantify the rate(s) of speciation and extinction for the studied species. To do this one commonly assumes a birth–death process with constant rates. However, the development of formal statistical tests whether a given tree comes from a given branching process model is an open area of research (see the still relevant “Work remaining” part at the end of Ch. in Felsenstein, 2004). The reason for the apparent lack of widespread use of such tests (but see Blum and François, 2005) could be the lack of a commonly agreed on test statistic. This is as a tree is a complex object and there are multiple ways in which to summarize it in a single number.
One proposed way of summarizing a tree is through indices that quantify how balanced it is, i.e. how close is it to a fully symmetric tree. Two such indices have been with us for many years now: Sackin’s (Sackin, 1972) and Colless’ (Colless, 1982). Alternatively, McKenzie and Steel (2000) proposed to measure balance by counting cherries on the tree and they showed that after appropriate centring and scaling, this index converges to the standard normal distribution (for examples of other indices see Ch. in Felsenstein, 2004).
Recently, a new balance index was proposed—the cophenetic index (Mir et al., 2013). The work here is inspired by private communication with evolutionary biologist Gabriel Yedid (current affiliation Nanjing Agricultural University, Nanjing, China) concerning the usage of the cophenetic index for significance testing of whether a given tree is consistent with the pure–birth process. He noticed that simulated distributions of the index have much heavier tails than those of the normal and t distributions and hence, comparing centred and scaled cophenetic indices with the usual Gaussian or t quantiles is not appropriate for significance testing. It would lead to a higher false positive rate—rejecting the null hypothesis of no extinction when a tree was generated by a pure–birth process.
Our aim here is to propose an approach for working analytically with the cophenetic index, especially to improve hypothesis tests for phylogenetic trees, i.e. how to recognize if the tree is out of the “Yule zone” (Yang et al., 2017). We show that there is a relationship between the cophenetic index and the Quicksort algorithm. This suggests that the methods exploring (e.g. Fill and Janson, 2000, 2001; Janson, 2015) the limiting distribution of the Quicksort algorithm can be an inspiration for studying analytical properties of the cophenetic index.
The paper is organized as follows. In Section 2 we formally define the cophenetic index (for trees with and without branch lengths) and present the most important results of the manuscript. We define an associated submartingale that converges almost surely and in (Thm. 2.4), propose an elegant representation (Thm. 2.7) and a very promising approximation (Def. 2.8). Afterwards in Section 3, we show that in the discrete setting the limit law of the normalized cophenetic index is a contraction–type distribution. Based on this we propose alternative approximations to the limit law of the normalized (with branch lengths) cophenetic index. In Section 4 we describe heuristic algorithms to simulate from these limit laws, show simulated quantiles, explore the power of the cophenetic index to recognize deviations from the Yule tree (comparing with Sackin’s and Colless’ indices’ powers), and apply the indices to example empirical data. In Section 5 we prove the claims presented in Section 2 alongside other supporting results. Then, in Section 6 we study the second order properties of this decomposition and conjecture a Central Limit Theorem (CLT, Rem. 6.10). We end the paper with Section 7 by describing alternative representations of the cophenetic index.
2 The cophenetic index and summary of main results
Mir et al. (2013) recently proposed a new balance index for phylogenetic trees.
Definition 2.1** (Mir et al. (2013))**
For a given phylogenetic tree on tips and for each pair of tips let be the number of branches from the root to the most recent common ancestor of tips and . We then the define the discrete cophenetic index as
[TABLE]
Mir et al. (2013) show that this index has a better resolution than the “traditional” ones. In particular the cophenetic index has a range of values of the order of while Colless’ and Sackin’s ranges have an order of . Furthermore, unlike the other two previously mentioned, makes mathematical sense for trees that are not fully resolved (i.e. not binary).
In this work we study phylogenetic trees with branch lengths and hence consider a variation of the cophenetic index.
Definition 2.2
For a given phylogenetic tree on tips and for each pair of tips let be the time from the most recent common ancestor of tips and to the root/origin (depending on the tree model) of the tree. We then define the continuous cophenetic index as
[TABLE]
Remark 2.3
In the original setting, when the distance between two nodes was measured by counting branches, Mir et al. (2013) did not consider the edge leading to the root. In our work here, where our prime concern is with trees with random branch lengths, we include the branch leading to the root. This is not a big difference, one just has to remember to add to each distance between nodes the same exponential (—parametrization by the rate) random variable (see Section 5 for description of the tree’s growth).
The results of the present manuscript are built around a scaled version of the cophenetic index which is an almost surely and convergent submartingale. We first introduce some notation. Let be the –algebra containing all the information on the Yule with tips tree and define
[TABLE]
Below we present the main results concerning the cophenetic index, leaving the proofs and supporting theorems for Section 5.
Theorem 2.4
Consider a scaled cophenetic index
[TABLE]
* is a positive submartingale that converges almost surely and in to a finite first and second moment random variable.*
Definition 2.5
For let us define as the indicator random variable taking the value of if a randomly sampled pair of species coalesced at the –th (counting from the origin of the tree) speciation event.
We know (e.g. Bartoszek and Sagitov, 2015b; Stadler, 2009; Steel and McKenzie, 2001) that
[TABLE]
Definition 2.6
For let us introduce the random variable
[TABLE]
Theorem 2.7
* can be represented as*
[TABLE]
where are i.i.d. random variables.
Definition 2.8
Define the random variable as
[TABLE]
where are i.i.d. random variables.
Remark 2.9
Despite the apparent elegance, it is not straightforward to derive a Central Limit Theorem (CLT) or limit statements concerning from the representation of Eq. (3). Initially one could hope (based on “typical” results on limits for randomly weighted sums, e.g. Thm. of Rosalsky and Sreehari, 1998) that could converge a.s. to a random variable that has the same limiting distribution as .
Similarly, as in the proof of Thm. 2.4 in Section 5, because , we have that is an bounded submartingale
[TABLE]
Hence, converges almost surely. Figure 1 can easily mislead one to believe in the equality of the limiting distributions of and . However, in Thm. 6.8 we can see that and convergence to different limits. Therefore, and cannot converge in distribution to the same limit. However, as we shall see in Section 4, provides a reasonable approximation (and importantly extremely cheap, in terms of computational time and memory) to in the sense of their distributions.
Definition 2.10
We naturally define the scaled discrete cophenetic index as
[TABLE]
Theorem 2.11
* is an almost surely and convergent submartingale.*
The applied reader will be most interested in how the results here can be practically used. As written already in the Introduction balance indices are often used to provide a single–number summary of the tree’s shape. Such statistics can be then used e.g. to test if the tree is consistent with some null model (here the Yule model). Naturally, there has been extensive work on using different balance indices for significance testing (e.g. Agapow and Purvis, 2002; Blum and François, 2005; Yang et al., 2017). However previous works nearly always worked with indices that only considered the topology and often obtained the rejection regions through direct simulations.
Unfortunately, looking only at the tree’s topology will not allow for distinguishing between some models. In particular (as seen in Tab. 3) there is no difference (from the topological indices perspective) between a Yule tree, a constant rate birth–death tree and a coalescent tree. Hence, a temporal index that also takes into account the branch lengths should be used (as indicated in the “Work remaining” section at the end of Ch. in Felsenstein, 2004). A statistic based on performs significantly better (but in these cases still leaves a lot to be desired). However, shows it true usefulness when employed to distinguish a biased speciation (Blum and François, 2005) from a Yule model. Blum and François (2005) indicated that there is a regime where topological indices fail completely. Table 3 shows that in this setup (and also certain others) the temporal index in superior in recognizing the deviation from the Yule tree.
Directly simulating a tree from a null model (Yule here) and then calculating the index will of course give a sample from the correct null distribution. However, this approach is costly both in terms of time and memory. Therefore, if theoretical results that provide equivalent, asymptotic or approximate representations of the index’s law are available they could speed up any study by orders of magnitude. In fact this is clearly visible in Tab. 1, calculating the cophenetic index directly from a sample of simulated pure–birth trees is over times slower than considering . Even more dramatically one can obtain a sample from an approximation to the equivalent representation of the asymptotic distribution of (after normalization) nearly times faster than directly sampling the discrete cophenetic index.
In Alg. 1 we describe how the presented here approach can be used for significance testing. Then, in Section 4 we discuss in detail the required computational procedures, present simulation results concerning the power of the tests and apply the tests to empirical data. Preceding this computational Section is a characterization of the limit distribution of (normalized) and another proposal to approximate the limit of (normalized) . This section justifies the described simulation algorithms in Section 4.
Theorem 2.12
A random variable with subscript NRE (no root–edge) indicates that this random variable comes from a tree lacking the edge leading to the root.
[TABLE]
[TABLE]
Proof The proof of the expectation part is due to Mir et al. (2013); Sagitov and Bartoszek (2012). The variance of is due to Cardona et al. (2013); Mir et al. (2013). The variance of is a consequence of the lemmata and theorems presented in Section 6. When the root edge is not included, then we have to decrease the expectation by . This is due to each pair of tips “having” the root edge included in the cophenetic distance between them. In the case of branch lengths, the expectation of the root edge, distributed as , is one. Without a root edge for the same reason the variance of has to be decreased by . In the discrete case the root edge has a deterministic length of and hence no effect on the variance.
3 Contraction–type limit distribution
Even though the representation of Eq. (3) is a very elegant one, it is not obvious how to derive asymptotic properties of the process from it (compare Section 6). We turn to considering the recursive representation proposed by Mir et al. (2013)
[TABLE]
where and are the number of left and right daughter tip descendants. Obviously .
From Eq. (8) we will be able to deduce the form of the limit of the process. In the case with branch lengths we attempt to approximate the cophenetic index with the following contraction–type law
[TABLE]
where , are independent random variables (we index with the mean to avoid confusion with , Section 5, the time between the second and third speciation event which is also distributed). These are the branch lengths leading from the speciation point. The rationale behind the choice of distribution is that a randomly chosen internal branch of a conditioned Yule tree with rate is distributed (Cor. and Thm. Stadler and Steel, 2012). This is of course an approximation, as we cannot expect that the laws of the branch lengths with the depth of the recursion should become indistinguishable from the law of the average branch. In fact, we should expect that the law of Eq. (9) has to depend on , i.e. the level of the recursion. For larger the branches have distributions concentrated on smaller values, e.g. compare the randomly sampled root adjacent branch length law (Thm. Stadler and Steel, 2012) with the law of the average branch length.
However, as we shall see simulations indicate that approximating with the average law still could still yield acceptable heuristics, but not as good as the approximation by . We use the notation , to differentiate from , where the root branch is included, i.e.
[TABLE]
Define now
[TABLE]
and using Eqs. (6) and (7) we obtain the following recursions
[TABLE]
and
[TABLE]
The process is related to the process as
[TABLE]
In the continuous case we do not have an exact equality, we rather hope for
[TABLE]
in some sense of approximation. Hence, knowledge of the asymptotic behaviour of , will immediately give us information about , in the obvious way
[TABLE]
The processes , look very similar to the scaled recursive representation of the Quicksort algorithm (e.g. Rösler, 1991). In fact, it is of interest that, just as in the present work, a martingale proof first showed convergence of Quicksort (Régnier, 1989), but then a recursive approach is required to show properties of the limit. The random variable weakly and as weak convergence is preserved under continuous transformations (Thm. , p. Grimmett and Stirzaker, 2009) we will have weakly. Therefore, we would expect the almost sure limits to satisfy the following equalities in distribution (remembering the asymptotic behaviour of the expectations)
[TABLE]
and
[TABLE]
where is uniformly distributed on , , and are identically distributed random variables, so are , and , and , , and are independent. Following Rösler (1991)’s approach it turns out that the limiting distributions do satisfy the equalities of Eqs. (10) and (11).
Let be the space of distributions with zero first moment and finite second moment. We consider on the Wasserstein metric
[TABLE]
Theorem 3.1
Let and assume that , , and are all independent. Define transformations , as
[TABLE]
and
[TABLE]
respectively. Both transformations and are contractions on and converge exponentially fast in the –metric to the unique fixed points of and respectively.
Remark 3.2
The proof of Thm. 3.1 is the same as Rösler (1991)’s proof of his Thm. . However, compared to the Quicksort algorithm (Rösler, 1991) we will have a upper bound on the rate of decay instead of . This speed–up should be expected as we have and instead of and . Thm. 3.1 can also be seen as a consequence of Rösler (1992)’s more general Thms. and . The rate of convergence is also a consequence of the general contraction lemma (Lemma 1, Rösler and Rüschendorf, 2001).
Now, using Lemmata 7.1, 7.2 (their proofs in Appendix B: Counterparts of Rösler (1991)’s Prop. for the cophenetic index differ only in detail from the proof of Prop. in Rösler, 1991) and arguing in the same way as Rösler (1991) did in his Section , especially his proof of his Thm. we obtain that and converge in the Wasserstein –metric to and whose laws are fixed points of and respectively. A minor point should be made. Here, we will have instead of in a counterpart of Rösler (1991)’s Prop .
Remark 3.3
One may directly obtain from the recursive representation that , and . We can therefore, see that in the discrete case the variance agrees. However, in the continuous case we can see that it slightly differs
[TABLE]
Remark 3.4
One can of course calculate what the mean and variance of , should be so that and . We should have and . This, in particular, means that these branch lengths cannot be exponentially distributed. We therefore, also experimented by drawing , from a gamma distribution with rate equalling and shape equalling . However, this significantly increased the duration of the computations but did not result in any visible improvements in comparison to Tab. 1.
4 Significance testing
4.1 Obtaining the quantiles
Algorithm 1 requires knowledge of the quantiles of the underlying distribution in order to define the rejection region. Unfortunately, an analytical form of the density of any scaled cophenetic index is not known so one will have to resort to some sort of simulations to obtain the critical values. Directly simulating a large number of pure–birth trees can take an overly long time, measured in minutes (on a modern machine with a large amount of memory, or hours on an older one). Fortunately, the cophenetic index can be calculated in time (Cor. 3 Mir et al., 2013) and such a tree–traversing algorithm was employed to obtain and . On the other hand, the suggestive (but wrong) approximations of Eq. (4) and contraction limiting distributions Eqs. (10) and (11) are significantly faster to simulate, see Tab. 1.
Simulating from the approximate Eq. (4) is straightforward. One simply draws independent random variables. Simulating random variables satisfying Eqs. (10) and (11) is more involved and it may be possible to develop an exact rejection algorithm (cf. Devroye et al., 2000). Here, we choose simple, approximate but still effective, heuristics in order to demonstrate the usefulness of the approach for significance testing.
We now describe algorithms (Algs. 2 and 3) for simulating from a more general distribution, , that satisfies
[TABLE]
where , are independent, , is some random vector, and for some appropriate that depends on ’s dimension. Of course in our case here we have , , ,
[TABLE]
and
[TABLE]
for , respectively. Of course, , are independent and distributed. If one considers also the root edge, then to the simulated random variable one needs to add when simulating or appropriately if one considers .
The recursion of Alg. 3 for a given realization of and random variables can be directly solved. However, from numerical experiments implementing Alg. 3 iteratively seemed computationally ineffective.
In Tab. 1 we report on the simulations from the different distribution. For each distribution we draw a sample of size and repeat this times. We compare the quantiles from the different distributions. We can see that the approximation of for is a good one and can be used when one needs to work with the distribution of the cophenetic index with branch lengths. In the case of the discrete cophenetic index we have found an exact limit distribution which is a contraction–type distribution. Therefore, one can relatively quickly simulate a sample from it without the need to do lengthy simulations of the whole tree and then calculations of the cophenetic index. Unfortunately, this contraction approach does not seem to give such good results in the Yule tree with branch lengths case. We used an approximation when constructing the contraction. Instead of taking the law of the length of two daughter branches, we took the law of an random internal branch. This induces a difference between the tails of the distributions that is clearly visible in the simulations. Even at the second moment level there is a difference. We calculated (Thm. 6.8) that , while . Therefore, the approximation by seems better already at the second moment level. Generally if one cannot afford the time and memory to simulate a large sample of Yule tree, simulating values seems an attractive option, as the discrepancy between the two distributions seems small.
In Fig. 2 we compare the density estimates of (scaled and centred) both continuous and discrete branches cophenetic indices and their respective contraction–type limit distributions. The density estimates generally agree but we know from Tab. 1 that for this is only an approximation. We simulated Yule trees and hence we report only the quantiles between and . Quantiles further out in the tails seemed less accurate and hence are not included in the table. Similarly, we can see less correspondence between the different estimates of kurtosis. This statistic relies on fourth moments and hence is more sensitive to the tails. On the other hand we can see much greater Monte Carlo error for the kurtosis in all simulations, including the setup where the values are extracted directly from Yule trees. The values for seem more similar to values from the Yule tree. We should expect this as here we have shown an exact limit distribution.
An overall assessment of the quantiles is given by the root–mean–square–error (RMSE) row in Tab. 1. We consider the quantiles at the , , , , , , , , , levels. The RMSE is defined as
[TABLE]
with dummy levels and . The normalizes the whole mean–square–error. We only look at the error at the tails, so we correct by the fraction of the distributions’ support that we consider. As a proxy for the true quantiles we take the pooled values (as explained in Tab. 1) from the “Yule columns”. The index runs over the repeats of the simulations.
The RMSE, when using , seems to be on the level of the RMSE of the “direct simulations”. has an error of about twice the size (both simulation methods). Looking at one can see that the RMSE is exactly on the level of the “Yule column’s” RMSE. This is even though we used a recursion of level , while an exact match of distributions should take place in the limit (infinite depth recursion). However, the rapid, exponential convergence of the contraction seems to make any differences invisible, already at this recursion level.
[FIGURE:]
4.2 Power of the tests
For a given test statistic to be useful one also needs to know its power, the ability to reject the null hypothesis (here Yule tree) when a given alternative one is true. For example, balance indices based only on topology like Sackin’s, Colless’ or cannot be expected to differentiate between any trees that are generated by different constant rate birth–death processes or by the coalescent. The rationale behind this is that the topologies induced by the contemporary species (i.e. we forget about lineages leading to extinct ones) are stochastically indistinguishable no matter what the death or birth rate is (Thm. 2.3, Cor. 2.4 Gernhard, 2008). Similarly, regarding the coalescent at the bottom of their p. Steel and McKenzie (2001) write “, one has the coalescent model [1,18,19]. In this model one starts with objects, then picks two at random to coalesce, giving objects. This process is repeated until there is only a single object left. If this process is reversed, starting with one object to give objects, then it is equivalent to the Yule model. Note that in the coalescent model there is commonly a probability distribution for the times of coalescences, but in the Yule model we ignore this element.” To differentiate between such trees one needs to take into consideration the branch lengths. Here we compare the power of the Sackin’s, Colless’, and indices at the significance level.
The null hypothesis is always that the tree is generated by a pure–birth process with rate . The alternative ones are birth–death processes , death rate , and using the TreeSim package, coalescent process (ape’s rcoal() function Paradis et al., 2004) and the biased speciation model for
[TABLE]
We also simulate a pure–birth process to check if the significance level is met. All trees were simulated with an root edge.
The so–called biased speciation model with parameter is the tree growth model as described by Blum and François (2005). In their words, “Assume that the speciation rate of a specific lineage is equal to . When a species with speciation rate splits, one of its descendant species is given the rate and the other is given the speciation rate where is fixed for the entire tree. These rates are effective until the daughter species themselves speciate. Values of close to [math] or yield very imbalanced trees while values around lead to over–balanced phylogenies.” We simulated such trees with in–house R code.
The quantiles of Sackin’s and Colless’ indices were obtained using Alg. 3. It is known (Eqs. 2 and 3 Blum and François, 2005; Blum et al., 2006) that after normalization (centring by expectation and dividing by ) in the limit they satisfy a contraction–type distribution of the form of Eq. (14), i.e.
[TABLE]
for . The function takes the form
[TABLE]
in Sackin’s case and
[TABLE]
in Colless’ case. It particular, studying the limit of Sackin’s index is equivalent to studying the Quicksort distribution (Blum and François, 2005). We can immediately see the main qualitative difference, the limit of the normalized cophenetic index has the square in , in the “recursion part” while Sackin’s and Colless’ have , .
Using repeats of Alg. 3 with recursion depth we obtained the following sets of quantiles , , , and , , respectively for the normalized Sackin’s and Colless’ indices.
Under each model we simulated trees conditioned on contemporary tips. We then checked if the tree was outside the “Yule zone” (Yang et al., 2017) by the procedure described in Alg. 1. We calculated the normalized Sackin’s, Colless’, discrete and continuous cophenetic indices (normalizations from Thm. 2.12). The functions sackin.phylo() and colless.phylo() of the phyloTop (Kendall et al., 2016) R package were used while the two cophenetic indices were calculated using a linear time in–house R implementation based on traversing the tree (Cor. 3 Mir et al., 2013). Two tests were considered, a two–sided one and a right–tailed one. For the discrete cophenetic index the quantiles from the simulation by Alg. 3 were considered, for the continuous those from (Tab. 1). The power is then estimated as the fraction of times the null hypothesis was rejected and represented in Tab. 3 by the corresponding Type II error rates. For the Yule tree simulation we can see that the significance level is met. All simulated trees are independent of the trees used to obtain the values in Tab. 1 and quantiles of Sackin’s and Colless’ indices. Hence, they offer a validation of the rejection regions. We summarize the power study in Tab. 3.
[FIGURE:]
As indicated in Alg. 1 one should first “correct” the tree for the speciation rate, when using the cophenetic index with branch lengths. The distributional results derived here on the cophenetic indices are for a unit speciation rate Yule tree. For a mathematical perspective this is not a significant restriction. If one has a pure–birth tree generated by a process with speciation rate , then multiplying all branch lengths by will make the tree equivalent to one with unit rate. Hence, all the results presented here are general up to a multiplicative constant. However, from an applied perspective the situation cannot be treated so lightly. For example, if we used the cophenetic index with branch lengths from a Yule tree with a very large speciation rate, then we would expect a significant deviation. However, unless one is interested in deviations from the unit speciation rate Yule tree, this would not be useful. Hence, one needs to correct for this effect. If the tree did come from a Yule process, then an estimate, , of the speciation rate by maximum likelihood can be obtained. For example, in the work here we used ape’s yule() function. Then, one multiplies all branch lengths in the tree by and calculates the cophenetic index for this transformed tree. It is important to point out that is only an estimate and hence a random variable. The effects of the this source of randomness on the limit distribution deserve a separate, detailed study. Balance indices that do not use branch lengths do not suffer from this issue but on the other hand miss another aspect of the tree—proportions between branch lengths that are non–Yule like.
The power analysis presented in Tab. 3 generally agrees with intuition and the power analysis done by Blum and François (2005). The first row shows that for all tests and statistics the significance level is approximately kept. Then, in the next three rows (coalescent and birth–death process) all topology based indices fail completely (the power is at the significance level). This is completely unsurprising as the after one removes all speciation events (with lineages) leading to extinct species from a birth–death tree, the remaining tree is topologically equivalent to a pure birth tree. The same is true for the coalescent model, its topology is identical in law to the Yule tree’s one. The cophenetic index with branch lengths has a high Type II error rate but is still better, than the topological indices. However, when one “corrects for ” this index manages to nearly ( trees were not rejected by the two–sided test) perfectly reject the coalescent model trees.
Power for the biased speciation model follows the same pattern as Blum and François (2005) observed. When imbalance is evident, , all ( uncorrected) tests were nearly perfect (two–sided discrete cophenetic is an exception). However, the correction significantly worsened the ability of to detect deviations. As imbalance decreased so did the power of the topological indices. For overbalanced trees one–sided tests failed, two sided worked (just as Blum and François, 2005, observed). The cophenetic index with branch lengths (without correction), that does not consider only the topology, was able to successfully reject the pure–birth tree for all (with only minimal Type II error for in the two–sided test case). Interestingly, ’s (both corrected and uncorrected) power seems invariant with respect to . These results are especially promising as seems to be an index that functions significantly better in the difficult, , regime, even after correcting.
At this stage we can point out that a normal approximation to the cophenetic indices’ limit distribution is not appropriate. When doing the above power study we observed that when using the quantiles of the standard normal distribution the right–tailed test based on rejects of Yule trees, based on rejects of Yule trees, two–sided test rejects of Yule trees and two–sided test rejects of Yule trees. The Type I error rates of the two–sided tests are within the observed Monte Carlo errors (in Tab. 3) but the right–tailed tests’ Type I error are evidently inflated. This confirms that the right tail of the scaled cophenetic index is much heavier than normal.
In short the power study indicates that the cophenetic index with branch lengths should be considered as an option to detect deviations from the Yule tree. This is because it is able to use information from two sources—the topology and time (a needed direction of development, as indicated in Ch. of Felsenstein, 2004). Actually, this is evident in the decomposition of Eq. (3). The s describe the topology and the s branch lengths. With more information a more powerful testing procedure is possible. Deviations that are not topologically visible, e.g. biased speciation in the regimes, are now detectable. To use one should correct for the effects of the speciation rate, as otherwise one merely detects deviations from the unit rate Yule tree. This correction is a mixed blessing. It can help or hinder detection.
4.3 Examples with empirical phylogenies
It is naturally interesting to ask how do the indices behave for phylogenies estimated from sequence data. Comparing a database of phylogenies, like TreeBase (http://www.treebase.org), with yet another index’s distribution under the Yule model should not be expected to yield interesting results. The Yule model has been indicated as inadequate to describe the collection of TreeBase’s trees (e.g. Blum and François, 2007). Therefore, we choose a particular study that estimated a tree and also reported a collection of posterior trees. Sosa et al. (2016a) is a recent work, providing all trees from BEAST’s (Drummond and Rambaut, 2007) output, well suited for such a purpose. Sosa et al. (2016a) estimate the evolutionary relationships between a set of tree ferns species. They report a posterior set of phylogenies (Sosa et al., 2016b).
In Tab. 4 we look what percentage of the trees from the posterior was accepted as being consistent with the Yule tree by the various tests and indices. It can be seen that the discrete cophenetic index has a high acceptance rate. The continuous one, which also takes into account branch lengths did not accept a single tree. However, this is lost when one corrects for the speciation rate (first ape’s multi2di() was used to make the trees binary ones). Most tests and indices rejected the Yule tree for the maximum likelihood estimate of the phylogeny with some exceptions. The two–sided discrete cophenetic index test did not reject the null hypothesis of the pure–birth tree. Also after correcting for the speciation rate (estimated at ), neither test based on the continuous cophenetic index rejected the Yule tree. Therefore, one should conclude (based on the “topological balances”) that the Yule tree null hypothesis can be rejected for this clade of plants.
We also followed Blum and François (2005) in looking at Yusim et al. (2001)’s phylogeny of the human immunodeficiency virus type 1 (HIV–1) group M gene sequences, available in the ape R package. The phylogeny consists of tips and Blum and François (2005) could not reject the null hypothesis of the pure–birth tree (using Sackin’s index amongst others). After pruning the tree to keep “only the old internal branches that corresponded to the oldest ancestors” they were able to reject the Yule tree. They conclude that the “results probably indicate a change in the evolutionary rate during the evolution which had more impact on cladogenesis during the early expansion of the virus.” Repeating their experiment we find that only the two versions of the cophenetic index point to a deviation but only in the two–sided test (see Tab. 5). Based only on the ’s test and that it conflicts with the conclusions of Sackin’s and Colless’ one should not draw any conclusions. However, as ’s test indicates a deviation, we can be inclined to reject the null hypothesis of the Yule tree. This is further strengthened by the fact that the significance remains after the correction for . Even though the topology as a whole seems consistent with the pure–birth tree the branch lengths are not. The fact that only the two–sided test rejected the Yule tree indicates that the HIV phylogeny is over–balanced in comparison to a pure–birth tree. In fact, in the biased speciation model tree over–balance is observed for values of close to (Blum and François, 2005). Such trees have a declining speciation rate as they grow and hence this supports Blum and François (2005)’s aforementioned explanation.
5 Almost sure behaviour of the cophenetic index
We study the asymptotic distributional properties of for the pure–birth tree model using techniques from our previous papers on branching Brownian and Ornstein–Uhlenbeck processes (Bartoszek, 2014; Bartoszek and Sagitov, 2015a, b; Sagitov and Bartoszek, 2012). We assume that the speciation rate of the tree is . The key property we will use is that in the pure–birth tree case the time between two speciation events, and (the first speciation event is at the root), is distributed, as the minimum of random variables. We furthermore, assume that the tree starts with a single species (the origin) that lives for time and then splits (the root of the tree) into two species. We consider a conditioned on contemporary species tree. This conditioning translates into stopping the tree process just before the speciation event, i.e. the last interspeciation time is distributed. We introduce the notation that is the height of the tree, is the time to coalescent of two randomly selected tip species and is the time between speciation events and (see Fig. 3 and Bartoszek and Sagitov, 2015b; Sagitov and Bartoszek, 2012).
Theorem 5.1
The cophenetic index is an increasing sequence of random variables, and has the recursive representation
[TABLE]
where is an indicator random variable whether tip split at the –th speciation event.
Proof From the definition we can see that
[TABLE]
where is the time to coalescent of tip species and . We now develop a recursive representation for the cophenetic index. First notice that when a new speciation occurs all coalescent times are extended by , i.e.
[TABLE]
where the “lone” is the time to coalescent of the two descendants of the split tip. The vector consists of [math]s and exactly one (a categorical distribution with categories all with equal probability). For each the marginal probability that is is . We rewrite
[TABLE]
and then obtain the recursive form
[TABLE]
Obviously, .
Proofof Theorem 2.4. Obviously
[TABLE]
and
[TABLE]
Hence, is a positive submartingale with respect to . Notice that
[TABLE]
Then, using the general formula for the moments of (Appendix A, Bartoszek and Sagitov, 2015b), we see that
[TABLE]
Hence, and are and by the martingale convergence theorem converges almost surely and in to a finite first and second moment random variable.
Corollary 5.2
* has finite third moment and is convergent.*
Proof We first recall the is positive. Using the general formula for the moments of again we see
[TABLE]
This implies that and hence convergence and finiteness of the third moment.
Remark 5.3
Notice that we (Appendix A, Bartoszek and Sagitov, 2015b) made a typo in the general formula for the cross moment of
[TABLE]
The should not be there, it will cancel with the from the derivative of the Laplace transform.
Proofof Theorem 2.7. We write as
[TABLE]
where are i.i.d. random variables.
Remark 5.4
We notice that we may equivalently rewrite
[TABLE]
The above and Eq. (3) are very elegant representations of the cophenetic index with branch lengths. They explicitly describe the way the cophenetic index is constructed from a given tree.
Proofof Theorem 2.11. The argumentation is analogous to the proof of Thm. 2.4 by using the recursion
[TABLE]
where is the number of nodes on the path from the root (or appropriately origin) of the tree to tip , (see also Bartoszek, 2014, esp. Fig. A.). An alternative proof for almost sure convergence can be found in Section 7.1.
6 Second order properties
In this Section we prove a series of rather technical Lemmata and Theorems concerning the second order properties of , and . Even though we will not obtain any weak limit, the derived properties do give insight on the delicate behaviour of and also show that no “simple” limit, e.g. Eq. (4), is possible. To obtain our results we used Mathematica for Linux x (–bit) running on Ubuntu LTS to evaluate the required sums in closed forms. The Mathematica code is available as an appendix to this paper.
Lemma 6.1
[TABLE]
Proof
[TABLE]
The following lemma is an obvious consequence of the definition of .
Lemma 6.2
For
[TABLE]
Lemma 6.3
[TABLE]
Proof Obviously
[TABLE]
We notice (as Bartoszek and Sagitov, 2015b; Bartoszek, 2016, in Lemmata and respectively) that we may write
[TABLE]
where , are two independent copies of , i.e. we sample a pair of tips twice and ask if both pairs coalesced at the –th speciation event. There are three possibilities, we (i) drew the same pair, (ii) drew two pairs sharing a single node or (iii) drew two disjoint pairs. Event (i) occurs with probability , (ii) with probability and (iii) with probability . As a check notice that . In case (i) , hence writing informally
[TABLE]
To calculate cases (ii) and (iii) we visualize the situation in Fig. 4 and recall the proof of Bartoszek and Sagitov (2015b)’s Lemma . Using Mathematica we obtain
[TABLE]
Similarly for case (iii)
[TABLE]
We now put this together as
[TABLE]
and we obtain (through Mathematica)
[TABLE]
Lemma 6.4
For
[TABLE]
Proof Obviously
[TABLE]
We notice that
[TABLE]
where , are the indicator variables if two independently sampled pairs coalesced at speciation events respectively. There are now two possibilities represented in Fig. 5 (notice that since the counterpart of event (i) in Fig. 4 cannot take place). Event (ii) occurs with probability and (iii) with probability . Event (iii) can be divided into three “subevents”.
Again we recall the proof of Bartoszek and Sagitov (2015b)’s Lemma and we write informally for (ii) using Mathematica
[TABLE]
In the same way for the subcases of (iii)
[TABLE]
We now put this together as
[TABLE]
and we obtain
[TABLE]
Theorem 6.5
[TABLE]
Proof We immediately have
[TABLE]
Theorem 6.6
[TABLE]
Proof We immediately may write using Lemmata 6.3, 6.4 and Mathematica
[TABLE]
Theorem 6.7
For we have
[TABLE]
Proof Again using Lemmata 6.3, 6.4, Mathematica and the fact that
[TABLE]
Theorem 6.8
[TABLE]
Proof We use Mathematica to first calculate
[TABLE]
For the second we again use Mathematica and the fact that the s are i.i.d. .
[TABLE]
For the third equality we use Mathematica and the fact that for independent families and of random variables we have
[TABLE]
As the s are i.i.d. we use Mathematica to obtain
[TABLE]
For the fourth equality we use the same properties and pair–wise independence of s.
[TABLE]
It is worth noting that the above Lemmata and Theorems were confirmed by numerical evaluations of the formulae and comparing these to simulations performed to obtain Fig. 1. As a check also notice that, as implied by variance properties,
[TABLE]
Theorem 6.9
[TABLE]
Proof Using the limit for the variance of (Thm. 6.6) and the independence of the s we have
[TABLE]
Now from Thms. 6.6 and 6.5 we have
[TABLE]
Plugging this in (and using Mathematica)
[TABLE]
Remark 6.10
Simulations presented in Fig. 6 and Thm. 6.9 suggest a different possible CLT, namely
[TABLE]
We sum over as and for all . It would be tempting to take the distribution to be a normal one. However, we should be wary after Rem. 2.9 and Fig. 1 that for our rather delicate problem even very fine simulations can indicate incorrect weak limits. It remains to study the variance of the conditional variance in Eq. (26). It is not entirely clear if this variance of the conditional variance will converge to [math]. Hence, it remains an open problem to investigate the conjecture of Eq. (27).
7 Alternative descriptions
7.1 Difference process
Let us consider in detail the families of random variables and . Obviously is times the number of pairs that coalesced after the speciation event for a given Yule tree. Denote
[TABLE]
As going from to means a new speciation event and coalescent at this new th event, then
[TABLE]
We also know by previous calculations that
[TABLE]
Let denote the number of newly introduced coalescent events after the –one when we go from to species. Obviously . Then, we may write
[TABLE]
Now,
[TABLE]
Therefore, for every , almost surely as it is a positive random variable whose expectation goes to [math]. However, is bounded by , as it can be understood in terms of the conditional (on tree) cumulative distribution function for the random variable —at which speciation event did a random pair of tips coalesce, i.e. for all
[TABLE]
Therefore, as is bounded by and the difference process
[TABLE]
goes almost surely to [math] we may conclude that converges almost surely to some random variable . In particular, this implies the almost sure convergence of to a limiting random variable . Furthermore, as and are both we may conclude that also converges almost surely. This means that the discrete version (all , corresponding to ) of the cophenetic index converges almost surely (compare with Thm. 2.11).
7.2 Polyá urn description
The cophenetic index both in the discrete and continuous version has the following Polyá urn description. We start with an urn filled with balls. Each ball has a number painted on it, [math] initially. At each step we remove a pair of balls, say with numbers and and return a ball with the number painted on it. We stop when there is only one ball, it will have value . Denote as the value painted on the –th ball in the –th step when we initially started with balls. Then we can represent the cophenetic index as
[TABLE]
Acknowledgments
I was supported by the Knut and Alice Wallenberg Foundation and am now by the Swedish Research Council (Vetenskapsrådet) grant no. –. I am grateful to the Barcelona Graduate School of Mathematics (BGSMath) for sponsoring the Workshop on Algebraical and Combinatorial Phylogenetics which significantly contributed to the development of my work. I would like to thank the whole Computational Biology and Bioinformatics Research Group of the Balearic Islands University for hosting me on multiple occasions, many discussions and suggestions on phylogenetic indices. My visits to the Balearic Islands University were partially supported by the the G S Magnuson Foundation of the Royal Swedish Academy of Sciences (grants no. MG–, MG–) and The Foundation for Scientific Research and Education in Mathematics (SVeFUM). I would like to acknowledge Gabriel Yedid for numerous discussions on the distribution of the cophenetic index and sharing his cophenetic index simulation R code. I am grateful to Cecilia Holmgren and Svante Janson for pointing me to the works on contraction–type distributions and many discussions. I would furthermore like to acknowledge Wojciech Bartoszek, Sergey Bobkov, Joachim Domsta, Serik Sagitov, Mike Steel for helpful comments and discussions related to this work. I am indebted to two anonymous reviewers, an anonymous editor and Haochi Kiang for careful reading of an earlier version of the manuscript and comments significantly improving it.
Appendix A: Mathematica code for Section 6
Appendix B: Counterparts of Rösler (1991)’s Prop. for the cophenetic index
Lemma 7.1
Define for
[TABLE]
and for , then
[TABLE]
Proof Writing out
[TABLE]
Therefore, assuming that
[TABLE]
If , we notice that and directly obtain
[TABLE]
Lemma 7.2
Define for ,
[TABLE]
and for ,
[TABLE]
then
[TABLE]
where is a positive random variable that converges to [math] almost surely with expectation decaying as and second moment as .
Proof Similarly, as in the proof of Lemma 7.1 we write out
[TABLE]
We denote and notice that it converges almost surely to [math] with . Now, assuming that
[TABLE]
If , we notice that and directly obtain
[TABLE]
Therefore, if we now denote
[TABLE]
we obtain the statement of the Lemma.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Agapow and Purvis (2002) P.–M. Agapow and A. Purvis. Power of eight tree shape statistics to detect nonrandom diversification: a comparison by simulation of two models of cladogenesis. Syst. Biol. , 51(6):866–872, 2002.
- 2Bartoszek (2014) K. Bartoszek. Quantifying the effects of anagenetic and cladogenetic evolution. Math. Biosci. , 254:42–57, 2014.
- 3Bartoszek (2016) K. Bartoszek. A central limit theorem for punctuated equilibrium. Ar Xiv e-prints , 2016.
- 4Bartoszek and Sagitov (2015 a) K. Bartoszek and S. Sagitov. A consistent estimator of the evolutionary rate. J. Theor. Biol. , 371:69–78, 2015 a.
- 5Bartoszek and Sagitov (2015 b) K. Bartoszek and S. Sagitov. Phylogenetic confidence intervals for the optimal trait value. J. App. Prob. , 52:1115–1132, 2015 b.
- 6Blum and François (2005) M. G. B. Blum and O. François. On statistical tests of phylogenetic tree imbalance: The Sackin and other indices revisited. Math. Biosci. , 195:141–153, 2005.
- 7Blum and François (2007) M. G. B. Blum and O. François. Which random processes describe the Tree of Life? A large–scale study of phylogenetic tree imbalance. Syst. Biol , 55(4):685–691, 2007.
- 8Blum et al. (2006) M. G. B. Blum, O. François, and S. Janson. The mean, variance and limiting distribution of two statistics sensitive to phylogenetic tree balance. Ann. Appl. Probab. , 16(4):2195–2214, 2006.