On the convergence of the maximum likelihood estimator for the transition rate under a 2-state symmetric model

TL;DR

This paper proves that the maximum likelihood estimator for the transition rate in a 2-state symmetric trait evolution model converges to the true value as the number of species increases, under certain regularity conditions.

Contribution

It establishes the convergence and provides bounds for the MLE of the transition rate in a large-tree limit for a 2-state symmetric model, extending previous results to correlated trait data.

Findings

MLE converges to the true transition rate under regularity conditions

Provides an upper bound for the estimation error

Results apply to various practical tree models such as Yule and coalescent trees

Abstract

Maximum likelihood estimators are used extensively to estimate unknown parameters of stochastic trait evolution models on phylogenetic trees. Although the MLE has been proven to converge to the true value in the independent-sample case, we cannot appeal to this result because trait values of different species are correlated due to shared evolutionary history. In this paper, we consider a $2$-state symmetric model for a single binary trait and investigate the theoretical properties of the MLE for the transition rate in the large-tree limit. Here, the large-tree limit is a theoretical scenario where the number of taxa increases to infinity and we can observe the trait values for all species. Specifically, we prove that the MLE converges to the true value under some regularity conditions. These conditions ensure that the tree shape is not too irregular, and holds for many practical scenarios such as trees with bounded edges, trees generated from the Yule (pure birth) process, and trees generated from the coalescent point process. Our result also provides an upper bound for the distance between the MLE and the true value.

Figures

Equations

$${\mathbf{P}}_{\mu }(t)=\left(\begin{array}{cc}\frac{1}{2}+\frac{1}{2}{e}^{-2\mu t}& \frac{1}{2}-\frac{1}{2}{e}^{-2\mu t}\\ \frac{1}{2}-\frac{1}{2}{e}^{-2\mu t}& \frac{1}{2}+\frac{1}{2}{e}^{-2\mu t}\end{array}\right).$$

$$P_{\mathbb{T},\mu }(\mathbf{Y}):=\mathrm{I}\mathrm{P}(\mathbf{Y}|\mathbb{T},\mu )=\frac{1}{2}\sum _y\left(\prod _{(u,v)\in E}[{\mathbf{P}}_{\mu }(d_{uv}){]}_{y_u{y}_v}\right)$$

$$R_{\mathbb{T},\mu }(\mathbf{Y})={\ell }_{\mathbb{T},{\mu }^{\ast }}(\mathbf{Y})-{\ell }_{\mathbb{T},\mu }(\mathbf{Y}).$$

$$\text{KL}(P_{\mathbb{T},{\mu }^{\ast }}\Vert P_{\mathbb{T},\mu })=\mathbb{E}[{\ell }_{\mathbb{T},{\mu }^{\ast }}(\mathbf{Y})-{\ell }_{\mathbb{T},\mu }(\mathbf{Y})]=\mathbb{E}[R_{\mathbb{T},\mu }(\mathbf{Y})]$$

$${\ell }_{{\mathbb{T}}^{\prime },\mu }(\mathbf{Y})={\ell }_{\mathbb{T},\mu }(\mathbf{Y})$$

$$|{\ell }_{\mathbb{T},\mu }(\mathbf{Y})-{\ell }_{{\mathbb{T}}_1,\mu }({\mathbf{Y}}_1)-{\ell }_{{\mathbb{T}}_2,\mu }({\mathbf{Y}}_2)|\le \max \left\{\log 2,\log \frac{1}{1-e^{-2\mu d}}\right\}.$$

$$\left|\frac{1}{n}{\ell }_{\mathbb{T},{\mu }_1}(\mathbf{Y})-\frac{1}{n}{\ell }_{\mathbb{T},{\mu }_2}(\mathbf{Y})\right|\le C|{\mu }_1-{\mu }_2|\forall {\mu }_1,{\mu }_2\in [\underline{\mu },\overline{\mu }],\forall \mathbb{T}$$

$$\hat{\mu }=\underset{\mu \in [\underline{\mu },\overline{\mu }]}{\mathrm{argmax}}{\ell }_{\mathbb{T},\mu }(\mathbf{Y}).$$

$$f(x)=\max \left\{\log 2,\log \frac{1}{1-e^{-2x}}\right\}\text{ and }{S}_{\mathbb{T},\mu }=\sum _{e\in E(\mathbb{T})}f(\mu t_e{)}^2.$$

$$|\hat{\mu }-{\mu }^{\ast }|\le C_{\delta ,\alpha ,\underline{d},\overline{d},\underline{\mu },\overline{\mu }}\left(\frac{\sqrt{\log n}}{n^{(2-\alpha )/6}}\right)$$

