Convergence of maximum likelihood supertree reconstruction

1.   Introduction

High-throughput sequencing is making large collections of sequences available to researchers at a low cost. These genomic data represent a broad spectrum of life and motivates studies of the problem of reconstructing large phylogenetic trees using statistical methods. Those data, however, also come from different sources, cover different genomic regions on which the evolutionary processes happen very differently, and may not be collected on the same set of species. Thus, combining trees on different, overlapping sets of species into a "supertree" has become a popular approach for reconstructing large species trees. Over the years, several methods of supertree reconstruction have been developed ^[8], and analyses of supertrees continue to play more and more important roles in the search for answers of many fundamental evolutionary questions.

While supertree methods have been of great interest in phylogenetics, little is known about their theoretical properties, especially in non-asymptotic settings. In cases when the individual trees are gene trees and the set of taxa are the same across the input trees, several statistical consistent methods have been derived ^{[5,7,11,14,16,17,18,20,29]}. However, most proofs of statistical consistency have been analyzed under the condition that each of the gene trees can be estimated accurately, and do not provide a guarantee of robustness to gene tree estimation errors ^[25]. This is a critical concern because when the species tree is constructed using sequences taken from large genomic regions, those regions have a high chance of involving some recombination, which violates one of the main assumptions of the multi-species coalescent model. On the other hand, limiting analyses to short regions increase gene tree estimation error (GTEE) and various summary methods had impaired accuracy when the error was high ^[10,19]. It is later proved in ^[24] that when the sequence length of each locus is bounded and the gene tree cannot be estimated reliably, most summary methods that estimate the species tree by combining gene trees are not statistically consistent.

Because some coalescent-based summary methods sometimes produce less accurate estimates than concatenation ^[3,19,22], seemingly as a result of GTEE, the question of whether provable guarantees can be established in the presence of GTEE (for both species trees and supertree) naturally arises. In the context of species tree estimation, ^[25] establish statistical consistency of the Rooted Triplet Consensus method and the Maximum Pseudo-likelihood for Estimating Species Trees method ^[18] under GTEE and provide bounds on the sampling complexity for these methods to construct the correct species tree with high probability.

The maximum likelihood (ML) supertree method is proposed by ^[28] based on a probability model that permits "errors" in gene tree topologies and allows the species tree to be estimated even if there is topological conflict amongst gene trees. ^[28] shows that ML estimate of the species tree is topologically consistent under fairly general conditions and also shows that the method of Matrix Representation with Parsimony ^[2] may be inconsistent under these same conditions. However, while the ML estimate is topologically consistent, no results on the convergence rate of the estimator are obtained. In this paper, we propose a new analytic approach to study the convergence of the ML supertree method. Based on embedding trees into a metric space, we establish the conditions for which the convergence rate of the ML supertree can be obtained. We verify these conditions for the popular exponential error model of gene trees, thus obtain a polynomial convergence rate for the ML supertree to the true species tree under this model.

2.   Mathematical framework