$$\mathrm{I}\mathrm{P}({\mathcal{C}}_{ij}^2=y|{\mathcal{A}}_{ij}=0)=\mathrm{I}\mathrm{P}({\mathcal{C}}_{ij}^2=y|{\mathcal{A}}_{ij}=1)=\{\begin{array}{ll}\frac{1}{2}(1+e^{-2\mu d_{ij}}),& \text{mbox}ify=0\\ \frac{1}{2}(1-e^{-2\mu d_{ij}}),& \text{mbox}ify=1.\end{array}$$

$$\mathrm{I}\mathrm{P}({\mathcal{C}}_{i_1{j}_1}^2=y_1,{\mathcal{C}}_{i_2{j}_2}^2=y_2,\dots ,{\mathcal{C}}_{i_m{j}_m}^2=y_m)$$

$$=\mathbb{E}\left(\mathrm{I}\mathrm{P}({\mathcal{C}}_{i_1{j}_1}^2=y_1,{\mathcal{C}}_{i_2{j}_2}^2=y_2,\dots ,{\mathcal{C}}_{i_m{j}_m}^2=y_m|\{{\mathcal{A}}_{i_k{j}_k}{\}}_{k=1}^m)\right)$$

$$=\mathbb{E}\left(\prod _{k=1}^m\mathrm{I}\mathrm{P}\left({\mathcal{C}}_{i_k{j}_k}^2=y_k|\{{\mathcal{A}}_{i_k{j}_k}{\}}_{k=1}^m\right)\right)=\mathbb{E}\left(\prod _{k=1}^m\mathrm{I}\mathrm{P}\left({\mathcal{C}}_{i_k{j}_k}^2=y_k|{\mathcal{A}}_{i_k{j}_k}\right)\right)$$

$$=\mathbb{E}\left(\prod _{k=1}^m\mathrm{I}\mathrm{P}\left({\mathcal{C}}_{i_k{j}_k}^2=y_k\right)\right)=\prod _{k=1}^m\mathrm{I}\mathrm{P}\left({\mathcal{C}}_{i_k{j}_k}^2=y_k\right),$$

$$\mathbb{E}[R_{\mathbb{T},\mu }]\ge C_{\underline{d},\overline{d},\overline{\mu }}n|\mu -{\mu }^{\ast }{|}^2$$

$$\text{TV}(Q_{k,{\mu }^{\ast }}\Vert Q_{k,\mu })$$

$$=|e^{-2{\mu }^{\ast }{d}_{i_k{j}_k}}-e^{-2\mu d_{i_k{j}_k}}|\ge C_{\underline{d},\overline{d},\overline{\mu }}|{\mu }^{\ast }-\mu |.$$

$$\mathbb{E}[R_{\mathbb{T},\mu }]$$

$$\ge 2\sum _{k=1}^m[\text{TV}(Q_{k,{\mu }^{\ast }}\Vert Q_{k,\mu }){]}^2\ge m{C}_{\underline{d},\overline{d},\overline{\mu }}|{\mu }^{\ast }-\mu {|}^2.$$

$$\mathrm{Var}[u(X)]\le \omega \mathrm{Var}[v(X)]+\frac{2\omega }{\omega -1}{c}^2.$$

$$f(x)=\max \left\{\log 2,\log \frac{1}{1-e^{-2x}}\right\},S_{\mathbb{T},\mu }=\sum _{e\in E(\mathbb{T})}f(\mu t_e{)}^2,$$

$$\mathrm{Var}{R}_{\mathbb{T},\mu }\le C_{\alpha ,\underline{\mu },\overline{\mu }}{n}^{\alpha }+\frac{8{\omega }_{\alpha ,\beta }}{{\omega }_{\alpha ,\beta }-1}{n}^{\beta }{S}_{\mathbb{T},\mu }$$

$$R_{\mathbb{T},\mu }\le 2C|{\mu }^{\ast }-\mu |\le 2C|\overline{\mu }-\underline{\mu }|$$

$$\mathrm{Var}{R}_{\mathbb{T},\mu }=\mathrm{Var}{R}_{{\mathbb{T}}^{\prime },\mu }=\mathrm{Var}(R_{{\mathbb{T}}^{\prime },\mu }|I).$$

$$\mathrm{Var}{R}_{\mathbb{T},\mu }=\mathrm{Var}(R_{{\mathbb{T}}^{\prime },\mu }|I)\le \omega \mathrm{Var}(R_{{\mathbb{T}}_1,\mu }+R_{{\mathbb{T}}_2,\mu }|I)+\frac{8\omega }{\omega -1}f(\mu t_e{)}^2.$$

$$\mathrm{Var}{R}_{\mathbb{T},\mu }$$

$$\le \omega C_{\alpha }\frac{2^{\alpha +1}}{3^{\alpha }}{n}^{\alpha }+\frac{8\omega }{\omega -1}\left[\omega \frac{2^{\beta }}{3^{\beta }}(S_{{\mathbb{T}}_1}+S_{{\mathbb{T}}_2})+f(\mu t_e{)}^2\right]n^{\beta }.$$

$$S_{\mathbb{T}}-[S_{{\mathbb{T}}_1}+S_{{\mathbb{T}}_2}+f(\mu t_e{)}^2]=f(\mu t_{e_1}{)}^2+f(\mu t_{e_2}{)}^2-f[\mu (t_{e_1}+t_{e_2}){]}^2\ge 0$$

$$1<\omega \le \min \left\{{\left(\frac{3}{2}\right)}^{\beta },\frac{1}{2}{\left(\frac{3}{2}\right)}^{\alpha }\right\},$$

$$\mathrm{Var}{R}_{\mathbb{T},\mu }\le C_{\alpha }{n}^{\alpha }+\frac{8\omega }{\omega -1}{n}^{\beta }{S}_{\mathbb{T}},$$

$$\mathrm{Var}{R}_{\mathbb{T},\mu }\le C_{\alpha ,\gamma }{n}^{\alpha }$$

Full text

1 Introduction

The maximum likelihood estimator (MLE) is frequently used to estimate unknown parameters in trait evolution models on phylogenetic trees. To safely use this machinery, it is important to know that the MLE is consistent: that is, the estimate converges to the true value as we observe trait values for more species. However, traditional statistical theory only guarantees the consistency of the MLE when observations are independent and identically distributed. In contrast, trait values of biological species are not independent because species are related to each other according to a phylogenetic tree. Therefore, traditional consistency results are not directly applicable to trait evolution studies. For this paper we consider the scenario where only a single trait is observed and the number of species increases to infinity. This is different from the setting where many traits/characters are observed for the same set of species; this alternate setting typically assumes independence between traits.

Binary traits studied by comparative biologists come in various types, including morphological traits (whether a fruit fly has curly wings), behavioral traits (whether an antelope hides from predators), and geographical traits (whether a species is aquatic). Here we study the consistency property of the MLE for estimating the transition rate of a binary trait evolved along a phylogenetic tree according to a $2$-state symmetric Markov process. Under this $2$-state symmetric model, there is an unique rate of switching back and forth between the two states of a phenotype, which is the parameter of interest. We show that under two mild conditions, the MLE of this transition rate is consistent, that is, the estimate converges to the correct value as we observe the trait values of more species. These conditions ensure that edge lengths of the tree are not too small and the pairwise distances between leaves are neither too small nor too large. We verify that these conditions holds for many practical cases including trees with bounded edges, trees generated from the Yule process (Yule, 1925), and trees generated from the coalescent point process (Lambert and Stadler, 2013).

The consistency of the MLE does not always hold under trait evolution models. Indeed, for estimating the ancestral state of the $2$-state symmetric model, Li et al. (2008) point out that the MLE can be inconsistent. As a consequence, the estimate of the ancestral state may not be close to the true value no matter how many species have been sampled. Several efforts have also been made to investigate the consistency of the MLE under evolution models of a continuous trait. Ané (2008) points out that unlike traditional linear regression, the MLE of the coefficients of phylogenetic linear regression under the Brownian motion model can be inconsistent. Additionally, Sagitov and Bartoszek (2012) show that the sample mean is an inconsistent estimator for the ancestral state under this model when the tree is generated from the Yule process. For the Ornstein-Uhlenbeck model, Ho and Ané (2013) show that if the height of the phylogenetic tree is bounded as more observations are collected, then the MLE of the selective optimum is not consistent. Moreover, they discovered that in this scenario, no consistent estimator for the selective optimum exists. Recently, Ané et al. (2017) provide a necessary and sufficient condition for consistency of the MLE under the Ornstein-Uhlenbeck model. Although the problem of reconstructing the ancestral state under the $2$-state symmetric model has been studied extensively (Tuffley and Steel, 1997; Mooers and Schluter, 1999; Li et al., 2008; Mossel and Steel, 2014), it remains unknown whether the transition rate can be estimated consistently.

In this paper, we show that the MLE is a consistent estimator for the transition rate of the $2$-state symmetric model under simple conditions. We start with introducing the $2$-state symmetric model for binary traits and derive several statistical properties of this model in Section 2. In Section 3, we state two necessary conditions for the consistency of the MLE of the transition rate and provide a detailed proof for this result. Section 4 verifies these conditions for several practical scenarios and illustrates our result through a simulation.

2 Properties of the $2$-state symmetric model

Let $\mu$ be the transition rate of a $2$-state symmetric Markov process. Then, the transition probability matrix has an analytical form:

[TABLE]

Hence, the probability that the process switches state after $t$ unit of time is $\frac{1}{2}-\frac{1}{2}e^{-2\mu t}$.

In this paper, the term phylogenetic tree (or phylogeny) refers to a bifurcating rooted tree with leaves labeled by a set of taxon (species) names. We can reroot a phylogenetic tree by moving the root to another location along the tree (see Figure 1). The evolution of a binary trait along a tree is modeled using the $2$-state symmetric Markov process as follows. At each node in the phylogeny, the children inherit the trait value of their parent and the trait of each child evolves independently of one another.

Let $\mathbb{T}$ be a phylogenetic tree with $n$ leaves, and $\mathbf{Y}$ be the trait values at the leaves of $\mathbb{T}$. We assume that the ancestral state at the root $\rho$ of $\mathbb{T}$ follows a stationary distribution, which is a Bernoulli distribution with success probability $1/2$. Let $E$ be the set of all edges of $\mathbb{T}$. The joint probability distribution of $\mathbf{Y}$ is

[TABLE]

where $y$ ranges over all extensions of $\mathbf{Y}$ to the internal nodes of the tree, $y_{u}$ denotes the assigned state of node $u$ by $y$, $d_{uv}$ is the edge length of $(u,v)$, and $[\mathbf{P}_{\mu}(t)]_{kl}$ is the element at $k$-th row and $l$-th column of matrix $\mathbf{P}_{\mu}(t)$. We define the log-likelihood function as $\ell_{\mathbb{T},\mu}(\mathbf{Y})=\log P_{\mathbb{T},\mu}(\mathbf{Y})$. Throughout this paper, we denote the true transition rate with $\mu^{*}$ and assume that $\mu^{*}\in[\underline{\mu},\overline{\mu}]$ where $\underline{\mu},\overline{\mu}$ are two known positive numbers. Define

[TABLE]

Then the Kullback-Leibler divergence from $P_{\mathbb{T},\mu}$ to $P_{\mathbb{T},\mu^{*}}$ is

[TABLE]

where $\mathbb{E}$ is the expectation with respect to $P_{\mathbb{T},\mu^{*}}$.

Lemma 1**.**

Let $\mathbb{T}^{\prime}$ be the phylogenetic tree obtained by rerooting a tree $\mathbb{T}$, then

[TABLE]

for any trait values $\mathbf{Y}$ at the leaves.

This lemma, sometimes called the “pulley principle,” is a direct consequence of the fact that the ancestral state follows a stationary distribution (Felsenstein, 1981). We also have the following Lemmas:

Lemma 2**.**

Let $\mathbb{T}$ be a rooted tree with root $\rho$, and $\mathbf{Y}$ be the trait values at the leaves of $\mathbb{T}$ generated under the $2$-state symmetric model. Let $h$ be a function such that $h(\mathbf{Y})=h(\mathbf{1}-\mathbf{Y})$, where $\mathbf{1}$ denotes the all-ones vector. Then, $h(\mathbf{Y})$ and the trait value at $\rho$ are independent.

The following two lemmas concern the regularity of the log-likelihood function. In particular, Lemma 3 shows that the log-likelihood function of a binary tree can be bounded by the sum of the log-likelihood functions of its two subtrees.

Lemma 3**.**

Let $\rho_{1}$ and $\rho_{2}$ be two direct descendants of $\rho$, and $\mathbb{T}_{1}$ and $\mathbb{T}_{2}$ the two subtrees descending from them. Let $\mathbf{Y}_{1}$, $\mathbf{Y}_{2}$ be the observations at the leaves of $\mathbb{T}_{1}$ and $\mathbb{T}_{2}$. We have

[TABLE]

where $d$ is the tree distance between $\rho_{1}$ and $\rho_{2}$.

Lemma 4** (Uniform Lipschitz).**

There exists $C>0$ such that

[TABLE]

where $n$ denotes the number of leaves of $\mathbb{T}$.

The proofs of Lemmas 2, 3, and 4 are provided in the Appendix A. Henceforth we will also use $\rho$ to denote the trait value at the root by an abuse of notation.

3 Convergence of the MLE of the transition rate

The MLE of $\mu$ is defined as follows:

[TABLE]

In this section, we will state our main result regarding the consistency of $\hat{\mu}$. We will need to make two assumptions to ensure that the shape of $\mathbb{T}$ is not too irregular. Let us define

[TABLE]

where $t_{e}$ denotes the edge length of edge $e$.

Assumption 1**.**

$S_{\mathbb{T},\mu}=\mathcal{O}(n^{\gamma}$) for some $1\leq\gamma<2$. That is, there exists a universal constant $c<\infty$ such that $S_{\mathbb{T},\mu}<cn^{\gamma}$.

Assumption 2**.**

There exist $\Omega(n)$ pairs of leaves such that the paths connecting each pair are pairwise disconnected and their lengths are bounded in some fixed range $[\underline{d},\overline{d}]$. Here, $\Omega(n)$ denotes a quantity that is greater than $cn$ for a positive constant $c$ and all $n$.

Assumption 1 makes sure that the edge lengths of $\mathbb{T}$ are not too small. It is worth noticing that $1/(1-e^{-2x})\leq 1/x$ when $x$ is small enough. Therefore, this assumption holds when the smallest edge length is $\Omega(e^{-n(\gamma-1)/2})$. On the other hand, Assumption 2 guarantees that the pairwise distances between leaves of $\mathbb{T}$ do not vary too extremely. In Section 4, we will verify these assumptions for several common tree models. Although we employ rerooting in proofs below, Assumption 1 is for the original root of the tree.

Theorem 1**.**

Under Assumptions 1 and 2, for any $\delta>0$, $\forall\alpha\in\left(\max\left\{\frac{\log(2)}{\log(3/2)},\gamma\right\},2\right)$, there exists a constant $C_{\delta,\alpha,\underline{d},\overline{d},\underline{\mu},\overline{\mu}}>0$ such that

[TABLE]

with probability $1-\delta$.

From now on, we will use $\ell_{\mathbb{T},\mu}$ and $R_{\mathbb{T},\mu}$ as short for $\ell_{\mathbb{T},\mu}(\mathbf{Y})$ and $R_{\mathbb{T},\mu}(\mathbf{Y})$. Recall that $\ell_{\mathbb{T},\mu}$ is the log-likelihood function and $R_{\mathbb{T},\mu}=\ell_{\mathbb{T},\mu^{*}}-\ell_{\mathbb{T},\mu}$ where $\mu^{*}$ is the true value of the transition rate. The main ideas of the proof of Theorem 1 can be outlined as follows. For a fixed value of $\mu$, we can view the function $R_{\mathbb{T},\mu}$ as the evidence to distinguish between $\mu$ and $\mu^{*}$. We will prove that as the number of leaves approaches infinity, the information to distinguish $\mu$ and $\mu^{*}$, characterized by the KL divergence between $P_{\mathbb{T},\mu}$ to $P_{\mathbb{T},\mu^{*}}$ (which is equal to $\mathbb{E}(R_{\mathbb{T},\mu})$) increases linearly (Lemma 6), while the associated uncertainty, characterized by the variance of $R_{T,\mu}$, only increases sub-quadratically (Lemma 9). It is worth noting that while the results of Lemma 6 and Lemma 9 are straightforward for independent and identically distributed data, the analyses for phylogenetic traits are more complicated due to the correlations between trait values at the leaves. To overcome this issue, we use the independent phylogenetic contrasts, introduced in the next section, to obtain a lower bound on the information. On the other hand, Lemma 3 shows that the log-likelihood function of a binary tree can be bounded by the sum of the log-likelihood functions of its two subtrees, which allows us to exploit the sparse structure of the tree to derive an upper bound on uncertainty through an induction argument. Finally, we obtain a uniform bound on the difference between the log-likelihood functions and its expected values (Lemma 10), which enables us to derive an analysis of convergence of the MLE.

3.1 Lower bound on information

Phylogenetic contrasts

Letting $i$ and $j$ be two different species, we define $\mathcal{C}_{ij}=Y_{i}-Y_{j}$ to be a contrast between the two species. This is a popular notion introduced by Felsenstein (1985) for computing the likelihood function under the Brownian motion model. Ané et al. (2017) use independent contrasts to construct consistent estimators for the covariance parameters of the Ornstein-Uhlenbeck model. Here, we will introduce a notion of a squared-contrast $\mathcal{C}_{ij}^{2}$, which is simply the square of a contrast $\mathcal{C}_{ij}$, and show that squared-contrasts have the same independence properties as independent contrasts. Let $\mathcal{A}_{ij}$ be the state of the most recent common ancestor of $i$ and $j$, $d_{ij}$ the tree distance from $i$ to $j$, and $p_{ij}$ the path on the tree connecting $i$ and $j$. By symmetry of the model, we have

[TABLE]

Therefore, $\mathcal{C}_{ij}^{2}$ and $\mathcal{A}_{ij}$ are independent.

Lemma 5** (Sequence of contrasts).**

Let $\{\mathcal{C}_{i_{k}j_{k}}\}_{k=1}^{m}$ be a sequence of contrasts such that any two paths in $\{p_{i_{k}j_{k}}\}_{k=1}^{m}$ have no common node. Then $\{\mathcal{C}_{i_{k}j_{k}}^{2}\}_{k=1}^{m}$ are independent.

Proof.

We have, using $\mathbb{E}\!\left(\cdot\right)$ here to denote expectation over ancestral states,

[TABLE]

where the third equality comes from the Markov property. This completes the proof. ∎

The independent squared-contrasts allow us to derive the following lower bound on the information.

Lemma 6**.**

Under Assumption 2, there exists $C_{\underline{d},\overline{d},\overline{\mu}}>0$ such that

[TABLE]

for all $\mu\in[\underline{\mu},\overline{\mu}]$.

Proof.

Under Assumption 2, there exists a set of $m=\Omega(n)$ pairwise disjoint pairs (that is, the paths connecting each pair are pairwise disconnected) $(Y_{k,1},Y_{k,2})_{k=1}^{m}$. Consider the corresponding set of squared-contrasts $\mathcal{C}^{2}_{i_{k}j_{k}}=(Y_{k,1}-Y_{k,2})^{2}$. By Lemma 5, $(\mathcal{C}^{2}_{i_{k}j_{k}})_{k=1}^{m}$ are independent. Let $Q_{k,\mu}$ be the distribution of $\mathcal{C}^{2}_{i_{k}j_{k}}$ corresponding to parameter $\mu$. The total variation from $Q_{k,\mu^{*}}$ to $Q_{k,\mu}$ is

[TABLE]

By the data processing inequality (Theorem 9 in Van Erven and Harremoës, 2014), the fact that $(\mathcal{C}^{2}_{i_{k}j_{k}})_{k=1}^{m}$ are independent, and Pinsker’s inequalities, we have

[TABLE]

Note that $m=\Omega(n)$ which completes the proof. ∎

3.2 Upper bound on the uncertainty

The result of Section 3.1 indicates that as the number of leaves increases, the information to distinguish $\mu$ and $\mu^{*}$ also increases linearly. In order to prove that the MLE can successfully reconstruct the true parameter $\mu$, we need to ensure that such information is not confounded by the uncertainty associated with the log-likelihood function due to randomness in the data. In other words, we want to guarantee that as $n\to\infty$, the variance of $R_{\mathbb{T},\mu}$ only grows sub-quadratically, that is, there exist $\alpha<2$ and $C>0$ such that $\mathrm{Var}R_{\mathbb{T},\mu}\leq Cn^{\alpha}$. To obtain the bound, we need the following Lemma, the proof of which appears in Appendix A.

Lemma 7**.**

If two functions $u$ and $v$ satisfy $\|u-v\|_{\infty}\leq c$, then for all $\omega>1$ and random variables $X$,

[TABLE]

Recalling that,

[TABLE]

we have the following estimate that does not depend on Assumption 1.

Lemma 8**.**

For all $\alpha\in\left(\frac{\log(2)}{\log(3/2)},2\right)$ and $\beta\in(0,1)$, there exist $C_{\alpha,\underline{\mu},\overline{\mu}},\omega_{\alpha,\beta}>1$ such that

[TABLE]

for all trees $\mathbb{T}$ with $n\geq 2$ taxa and $\mu\in[\underline{\mu},\overline{\mu}]$.

Proof.

We will prove the result by induction in the number of taxa $n$. For $n=2$, we apply Lemma 4 to obtain

[TABLE]

where $C$ is the constant in Lemma 4. Thus $\mathrm{Var}R_{\mathbb{T},\mu}\leq 4C^{2}|\overline{\mu}-\underline{\mu}|^{2}$. Assume that the statement is valid for all $k<n$. We will prove that it is also valid for $k=n$. Now, let $\mathbb{T}$ be a bifurcating tree with $n$ taxa. Lipton and Tarjan (1979) show that there exists an edge $e=(I_{1},I_{2})$ of $\mathbb{T}$ such that if we reroot $\mathbb{T}$ at the middle of this edge, the two subtrees $\mathbb{T}_{1}$ and $\mathbb{T}_{2}$ stemming from $I_{1}$ and $I_{2}$ have no more than $2n/3$ leaves. Let $I$ be the middle point of edge $e$ and $\mathbb{T}^{\prime}$ be the tree obtained by rerooting $\mathbb{T}$ to $I$. Denote by $n_{1}$ and $n_{2}$ the number of leaves of $\mathbb{T}_{1}$ and $\mathbb{T}_{2}$. According to Lemma 1 and Lemma 2,

[TABLE]

By Lemma 3, we have $|R_{\mathbb{T}^{\prime},\mu}-R_{\mathbb{T}_{1},\mu}-R_{\mathbb{T}_{2},\mu}|\leq 2f(\mu t_{e})$. We apply Lemma 7 to obtain

[TABLE]

Let $\mathbf{Y}_{1}$ and $\mathbf{Y}_{2}$ be the observations at the leaves of $\mathbb{T}_{1}$ and $\mathbb{T}_{2}$ respectively. Since $\mathbf{Y}_{1}$ and $\mathbf{Y}_{2}$ are independent conditional on $I$, we have $\mathrm{Var}(R_{\mathbb{T}_{1},\mu}+R_{\mathbb{T}_{2},\mu}~{}|~{}I)=\mathrm{Var}(R_{\mathbb{T}_{1},\mu}~{}|~{}I)+\mathrm{Var}(R_{\mathbb{T}_{2},\mu}~{}|~{}I)$. By Lemma 2, we deduce that

[TABLE]

Recalling that the trees $I_{1}$ and $I_{2}$ have no more than $2n/3$ leaves, using the induction hypothesis for $\mathbb{T}_{1}$ and $\mathbb{T}_{2}$, we have

[TABLE]

Let $e_{1}$ and $e_{2}$ be the edges stemming from the root of $\mathbb{T}$. By definition, we have

[TABLE]

where the last inequality comes from the fact that $f$ is a non-increasing function. Thus, if we choose $\omega$ such that

[TABLE]

then

[TABLE]

which completes the proof. ∎

A combination of Lemma 8 and Assumption 1 gives rise to the desired bound:

Lemma 9** (Sub-quadratic upper bound on variance).**

Under Assumption 1, $\forall\alpha\in\left(\max\left\{\frac{\log(2)}{\log(3/2)},\gamma\right\},2\right)$, there exists $C_{\alpha,\gamma}>0$ such that

[TABLE]

for all $n\geq 2$ and $\mu\in[\underline{\mu},\overline{\mu}]$.

Proof.

By Lemma 8 and Assumption 1, we have

[TABLE]

Note that $1\leq\gamma<2$. So, $\forall\alpha\in\left(\max\left\{\frac{\log(2)}{\log(3/2)},\gamma\right\},2\right)$, we can choose $\beta=\alpha-\gamma$. ∎

3.3 Concentration bound

We are now ready to prove the concentration inequality for the 2-state symmetric model.

Lemma 10** (Concentration bound).**

Under Assumption 1, for any $\delta>0$ and $\forall\alpha\in\left(\max\left\{\frac{\log(2)}{\log(3/2)},\gamma\right\},2\right)$, there exists $C_{\delta,\alpha,\underline{\mu},\overline{\mu}}$ such that

[TABLE]

with probability at least $1-\delta$.

Proof.

Applying Chebyshev’s inequality, we obtain

[TABLE]

for any $\mu\in[\underline{\mu},\overline{\mu}]$. On the other hand, by Lemma 4, we have

[TABLE]

Therefore, if

[TABLE]

then

[TABLE]

for all

[TABLE]

Define

[TABLE]

We have

[TABLE]

The last inequality comes from Lemma 9. We complete the proof by picking

[TABLE]

∎

3.4 Proof of Theorem 1

From Lemma 6, we have

[TABLE]

where $E_{\mu^{*}}$ is the expectation with respect to $\bf P_{\mathbb{T},\mu^{*}}$.

Note that

[TABLE]

By Lemma 10, with probability $1-\delta$, we obtain

[TABLE]

which completes the proof.

4 Applications

In this section, we discuss applications of Theorem 1 for several practical scenarios. We first consider trees with edges of bounded length.

Theorem 2** (Trees with bounded edges).**

If edge lengths of $\mathbb{T}$ are bounded from below and above, then the MLE is consistent.

Proof.

Since all edge lengths are bounded from below, Assumption 1 holds with $\gamma=1$. Next, we will show that Assumption 2 is satisfied. We select $\lfloor n/2\rfloor$ pairs of leaves using the following procedure:

Pick a cherry as a pair. 2. 2.

Remove the cherry from the tree as well as the edge immediately above it. 3. 3.

Repeat step 1 and 2 until the tree has [math] or $1$ leaves.

Denote $\mathcal{C}=(i_{k},j_{k})_{k=1}^{m}$ be the set of pairs returned by the procedure where $m=\lfloor n/2\rfloor$. It is obvious that the paths connecting each pair are pairwise disjoint. Let $\underline{e}$ and $\overline{e}$ be the lower and upper bound of edge lengths, so we have a lower bound for distances between each pair $d_{i_{k}j_{k}}\geq 2\underline{e}$. Moreover,

[TABLE]

We will prove by contradiction that the set $\mathcal{B}=\{k:d_{i_{k}j_{k}}\leq 8\overline{e}\}$ has at least $n/8$ elements. Assume that $|\mathcal{B}|<n/8$. Note that

[TABLE]

We deduce that

[TABLE]

Therefore

[TABLE]

which contradicts (3). Therefore, $|\mathcal{B}|\geq n/8$.

∎

Next, we apply our result to trees generated from a pure-birth process (Yule, 1925). This is a classical tree-generating process which assumes that lineages give birth independently to one another with the same birth rate.

Theorem 3** (Yule process).**

If $\mathbb{T}$ is generated from a Yule process, then the MLE is consistent.

Proof.

Again, we only need to check Assumptions 1 and 2. Let $(t_{k})_{k=1}^{n}$ be the amount of time during which $\mathbb{T}$ has $k$ lineages. Then, $(t_{k})_{k=1}^{n}$ are independent exponential random variables with rate $k\lambda$ where $\lambda$ is the birth rate of the Yule process. Let $e_{\min}$ be the minimum edge length of $\mathbb{T}$. It is sufficient for us to show that $e_{\min}\geq e^{-n(\gamma-1)/2}$. Indeed, we have

[TABLE]

which converges to $1$ as $n\to\infty$. Hence, Assumption 1 is satisfied asymptotically with $\gamma>1$.

Define the age of an internal node to be its distance to the leaves. By Corollary 4 of Ané et al. (2017), there exist $0<c_{1}<c_{2}<+\infty$ such that the probability of there being $\Omega(n)$ internal nodes with age between $c_{1}$ and $c_{2}$ is $1-\mathcal{O}(n^{-1})$. Let $\mathcal{I}$ be the set of internal nodes which have age between $c_{1}$ and $c_{2}$, we select $\Omega(n)$ pairs of leaves as follows:

Start with an internal node $i\in\mathcal{I}$ that has the smallest age. Pick two of its descendants such that $i$ is their most recent common ancestor to form a pair. 2. 2.

Remove $i$ and its parent node from $\mathcal{I}$. 3. 3.

Repeat step 1 and 2 until $\mathcal{I}$ is empty.

This procedure is similar to the procedure for selecting contrasts in Ané et al. (2017, Lemma 3). Hence, it selects at least $|\mathcal{I}|/2$ pairs, which is of order $\Omega(n)$. Since these internal nodes have age between $c_{1}$ and $c_{2}$, Assumption 2 is satisfied.

∎

Finally, we consider the coalescent point process (CPP), which generates a random ultrametric tree with a given height $T$ (see Lambert and Stadler, 2013, for more details). Conditioning on the number of species $n$, the internal node ages $(t_{k})_{k=1}^{n-1}$ of the tree generated from the CPP are independent and identically distributed according to a probability distribution in $[0,T]$. The probability density function $\phi$ of this common distribution is called the coalescent density. We say that a CPP is regular if its common distribution is not a point mass at [math] and the coalescent density is bounded from above.

Theorem 4** (Coalescent point process).**

If $\mathbb{T}$ is generated from a regular coalescent point process, then the MLE is consistent.

Proof.

Denote the minimum edge length of $\mathbb{T}$ by $e_{\min}$. We have

[TABLE]

Therefore,

[TABLE]

By the results of Section 4 in Jammalamadaka and Janson (1986) (for which the details will be provided in the Appendix), we have

[TABLE]

for any $\epsilon>0$. We deduce that

[TABLE]

Since $\epsilon$ can be arbitrary small, we conclude

[TABLE]

On the other hand

[TABLE]

Note that $\phi$ is bounded from above by $M$ in $[0,T]$, thus

[TABLE]

where $M$ is the upper bound of $\phi$. Hence,

[TABLE]

Similarly,

[TABLE]

Therefore,

[TABLE]

Hence, Assumption 1 is satisfied.

Since the common distribution is not a point mass at [math], there exist a constant $c>0$ such that $\int_{c}^{T}{\phi}>0$. Let $\mathcal{I}$ be the set of internal nodes of $\mathbb{T}$ which have age between $c$ and $T$. Then, Assumption 2 holds if we can prove that $|\mathcal{I}|$ is of order $\Omega(n)$ because we can use the same argument as in the proof of Theorem 3. By strong law of large numbers, we have

[TABLE]

which completes the proof.

∎

5 Practical implication

The consistency property proved in Theorem 1 suggests that adding more taxa helps to significantly improve the accuracy of the MLE of the transition rate. It is worth noticing that this property does not hold in many scenarios (see Li et al., 2008; Ané, 2008; Ho and Ané, 2013, 2014). To illustrate our theoretical results, we perform the following simulation using the R package geiger (Harmon et al., 2007; Pennell et al., 2014). We simulate $100$ trees (the number of taxa varies from $50$ to $2000$) according to the Yule process with birth rate $\lambda=1$. For each tree, we simulate $100$ traits under the $2$-state symmetric model with the transition rate $\mu=0.5$. For each of these 100 traits, the MLE for the transition rate is computed separately. The result is summarized in Figure 2. We can see that the MLEs concentrate more and more around the true transition rate as the number of species increases.

6 Discussion

In this paper, we investigate the convergence of the MLE for the transition rate of a $2$-state symmetric model under regularity conditions on tree shape. These conditions ensure that edge lengths of the tree are not too small and the pairwise distances between leaves are not too extreme. For example, these conditions are satisfied when the edge lengths of the tree are bounded from above and below. We have also verified that trees generated from pure-birth process and coalescent point process satisfy these conditions. On the other hand, whether these sufficient conditions are also necessary conditions remains open.

Our results suggest that adding data in terms of the number of tips tends to improve the accuracy of the MLE of the transition rate in general. To demonstrate this result, we use simulations to confirm that the MLE of the transition rate is consistent when the tree is generated from the Yule process. It is worth noticing that for evolutionary data, the MLE for other quantities of interest in evolutionary studies can be inconsistent (Li et al., 2008; Ané, 2008; Ho and Ané, 2013, 2014). In these situations, adding more data is a waste of resources because it does not significantly improve the precision of the MLE. Therefore, it is important to study the consistency of the MLE in the context of trait evolution models.

Beside its practical biological application, the paper also investigates many theoretical properties of tree-generating processes. Specifically, in Theorem 3 and Theorem 4, we consider the problem of bounding the minimum edge length $e_{min}$ of a tree generated by Yule and coalescent point processes. The theorems show that generically, $e_{min}$ is bounded from below by a function of order $n^{-2}$ for the Yule process and $n^{-3}$ for the coalescent point process. Since the minimum edge length plays an important role in many phylogenetic problems, these results may be of independent interest.

Finally, we remark that if we have multiple independent traits which evolve according to the same 2-state Markov process, the distance from the MLE to the true value of the transition rate will decrease linearly with respect to the number of traits. This result comes from the fact that traditional properties of MLE can be applied because these traits are independent. Therefore, incorporating additional traits to the analysis is another way to improve the precision of the MLE.

Appendix A Proofs

A.1 Proof of lemma 2

Note that by symmetry, we have ${\mathrm{I}\!\mathrm{P}}(\mathbf{Y}=\mathbf{y}~{}|~{}\rho=0)={\mathrm{I}\!\mathrm{P}}(1-\mathbf{Y}=\mathbf{y}~{}|~{}\rho=1)$. We deduce that

[TABLE]

which completes the proof.

A.2 Proof of lemma 3

Denote $P^{(u)}_{v}={\mathrm{I}\!\mathrm{P}}(\mathbf{Y_{u}}~{}|~{}\mathbb{T}_{u},\mu,\rho_{u}=v)$ for $u\in\{0,1\}$, $v\in\{0,1\}$. We have

[TABLE]

Moreover

[TABLE]

Therefore

[TABLE]

A.3 Proof of lemma 4

Without loss of generality, we assume that $\mu_{1}<\mu_{2}$. By the mean value theorem, there exists $\tilde{\mu}_{uv}\in(\mu_{1},\mu_{2})$ for any $u,v\in\{0,1\}$ such that

[TABLE]

We observe that there exists a $C_{\underline{\mu},\overline{\mu}}>0$ such that

[TABLE]

Therefore,

[TABLE]

This implies that

[TABLE]

Note that

[TABLE]

By applying (5) for all $2n-3$ edges on the tree, we deduce that

[TABLE]

Hence,

[TABLE]

which validates the lemma.

A.4 Proof of lemma 7

For all $x,y$, we have $|u(x)-u(y)|\leq|v(x)-v(y)|+2c$. Let $Y$ be an independent and identically distributed copy of $X$, we have

[TABLE]

Note that for all $z,c\in\mathbb{R}$ and $\omega>1$,

[TABLE]

Therefore,

[TABLE]

A.5 Proof of Equation (4)

In order to establish Equation (4), we use the following Lemma.

Lemma 11** **(Remark 3.4 in Jammalamadaka and

Janson (1986)).

Let $X_{1},X_{2},\ldots,X_{n}$ be an i.i.d. sequence of random variables and $f_{n}(x,y)$ be an indicator function on $\mathbb{R}^{2}$ such that

[TABLE]

for some constant $\lambda>0$. Define $U_{n}=\sum_{1\leq i<j\leq n}{f_{n}(X_{i},X_{j})}$.

Then $U_{n}\to_{d}\text{Poisson}(\lambda)$.

We apply this Lemma with $f_{n}(x,y)=I\{|x-y|<r_{n}\}$ where $r_{n}=\epsilon/n^{2}$ for the sequence $t_{1},t_{2},\ldots,t_{n}$ of the coalescent point process. Note that by Equation (4.3) in Jammalamadaka and Janson (1986),

[TABLE]

for some constant $c>0$. On the other hand, we have

[TABLE]

Therefore, $n^{3}E[f_{n}(t_{1},t_{2})f_{n}(t_{1},t_{3})]\to 0$. Hence,

[TABLE]