In this paper, the term phylogenetic tree refers to a tree $\mathbb{T}$ with leaves labeled by a set of species. Each branch of $\mathbb{T}$ is associated with a non-negative branch length. A tree is said to be resolved if it is bifurcating and all branch lengths are positive. Given an unrooted phylogenetic tree $\mathbb{T}$ on a finite set $\mathscr{L}$ of species, any subset $\mathscr{L}'$ of $\mathscr{L}$ induces a phylogenetic tree on $\mathscr{L}'$, denoted $\mathbb{T} |_{ \mathscr{L}'}$, which is the subtree of $\mathbb{T}$ that connects the species in $\mathscr{L}'$ only.

We regard the tree space as a metric space $(\mathcal{T}, d)$ where $\mathcal{T}$ is the set of all phylogenetic trees with branch lengths bounded from above by a positive constant $g_0$ and $d$ is a continuous metric. For simplicity of presentation, we will assume that $d$ is the BHV distance on the set of trees with $n$ species ^[4], but note that the analysis of the paper can be extended to any continuous and locally-Euclidean distance, including the branch-score distance ^[15] and the AGPS distance ^[1].

To describe trees that are "near" to each other, the BHV distance use the class of nearest neighbor interchange (NNI) moves ^[23]. An NNI move is defined as a transformation that collapses an interior branch to zero and then expands the resulting degree 4 vertex into a branch in a different way. The BHV space models the set of trees $\mathcal{T}$ on $n$ species as a cubical complex consisting of a collection of orthants, each isomorphic to $\mathbb{R}_{\ge 0}^{2n-3}$. Each orthant of $\mathcal{T}$ corresponds uniquely to a tree topology, and the coordinates in each orthant parameterize the branch lengths for the corresponding tree. The adjacent orthants of the complex with the same dimension correspond to NNI-adjacent trees.

The BHV space is equipped with a natural metric distance: the shortest path lying in the BHV space between the points. If two points lie in the same orthant, this distance is the usual Euclidean distance. If two points are in different orthants, they can be joined by a sequence of straight segments, with each segment lying in a single orthant. We can then measure the length of the path by adding up the lengths of the segments. The distance between the two trees $\mathbb{T}_1$ and $\mathbb{T}_2$ on the BHV space is defined as the minimum of the lengths of such segmented paths joining the two points.

Throughout the paper, we assume that the evolution of the species has not involved reticulate processes, and there exists an underlying "true species tree" in $\mathcal{T}$, denoted by $\mathbb{T}^*$. Furthermore, $\mathbb{T}^*$ is a resolved tree with leaf set $\mathscr{L}^*$. In a supertree reconstruction problem, we observe a sequence of gene trees $(\mathbb{T}_1, \mathbb{T}_2, \ldots, \mathbb{T}_k)$, where $\mathbb{T}_i$ has leaf set $\mathscr{L}_i$, and wish to combine these trees into a phylogenetic tree $\hat {\mathbb{T}}_k$ on the union of the leaf sets

In this paper, we consider a tree-generating probability model $\left(\mathcal{P}_{ \mathbb{T}} \right)_{ \mathbb{T} \in \mathcal{T}}$ such that for each set of species $\mathscr{L} \subset \mathscr{L}^*$, $\mathcal{P}^{ \mathscr{L}}_{ \mathbb{T}}(\cdot)$ is a distribution of trees forming by species in $\mathscr{L}$. We assume that the observed gene trees $\{ \mathbb{T}_i\}_{i = 1}^k$ are independently distributed according to $\left \{ \mathcal{P}^{ \mathscr{L}_i}_{ \mathbb{T}^*}(\cdot) \right \}_{i = 1}^k$ respectively. In other words, the joint probability of the observed gene trees is

Remark 2.1. We note that the tree-generating probability model $\mathcal{P}_{ \mathbb{T}^*}$ may not necessarily be nested. That is, given two nested sets of leaves $\mathscr{L}_1 \subset \mathscr{L}_2 \subset \mathscr{L}^*$, although the probability $\mathcal{P}^{ \mathscr{L}_2}_{ \mathbb{T}^*}$ (which is used to generate trees with leaf set $\mathscr{L}_2$) also induces a natural distribution on the set of trees with leaf set $\mathscr{L}_1$, this probability $\mathcal{P}^{ \mathscr{L}_2}_{ \mathbb{T}^*} |_{ \mathscr{L}_1}$ may not be the same as $\mathcal{P}^{ \mathscr{L}_1}_{ \mathbb{T}^*}$.

Given the observed gene trees $(\mathbb{T}_1, \mathbb{T}_2, \ldots, \mathbb{T}_k)$, the ML supertree is defined as

where $\ell_k(\mathbb{T})$ denotes the log-likelihood function

To enable theoretical analyses of the ML supertree method, we make the following assumptions.

Assumption 2.1 (Weak covering property). The sequence of subsets of $\mathscr{L}^*$ satisfies the $\text{weak covering property}$: there exists $c_1 > 0, \gamma > 1/2$, and $K > 0$ such that for each subset $\mathscr{L}$ of taxa from $\mathscr{L}^*$ of size $4$,

Here, $| \mathcal{A}|$ denotes the number of elements in the set $\mathcal{A}$.

We note that the weak covering property ensures that as the number of gene trees increases, all quartets (subtrees with $4$ leaves) of $\mathbb{T}^*$ are visited with enough frequency to enable a reliable estimate for each of the quartets. Assumption $2.1$ is a direct generalization of the covering property introduced in ^[28], which can be obtained from Assumption $2.1$ by setting $\gamma = 1$.

Assumption 2.2 (Model identifiability). For all $\mathscr{L} \subset \mathscr{L}^*$, the distribution $\mathcal{P}^{ \mathscr{L}}_{ \mathbb{T}^*}(\cdot)$ is identifiable. That is, if $d(\mathbb{T}^* |_{ \mathscr{L}}, \mathbb{T} |_{ \mathscr{L}}) > 0$, then $\mathit{\text{KL}}(\mathcal{P}^{ \mathscr{L}}_{ \mathbb{T}^*}, \mathcal{P}^{ \mathscr{L}}_{ \mathbb{T}}) > 0$.

Assumption $2.2$ guarantees that it is at least possible to reconstruct the restriction $\mathbb{T}^*|_{ \mathscr{L}}$ of $\mathbb{T}^*$ to the subset $\mathscr{L}$ with a complete knowledge of $\mathcal{P}^{ \mathscr{L}}_{ \mathbb{T}^*}$. This assumption is similar to, but distinct from, the condition of basic centrality introduced by ^[28], which requires that for all subsets of $\mathscr{L} \subset \mathscr{L}^*$ of size 4,

for all trees $\mathbb{T}'$ on leaf set $\mathscr{L}$ that are different from $\mathbb{T}^* |_ \mathscr{L}$, and where $\eta > 0$:

● On one hand, the basic centrality condition implies that if $d(\mathbb{T}^* |_{ \mathscr{L}}, \mathbb{T} |_{ \mathscr{L}}) > 0$, the modes of the distributions $\mathcal{P}^{ \mathscr{L}}_{ \mathbb{T}^*}$ and $\mathcal{P}^{ \mathscr{L}}_{ \mathbb{T}}$ are different. From this, we can deduce that $\text{KL}(\mathcal{P}^{ \mathscr{L}}_{ \mathbb{T}^*}, \mathcal{P}^{ \mathscr{L}}_{ \mathbb{T}}) > 0$. Thus, our identifiability assumption is somewhat weaker than this condition. Assumption $2.2$ also does not impose any assumption on the family of distribution themselves and thus can be applied to a wider class of probabilistic models.

● On the other hand, the basic centrality condition only concerns leaf sets of size 4, while Assumption $2.2$ impose restrictions on leaf sets of all sizes. However, we note that conditions on subset of leaves (such as the basic centrality condition) only work under the implicit assumption that there are some connection between $\mathcal{P}^{ \mathscr{L}_2}_{ \mathbb{T}}|_{ \mathscr{L}_1}$ and $\mathcal{P}^{ \mathscr{L}_1}_{ \mathbb{T}}$ for $\mathscr{L}_1 \subset \mathscr{L}_2$. Since our framework does not assume any nested structure in the probability model, Assumption $2.2$ is more appropriate.

Finally, we impose the following regularity conditions on the tree-generating probability model:

Assumption 2.3. (Regularity)

(a) For all $\mathbb{T} \in \mathcal{T}$, $\mathscr{L} \subset \mathscr{L}^*$, and any tree $\mathbb{T}'$ with leaf set $\mathscr{L}$,

$i.$ $\mathcal{P}^{ \mathscr{L}}_{ \mathbb{T}}(\mathbb{T}') > 0$,

$ii.$ $\log \mathcal{P}^{ \mathscr{L}}_{ \mathbb{T}}(\mathbb{T}')$ is a locally-Lipchitz function with respect to $\mathbb{T}$, and the Lipchitz constant does not depend on $\mathbb{T}'$.

(b) There exist $c_2, c_3 > 0$, $m \ge 2$ such that for any leaf set $\mathscr{L}$ and tree $\mathbb{T}$, if $d(\mathbb{T}^*, \mathbb{T}) \le c_2$, then

Remark 2.2. If for each tree $\mathbb{T}'$ with leaf set $\mathscr{L}$, the probability density function $\mathcal{P}^{ \mathscr{L}}_{ \mathbb{T}}(\mathbb{T}')$ is an analytic function with respect to $\mathbb{T}$ in a neighborhood of $\mathbb{T}^*$, then Assumption $2.3$(b) holds.

Proof. Note that the function $\mathbb{E}_{\mathcal{P}_{ \mathbb{T}^*}}\left [\log \mathcal{P}^{ \mathscr{L}}_{ \mathbb{T}}(\mathbb{T}')\right]$ (with respect to $\mathbb{T}$) is analytic in a neighborhood $\mathcal{U}$ of $\mathbb{T}^*$. Define

For all $\mathbb{T} \in \mathcal{A}$, we have

By Assumption 2.2, we conclude that $\mathbb{T} |_ \mathscr{L} = \mathbb{T}^* |_ \mathscr{L}$ for all $\mathbb{T} \in \mathcal{A}$. Applying Łojasiewicz inequality [13, Theorem 1], we deduce that there exists $c_3 > 0$ and $m \ge 2$ such that for any $\mathbb{T}$ in the neighborhood,

This implies the result.

Throughout this paper, we will assume that Assumptions $2.1$, $2.2$, and $2.3$ hold. The main point of Assumption $2.1$ is to guarantee that every subset of $4$ species is sufficiently represented by the data. It is worth noticing that $4$ is the smallest number of species that carries the information of the topology of an unrooted tree. Therefore, if a set of $4$ species is only found in few trees, we will not be able to reconstruct their phylogenetic relationship. While Assumption $2.2$ is technical, its goal is to guarantee that there is a unique "true" supertree. Finally, Assumption $2.3$ ensures that the model has enough information to distinguish between the "true" supertree and any alternative with substantial accuracy.

In the next section, we will establish the following convergence rate of the ML supertree.

Theorem 1. Under Assumptions $2.1$, $2.2$, and $2.3$, for any $\delta > 0$, there exist constants $m > 0$, $C_{\delta} > 0$ and $K_\delta > 0$ such that for all $k \ge K_\delta$, we have

with probability at least $1-\delta$.

Theorem 1 establishes that under fairly mild regularity conditions, the ML supertree is consistent and has a polynomial convergence rate. Notably, the result holds for all identifiable family of locally-analytic distributions (Remark 2.2). The degree of this (polynomial) convergence rate depends on the sampling scheme of the leaves (characterized by the covering coefficient $\gamma$) and on the geometry of the probabilistic model (characterized by the constant $m$ in Assumption $2.3$(b)). We further note that Assumption $2.3$(b)) is only required for establishing the convergence rate of the ML estimator and the absence of this condition does not affect the proof of consistency.

Corollary 1. Under Assumptions $2.1$, $2.2$, and $2.3(a)$, the ML supertree is consistent.

Remark 2.3. Theorem $1$ is still valid if the universal Assumption $2.3$(b) is replaced by the following condition:

There exist $c_2, c_3 > 0$, $m \ge 2$ such that if $d(\mathbb{T}^*, \mathbb{T}) \le c_2$, then

Remark 2.4. We note that the bound of the convergence rate in Theorem $1$ is not sharp.

3.   Convergence of maximum likelihood supertree

To enable the analysis of convergence, we define

and

We refer to $R_{k}$ and $E_{\mathcal{P}_{ \mathbb{T}^*}}[R_{k}(\mathbb{T})]$ as the empirical risk function and the expected risk function, respectively.

An intuitive argument for the consistency of the ML estimator can be described as follows. As the number of gene trees increases, we have

with sufficiently high probability. Therefore, the ML estimator will converge to the optimal value of the risk function, which is attained at the true species tree $\mathbb{T}^*$. This simple argument is formalized by a lower bound of the expected risk function (Section 3.1) and a uniform concentration bound on the deviation of the empirical risk function and its expectation (Section 3.2).

3.1.   Lower bound of the expected risk

Lemma 3.1. There exist a neighborhood $\mathcal{U}$ of $\mathbb{T}^*$ and $c_4 > 0$ such that

Proof. We note that

Consider an arbitrary leaf set $\mathscr{L} \subset \mathscr{L}^*$ of size 4 and define $\mathcal{I} = \{i \le k: \mathscr{L} \subset \mathscr{L}_i\}$. By the weak covering property, we have $| \mathcal{I}| \ge c_1 k^{\gamma}$. Thus,

for all $\mathbb{T}$ in some neighborhood $U$ around $\mathbb{T}^*$. Here, the second inequality comes from Assumption $2.3$(b).

On the other hand, ^[9] (Lemma 6.2 (i)) proved that for some leaf set $\mathscr{L} \subset \mathscr{L}^*$ of size 4, we have

We deduce that

Lemma 3.2. For any leaf set $\mathscr{L}$ and any $x > 0$, there exists $C_x > 0$ such that

when $d(\mathbb{T}^*|_ \mathscr{L}, \mathbb{T}|_ \mathscr{L}) \ge x$.

Proof. Assume that there exists a sequence of tree $\{ \mathbb{T}_i\}$ such that $d(\mathbb{T}^*|_ \mathscr{L}, \mathbb{T}_i|_ \mathscr{L}) \ge x$ and $\text{KL}(\mathcal{P}^{ \mathscr{L}}_{ \mathbb{T}^*}, \mathcal{P}^{ \mathscr{L}}_{ \mathbb{T}_i}) \to 0$. Since the tree space is compact, we can extract a sub-sequence $\{ \mathbb{T}_{i_j}\}$ of $\{ \mathbb{T}_i\}$ that converges to some tree $\mathbb{T}_0$. We deduce that $\text{KL}(\mathcal{P}^{ \mathscr{L}}_{ \mathbb{T}^*}, \mathcal{P}^{ \mathscr{L}}_{ \mathbb{T}_0}) = 0$ and $\mathbb{T}^*|_ \mathscr{L} \ne \mathbb{T}_0|_ \mathscr{L}$. This contradicts Assumption $2.2$.

Lemma 3.3. Let $\mathcal{V}$ be a neighborhood of $\mathbb{T}^*$. There exist $c_{ \mathcal{V}} > 0$ such that for any $\mathbb{T} \not \in \mathcal{V}$, we have

Proof. By ^[9] (Lemma 6.2 (i)), there exists $\mathscr{L} \subset \mathscr{L}^*$ of size 4 such that

since $\mathbb{T} \not \in \mathcal{V}$. Define $\mathcal{I} = \{i \le k: \mathscr{L} \subset \mathscr{L}_i\}$, we note that for all $i \in \mathcal{I}$

which implies $\text{KL}(\mathcal{P}^{ \mathscr{L}_i}_{ \mathbb{T}^*}, \mathcal{P}^{ \mathscr{L}_i}_{ \mathbb{T}}) \ge C''_{ \mathcal{V}}$ for some $C''_{ \mathcal{V}} > 0$ by Lemma 3.2.

By the weak covering property, we have $| \mathcal{I}| \ge c_1 k^{\gamma}$. Thus,

which completes the proof.

3.2.   Uniform concentration bound

Lemma 3.4 (Concentration bound). For any $\delta > 0$, $k \ge 3$, there exists $c_5(\delta)$ such that

with probability at least $1-\delta$.

Proof. Since the functions $\log \mathcal{P}^{ \mathscr{L}_i}_{ \mathbb{T}}(\mathbb{T}_i)$ is locally Lipschitz with respect to $\mathbb{T}$ (Assumption $2.3$(a)) and number of species is finite, there exists $C_1 > 0$ such that

Since the tree space is compact, $C_1 d(\mathbb{T}, \mathbb{T}')$ is bounded by a constant $C_2$. Using Hoeffding's inequality ^[12], we obtain

On the other hand, we have $|R_{k}(\mathbb{T}) - R_k(\mathbb{T}')| \le C_1 d(\mathbb{T}, \mathbb{T}')$ and $|\mathbb{E}_{\mathcal{P}_{ \mathbb{T}^*}}\left [R_{k}(\mathbb{T}) \right] - \mathbb{E}_{\mathcal{P}_{ \mathbb{T}^*}}\left [R_{k}(\mathbb{T}') \right]| \leq C_1 d(\mathbb{T}, \mathbb{T}')$. Thus, if we define the events

and

then $\mathcal{B}(x, k, \mathbb{T}) \subset \mathcal{A}(x, k, \mathbb{T})$. Note that the total number of balls of radius $x/(4C_1\sqrt{k})$ required to cover the tree space is bounded above by

where $g_0$ is the upper bound of branch lengths and $C_n$ is a constant depending on the number of species. To obtain the desired inequality, we will chose $x$ such that

which can be done with $x = C(\delta, n, C_2, g_0) \log k$.

3.3.   Proof of Theorem 1

First, we establish that the ML estimator is consistent.

By Lemma $3.4$, we have

with probability at least $1-\delta$, since $\hat {\mathbb{T}}_k$ is the maximizer of the empirical risk function.

On the other hand, let $\mathcal{V}$ be a neighborhood of $\mathbb{T}^*$, by Lemma 3.3, we have

Since $\gamma > 1/2$, there exists $K_{\delta, \mathcal{V}}$ such that for $k \ge K_{\delta, \mathcal{V}}$, we have

We deduce that for $k \ge K_{\delta, \mathcal{V}}$, $\hat {\mathbb{T}}_k \in \mathcal{V}$ with probability at least $1 -\delta$. This proves that the ML estimator is consistent.

Next, we will derive the convergence rate of the ML supertree. From Lemma 3.1, there exist a neighborhood $\mathcal{U}$ of $\mathbb{T}^*$ and $c_4, m > 0$ such that

For $k \ge K_{\delta, \mathcal{U}}$, we have $\hat {\mathbb{T}}_k \in \mathcal{U}$ with probability at least $1 -\delta$. By Lemma 3.4, with probability $1 - 2\delta$, we obtain

Here,

because $\hat {\mathbb{T}}_k$ is the maximizer of $R_k$. This completes the proof.

4.   Applications: convergence under the exponential model

In this section, we consider the exponential model, a simple model of gene tree estimation errors in which the probability of observing a given tree decreases exponentially with its distance from the species tree. For a detailed discussion on its validity and applicability, we refer the readers to ^[28] and the references therein.

Suppose $d$ is some metric on the set of trees, in the exponential model, the probability of reconstructing any tree $\mathbb{T}'$ with a leaf set $\mathscr{L}$, when $\mathbb{T}^*$ is the generating tree is proportional to an exponentially decaying function of the distance from $\mathbb{T}'$ to $\mathbb{T}^*|_{ \mathscr{L}}$:

Here $\beta_{ \mathscr{L}}$ is a constant that depends only on the set of leaves $\mathscr{L}$, while $\alpha_{ \mathbb{T}^*, \mathscr{L}}$ is the normalizing constant to ensure that $\mathcal{P}^{ \mathscr{L}}_{ \mathbb{T}^*}$ is a density function.

In this section, we verify the identifiability and regularity conditions for the exponential model when $d$ is any continuous tree distance such that if $\mathbb{T}$ and $\mathbb{T}'$ have the same topology, then $d(\mathbb{T}, \mathbb{T}')$ is the Euclidean distance.

Theorem 2. Under the exponential model with Assumptions $2.1$, for any $\delta > 0$, $C_{\delta} > 0$ and $K_\delta > 0$ such that for all $k \ge K_\delta$,

with probability at least $1-\delta$, where

Proof. First, we note that for all $\mathbb{T}'$ with leaf set $\mathscr{L}$, $\mathcal{P}^{ \mathscr{L}}_{ \mathbb{T}}(\mathbb{T}')$ is a positive and locally Lipschitz function. To obtain an upper bound on the convergence rate for the ML supertree, we need to verify Assumption $2.2$ and Assumption $2.3(b)$/Remark $2.3$.

Identifiability. Since it is sufficient to prove either $\text{KL}(\mathcal{P}^{ \mathscr{L}}_{ \mathbb{T}^*}, \mathcal{P}^{ \mathscr{L}}_{ \mathbb{T}}) > 0$ or $\text{KL}(\mathcal{P}^{ \mathscr{L}}_{ \mathbb{T}}, \mathcal{P}^{ \mathscr{L}}_{ \mathbb{T}^*}) > 0$, we can assume that $\alpha_{ \mathbb{T}^*, \mathscr{L}} \geq \alpha_{ \mathbb{T}, \mathscr{L}}$.

We denote

and pick $c$ small enough such that $\beta_{ \mathscr{L}} d(\mathbb{T}', \mathbb{T} |_{ \mathscr{L}}) \ge 2c$. We note that in $A_c$,

Since $\alpha_{ \mathbb{T}^*, \mathscr{L}} \geq \alpha_{ \mathbb{T}, \mathscr{L}}$, we have

By Pinsker's inequality,

which establishes the identifiability of the exponential model. Here, $\mu$ is the Lebesgue measure.

Regularity. Consider a fixed tree $\mathbb{T} \in \mathcal{T}$ and leaf set $\mathscr{L}$, we first assume that $\alpha_{ \mathbb{T}^*, \mathscr{L}} \geq \alpha_{ \mathbb{T}, \mathscr{L}}$. If we define

then

Using the same argument as above, we have

It can be verified that there exists $C > 0$ such that for all $x \in [0, C]$, we have

Moreover since $d$ is Euclidean inside each orthant, if $d(\mathbb{T}, \mathbb{T}^*)$ is small enough, we have

Thus, let $\mathcal{W}$ be a neighborhood in the same topology of $\mathbb{T}^*$ such that

and both Eqs $(4.1)$ and $(4.2)$ are satisfied, we have

for some constant $C > 0$ independent of $\mathbb{T}$.

For the case when $\alpha_{ \mathbb{T}^*, \mathscr{L}} < \alpha_{ \mathbb{T}, \mathscr{L}}$ we define

and the argument proceeds similarly. This validates Assumption $2.3(b)$ for exponential model with $m = 4n-4$. However, if we use the regularity condition in Remark $2.3$, we can obtain the result with

5.   Application to yeast species-tree reconstruction

We apply our method to reconstruct the yeast species-tree from 106 genes of eight yeast species obtained from the $\texttt{R}$ package $\texttt{phangorn}$ ^[27]. Specifically, we reconstruct the gene-tree for each gene using UPGMA (unweighted pair group method with arithmetic mean) method implemented in $\texttt{phangorn}$. To replicate the setting of this paper, we randomly remove two species from each reconstructed gene-tree. In other words, the inputs are $106$ reconstructed gene-trees, each has $6$ species (see Figure 1 for examples).

We apply our method to estimate the full species-tree using the exponential model with BHV distance. The $\texttt{R}$ function $\texttt{dist.multiPhylo}$ in the package $\texttt{distory}$ ^[6] is used to compute the BHV distance. The reconstructed species-tree (see Figure 2) is consistent with the well-known result of ^[26]. This confirms the applicability of our method to reliably estimate supertrees.

6.   Discussion and conclusions

In this paper, we propose a novel analytic approach to analyze the convergence of the ML supertree method. Instead of focusing on reconstructing the correct discrete topology of the species tree as in previous studies [25,28, e.g.], we employ a continuous model of the tree space and analyze the ML estimator on this metric space, aiming at recovering both the topology and the branch lengths of the species tree. This framework enables us to use tools from statistical learning and information theory to establish the convergence rate of the ML estimator and at the same time, to weaken the conditions to obtain consistency and convergence of the estimator. Our weak covering property is an extension of the classical covering property ^[28] and provides a considerable relaxation on the sampling schemes for supertree estimation. Our identifiability condition is also more intuitive and generalizable than the well-known basic centrality condition and does not impose constraints on the shape of the probabilistic model of gene tree estimation errors. Our information-theoretical approach to analyze statistical estimator on tree spaces is of independent interest and can be extended to other problems in phylogenetics.

There are several avenues for future directions for this work. The first direction is extending our results to other practical models of phylogenetic errors, including the multiple-coalescent model (along with a detailed model of the effects of short sequence length on the accuracy in estimating the individual gene trees). Second, while our result provides a polynomial bound on the convergence rate, the power of the convergence (characterized by the geometric constant $m$) is not sharp. A sharper bound of the convergence rate would be of great interest to the field (from both theoretical and applied perspective) and would require further understanding of the tree-generating probabilistic model.

Acknowledgments

LSTH was supported by startup funds from Dalhousie University, the Canada Research Chairs program, the NSERC Discovery Grant RGPIN-2018-05447, and the NSERC Discovery Launch Supplement DGECR-2018-00181. VD was supported by a startup fund from University of Delaware and National Science Foundation grant DMS-1951474.

Conflict of interest

The author declares no conflict of interest in this paper.