---
abstract: |
    Diffusion processes on trees are commonly used in evolutionary biology to model the joint distribution of continuous traits, such as body mass, across species. Estimating the parameters of such processes from tip values presents challenges because of the intrinsic correlation between the observations produced by the shared evolutionary history, thus violating the standard independence assumption of large-sample theory. For instance Ho and Ané [@HoAne13] recently proved that the mean (also known in this context as selection optimum) of an Ornstein-Uhlenbeck process on a tree cannot be estimated consistently from an increasing number of tip observations if the tree height is bounded. Here, using a fruitful connection to the so-called reconstruction problem in probability theory, we study the convergence rate of parameter estimation in the unbounded height case. For the mean of the process, we provide a necessary and sufficient condition for the consistency of the maximum likelihood estimator (MLE) and establish a phase transition on its convergence rate in terms of the growth of the tree. In particular we show that a loss of $\sqrt{n}$-consistency (i.e., the variance of the MLE becomes $\Omega(n^{-1})$, where $n$ is the number of tips) occurs when the tree growth is larger than a threshold related to the phase transition of the reconstruction problem. For the covariance parameters, we give a novel, efficient estimation method which achieves $\sqrt{n}$-consistency under natural assumptions on the tree. Our theoretical results provide practical suggestions for the design of comparative data collection.

    #### Keywords

    Ornstein-Uhlenbeck, phase transition, evolution, phylogenetic, consistency, maximum likelihood estimator.
author:
- 'Cécile Ané[^1]'
- 'Lam Si Tung Ho[^2]'
- 'Sebastien Roch[^3]'
title: 'Phase transition on the convergence rate of parameter estimation under an Ornstein-Uhlenbeck diffusion on a tree'
---

Introduction
============

Analysis of data collected from multiple species presents challenges because of the intrinsic correlation produced by the shared evolutionary history. This dependency structure can be modeled by assuming that the traits of interest evolved along a phylogeny according to a stochastic process. Two commonly used processes for continuous traits, such as body mass, are Brownian motion (BM) and the Ornstein-Uhlenbeck (OU) process. BM is used to model neutral evolution, with no favored direction (see e.g. [@felsenstein1985phylogenies]). On the other hand, the OU process can account for natural selection using two extra parameters: a “selection optimum” $\mu$ towards which the process is attracted and a “selection strength” $\alpha$ [@hansen1997stabilizing]. The OU process has a stationary distribution, which is Gaussian with mean $\mu$ and variance $\gamma = {\sigma^2/2 \alpha}$. The presence of natural selection can be detected by testing whether $\alpha>0$ (e.g. [@harmon-etal10]). Changes in $\mu$ across different groups of organisms are used to correlate changes in selection regime with changes in behavior or environmental conditions (see e.g. [@butler2004phylogenetic; @bartoszek2012multivariate]). For instance, the optimal body size $\mu$ might be different for terrestrial animals than for birds and bats. In practice, $\mu$, $\alpha$ and the infinitesimal variance $\sigma^2$ (or stationary variance $\gamma$) are estimated from data on extant species. In other words, only data at the tips of the tree are available. The process at internal nodes and edges is unobserved. Also, the tree is reconstructed independently from external and abundant data, typically from DNA sequences. In practice there can be some uncertainty about a few nodes in the tree, but we assume here that the tree is known without error.

The OU process on a tree has been used extensively in practice (see e.g. [@butler2004phylogenetic; @cooperPurvis10; @brawand_etal11; @rohlfs2014modeling]), but very few authors have studied convergence rates of available estimators. Recently Ho and Ané [@HoAne13] showed that if the tree height is bounded as the sample size goes to infinity, no estimator for $\mu$ can ever be consistent. This is because $\mu$ is not “microergodic": the distribution $P_\mu$ of the whole observable process $(Y_i)_{i\geq 1}$ at the tips of the tree is such that $P_{\mu_1}$ and $P_{\mu_2}$ are not orthogonal for any values $\mu_1\neq\mu_2$, if the tree height is bounded. This boundedness assumption does not hold for common models of evolutionary trees however, such as the pure-birth (Yule) process [@yule1925mathematical]. We consider here the case of an unbounded tree height. We study the consistency and convergence rates of several estimators, including some novel estimators, using tools from the literature on the reconstruction problem in probability theory. In particular we relate the convergence rates of these estimators to the growth rate of the phylogeny. This connection is natural given that the growth rate (and the related branching number) is known to play an important role in the analysis of a variety of stochastic processes including random walks, percolation and ancestral state reconstruction on trees [@peres1999climb]. In particular we leverage a useful characterization of the variance of linear estimators in terms of electrical networks.

#### Main results

We present the asymptotic properties of two common estimators for $\mu$: the sample mean and the maximum likelihood estimator (MLE). Conditional on the tree, the MLE ${\hat{\mu}_{\mathrm{ML}}}$ is known to be the best linear unbiased estimator for $\mu$ assuming that $\alpha$ is known. (The assumption of known $\alpha$ is proved not to be restrictive for our convergence rate results if $\alpha$ can be well estimated.) In fact, we give an example when ${\hat{\mu}_{\mathrm{ML}}}$ performs significantly better than the sample mean, which is not consistent in that particular case. In one of our main results, we identify a necessary and sufficient condition for the consistency of ${\hat{\mu}_{\mathrm{ML}}}$. We also derive a phase transition on its convergence rate, which drops from $\sqrt{n}$-consistency (i.e. the variance is $O(n^{-1})$) to a lower rate, $n$ being the number of samples (i.e. tip observations). This phase transition depends on the growth rate of the tree. Tree growth measures the rate at which new leaves arise as the tree height increases (see Section \[sec:def\] for a formal definition). Roughly, when the growth rate is below $2\alpha$, we show that $\sqrt{n}$-consistency holds. This is intuitive as a lower growth rate means lower correlations between the leaf states. On the other hand, when the growth rate is above $2\alpha$ implying a sample size $n \gg e^{2\alpha T}$, i.e. when the tree is sufficiently “bushy,” then the “effective sample size" is reduced to ${n^{\mathrm{eff}}}= e^{2\alpha T}$ and the $\sqrt n$-consistency of ${\hat{\mu}_{\mathrm{ML}}}$ is lost. We also provide novel, efficient estimators for the other two parameters, $\alpha$ and $\gamma$, which achieve $\sqrt{n}$-consistency and do not require the knowledge of $\mu$. Interestingly, the $\sqrt{n}$-consistency in this case is not affected by growth rate, unlike the case of the MLE for $\mu$. Our results lead to a practical method to assess whether additional species are informative or not, thus helping researchers to avoid wasting money and effort. Section \[sec:app\] presents simulations to illustrate these suggestions. Our main results are stated formally and further discussed in Section \[sec:def\], after necessary definitions. Their proofs are found in Section \[sec:proofs\].

#### Related work

Bartoszek and Sagitov [@bartoszek2012phylogenetic] obtained a corresponding phase transition for the convergence rate of the sample mean to estimate $\mu$, assuming a Yule process for the tree. Phase transitions for the convergence rate of some U-statistics have also been obtained for the OU model when the tree follows a supercritical branching process [@adamczak2011clt; @adamczak2011u]. A main difference between these studies and our work is that we assume that the tree is known. Even though tree-free estimators are the only practical options when the tree is unknown, this situation is now becoming rare due to the ever-growing availability of sequence data for building trees. For instance Crawford and Suchard [@crawford2013diversity] acknowledge that “as evolutionary biologists further refine our knowledge of the tree of life, the number of clades whose phylogeny is truly unknown may diminish, along with interest in tree-free estimation methods.”

As we mentioned, related phase transitions have been obtained for other processes on trees. For instance, the growth rate of the tree determines whether the state at the root can be reconstructed better than random for a binary symmetric channel on a binary tree (see e.g. [@EvKePeSc:00] and references therein). In a recent result, Mossel and Steel [@mossel2014majority] established a transition for ancestral state reconstruction by majority rule for the binary symmetric model on a Yule tree at the same critical point as above. Note that majority rule is a tree-free estimator like the sample mean in [@bartoszek2012phylogenetic], but adapted to discrete traits. In the context of the OU model, Mossel et al. [@mossel2013robust] obtained a phase transition for estimating the ancestral state at the root, with the same critical growth rate we derive in our results.

Definitions and statements of results {#sec:def}
=====================================

In this section, we state formally and further explain our main results. First, we define our model and describe the setting in which our results are proved.

Model
-----

Our main model is a stochastic process on a species tree ${\mathbb{T}}$. Let ${\mathbb{T}}= (\mathscr{E},\mathscr{V})$ be a finite tree with leaf set ${\mathscr{L}}= \{1,\ldots,n\}$ and root $\rho$. The leaves typically correspond to extant species. We think of the edges of ${\mathbb{T}}$ as being oriented away from the root. To each edge (or branch) $b \in \mathscr{E}$ of the tree is associated a positive length $|b| > 0$ corresponding to the time elapsed between the endpoints of $b$. For any two vertices $u,v \in \mathscr{V}$, we denote by $d_{uv}$ the distance between $u$ and $v$ in ${\mathbb{T}}$, that is, the sum of the branch lengths on the unique path between $u$ and $v$. We assume that the species tree is [*ultrametric*]{}, that is, that the distance from the root to every leaf is the same. It implies that, for any two tips $i,j \in {\mathscr{L}}$, $d_{ij}$ is twice the time to the most recent common ancestor of $i$ and $j$ from the leaves. We let $T$ be the height of ${\mathbb{T}}$, that is, the distance between the root and any leaf, and we define $t_{ij} = T - \frac{d_{ij}}{2}$. [*Throughout we assume that the species tree is known.*]{}

We consider an Ornstein-Uhlenbeck (OU) process on ${\mathbb{T}}$. That is, on each branch of ${\mathbb{T}}$, we have a diffusion $$dY_t = -\alpha(Y_t - \mu)dt + \sigma dB_t,$$ where $B_t$ is a standard Brownian motion (BM). In the literature on continuous traits, $Y_t$ is known as the response variable, $\mu$ is the selection optimum, $\alpha > 0$ is the selection strength, $\sigma > 0$ is the scale parameter of the Brownian motion. We assume that the root value follows the stationary Gaussian distribution ${\cal N}(\mu,\gamma)$, where $\gamma = \frac{\sigma^2}{2 \alpha}$. At each branching point, we run the process independently on each descendant edge starting from the value at the branching. Equivalently, the column vector of observations $\mathbf{Y} = (Y_\ell)_{\ell \in {\mathscr{L}}}$ at the tips of the tree are Gaussian with mean $\mu$ and variance matrix $\mathbf{\Sigma} = \gamma \mathbf{V_{\mathbb{T}}}$ where $$(V_{\mathbb{T}})_{ij}=e^{-\alpha d_{ij}}.$$ [*We assume throughout that $\alpha$, $\mu$ and $\sigma$ are the same on every branch of ${\mathbb{T}}$. We will specify below whether these parameters are known, depending on the context.*]{}

#### Parameter estimators

Our interest lies in estimating the parameters of the model, given ${\mathbb{T}}$, from a sample of $\mathbf{Y}$. In addition to proposing new estimators for $\alpha$ and $\sigma$, we study common estimators of $\mu$. In particular we consider the empirical average at the tips $
\overline{Y} = {\mathbf{1}}'\mathbf{Y} /n,
$ where ${\mathbf{1}}$ denotes the all-ones vector and $\mathbf{v}'$ denotes the transposes of a vector or matrix $\mathbf{v}$. Also, the MLE of $\mu$ [*given the tree and $\alpha$*]{} is $${\hat{\mu}_{\mathrm{ML}}}= ({\mathbf{1}}' \mathbf{V}_{\mathbb{T}}^{-1} {\mathbf{1}})^{-1}  {\mathbf{1}}' \mathbf{V}_{\mathbb{T}}^{-1} \mathbf{Y},$$ which is the well-known generalized least squares estimator for the linear regression problem $\mathbf{Y} =  \mu {\mathbf{1}}+ \boldsymbol{\varepsilon},$ where $\boldsymbol{\varepsilon}$ is multivariate normal with covariance matrix $\mathbf{\Sigma}$ (see e.g. [@Anderson:1984]). Note that the mean squared error is given by $${\mathrm{Var}}_{\mathbb{T}}[{\hat{\mu}_{\mathrm{ML}}}] 
= ({\mathbf{1}}' \mathbf{V}_{\mathbb{T}}^{-1} {\mathbf{1}})^{-2} {\mathbf{1}}' \mathbf{V}_{\mathbb{T}}^{-1}\mathbf{\Sigma} (\mathbf{V}_{\mathbb{T}}^{-1})' {\mathbf{1}}= \gamma ({\mathbf{1}}' \mathbf{V}_{\mathbb{T}}^{-1} {\mathbf{1}})^{-1}.
\label{eq:varmuhat}$$ We drop the ${\mathbb{T}}$ in ${\mathrm{Var}}_{\mathbb{T}}$ when the tree is clear from the context.

The estimators $\overline{Y}$ and ${\hat{\mu}_{\mathrm{ML}}}$ are both linear estimators. It is useful to think of the MLE in this context as an unbiased linear estimator minimizing the mean squared error (that is, a best linear unbiased estimator), which follows from the Gauss-Markov Theorem [@shao2003mathstat].

Asymptotic setting
------------------

Our results are asymptotic. Specifically, we consider sequences of trees ${\mathcal{T}}= ({\mathbb{T}}_k)_{k \geq 1}$ with [*fixed*]{} parameters $\alpha, \mu, \sigma$. For $k \geq 1$, let $n_k$ be the number of leaves in ${\mathbb{T}}_k$ and $T_k$ be the height of ${\mathbb{T}}_k$. As before, we denote the leaf set of ${\mathbb{T}}_k$ as ${\mathscr{L}}_k = [n_k]$.

\[assumption:unbounded\] Throughout we assume that $n_k \leq n_{k+1}$, $T_k \leq T_{k+1}$, and that $n_k \to +\infty$ and $T_k \to +\infty$ as $k \to +\infty$.

For such a sequence of trees and a corresponding sequence of estimators, say $X_k$, we recall various desirable asymptotic properties of $X_k$.

Let $(X_k)_k$ be a sequence of estimators for a parameter $x$. We say that $(X_k)_k$ is [*consistent for $x$*]{} if $X_k$ converges in probability to $x$, denoted as $|X_k - x| = o_p(1)$. For $\beta > 0$, we say that $(X_k)_k$ is [*$(n_k^{\beta})$-consistent for $x$*]{} if $(n_k^{\beta}(X_k-x))_k$ is bounded in probability, which we denote as $|X_k - x| = O_p(n_k^{-\beta})$.

We also recall the following notation. Let $(x_k)_k$ and $(y_k)_k$ be two sequences of real numbers. We let $y_k = O(x_k)$ if there exists $C_1 > 0$ such that $|y_k| \leq C_1 |x_k|$; $y_k = \Omega(x_k)$ if there exists $C_2 > 0$ such that $|y_k| \geq C_2 |x_k|$; and $y_k = \Theta(x_k)$ if $y_k = O(x_k)$ and $y_k = \Omega(x_k)$.

#### Growth

Our asymptotic results depend on how fast the tree grows. We first provide some intuition through a toy example.

\[ex:firstphase\] Let ${\mathbb{T}}_k$ be a star tree with $n_k$ leaf edges of length $T_k$ emanating from the root. By symmetry, ${\mathbf{1}}$ is an eigenvector of $\mathbf{\Sigma}$ with eigenvalue $\lambda_k = \gamma[1+(n_k-1)e^{-2\alpha T_k}]$. Hence, ${\mathbf{1}}$ is also an eigenvector of $\mathbf{\Sigma}^{-1}$ with eigenvalue $\lambda_k^{-1}$ and $ {\mathbf{1}}'\mathbf{\Sigma}^{-1}{\mathbf{1}}= n_k\lambda_{n_k}^{-1},$ so that ${\hat{\mu}_{\mathrm{ML}}}^{(k)} = \overline{Y}$ and $$\label{eq:varmlestar}
    {\mathrm{Var}}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}]  = \frac{\lambda_k}{n_k}      =  \gamma\left[e^{-2\alpha T_k} + \frac{1 - e^{-2\alpha T_k}}{n_k}\right].$$ If both $n_k$ and $T_k \to +\infty$, then ${\mathrm{Var}}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}] \to 0$ and the MLE (and $\overline{Y}$) is consistent for $\mu$. Furthermore, if $$\liminf_k \frac{2\alpha T_k}{\log n_k} > 1,$$ then $$n_k {\mathrm{Var}}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}] \leq \gamma[n_k e^{-2\alpha T_k} + 1] = O(1)$$ and the MLE is $\sqrt{n_k}$-consistent (by an application of Chebyshev’s inequality). On the other hand, if $$\liminf_k \frac{2\alpha T_k}{\log n_k} < 1,$$ then $$n_k {\mathrm{Var}}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}] \geq \gamma[n_k e^{-2\alpha T_k}],$$ which goes to $+\infty$ along a subsequence, and the MLE is [*not*]{} $\sqrt{n_k}$-consistent (using that ${\hat{\mu}_{\mathrm{ML}}}$ is unbiased and normally distributed).

To study more general trees, we use several standard notions of growth, which play an important role in random walks, percolation and ancestral state reconstruction on trees (see e.g. [@peres1999climb]).

\[def:growth\] The [*lower growth*]{} and [*upper growth*]{} of a tree sequence ${\mathcal{T}}$ are defined respectively as $${\underline{\Lambda}^\mathrm{g}}= \liminf_k \frac{\log n_k}{T_k},
\quad\mbox{and}\quad
{\overline{\Lambda}^\mathrm{g}}= \limsup_k \frac{\log n_k}{T_k}.$$ In case of equality we define the growth ${\Lambda^\mathrm{g}}= {\underline{\Lambda}^\mathrm{g}}= {\overline{\Lambda}^\mathrm{g}}$. (Note that our definition differs slightly from [@peres1999climb] in that we consider the “exponential rate” of growth.)

That is, for all $\epsilon > 0$, eventually $
e^{({\underline{\Lambda}^\mathrm{g}}- \epsilon)T_k} \leq
n_k \leq e^{({\overline{\Lambda}^\mathrm{g}}+ \epsilon)T_k},
$ and along appropriately chosen subsequences $
n_{k_j} \geq e^{({\overline{\Lambda}^\mathrm{g}}- \epsilon)T_{k_j}}
$ and $
n_{k'_{j}} \leq e^{({\underline{\Lambda}^\mathrm{g}}+ \epsilon)T_{k'_j}}.
$

We also need a stronger notion of growth. For a tree ${\mathbb{T}}$, thinking of the branches of ${\mathbb{T}}$ as a continuum of [*points*]{}, a cutset $\pi$ is a set of points of ${\mathbb{T}}$ such that all paths from the root to a leaf must cross $\pi$. Let $\Pi^k$ be the set of cutsets of ${\mathbb{T}}_k$.

The [*branching number*]{} of ${\mathcal{T}}$ is defined as $${\Lambda^\mathrm{b}}= \sup \left \{ \Lambda \geq 0\ :\ 
\inf_{k, \pi \in \Pi^k} 
\sum_{x \in \pi}{ e^{- \Lambda \delta_k(\rho,x)}  > 0 } \right \},$$ where $\delta_k(\rho,x)$ is the length of the path from the root to $x$ in ${\mathbb{T}}_k$.

Because the leaf set ${\mathscr{L}}_k$ forms a cutset, it holds that $${\Lambda^\mathrm{b}}\leq {\underline{\Lambda}^\mathrm{g}}\leq {\overline{\Lambda}^\mathrm{g}}.$$ Unlike the growth, the branching number takes into account aspects of the “shape” of the tree.

\[ex:first-transition-continued\] Consider again the setup of Example \[ex:firstphase\]. The infimum $$\inf_{\pi \in \Pi^k} 
\sum_{x \in \pi} e^{- \Lambda \delta_k(\rho,x)},$$ is achieved by taking $\pi = {\mathscr{L}}_k$ for every $k$. Hence ${\Lambda^\mathrm{b}}= {\underline{\Lambda}^\mathrm{g}}$. We showed in Example \[ex:firstphase\] that the MLE of $\mu$ given $\alpha$ is $\sqrt{n_k}$-consistent if ${\overline{\Lambda}^\mathrm{g}}< 2\alpha $, but not $\sqrt{n_k}$-consistent if ${\overline{\Lambda}^\mathrm{g}}> 2\alpha $.

Finally, we will need a notion of uniform growth.

Let ${\mathcal{T}}= ({\mathbb{T}}_k)_k$ be a tree sequence. For any point $x$ in ${\mathbb{T}}_k$, let $n_k(x)$ be the number of leaves below $x$ and let $T_k(x)$ be the distance from $x$ to the leaves. Then the [*uniform growth*]{} of ${\mathcal{T}}$ is defined as $${\Lambda^\mathrm{ug}}= \lim_{M\to+\infty} 
\sup_{k, x \in {\mathbb{T}}_k}\ \frac{\log n_k(x)}{T_k(x)\lor M}.$$ (The purpose of the $M$ in the denominator is to alleviate boundary effects.)

Statement of results {#sec:statementresults}
--------------------

We can now state our main results.

#### Results concerning the mean $\mu$

We first give a characterization of the consistency of the MLE of $\mu$. In words, the MLE sequence is consistent if, in the limit, we can find arbitrarily many descendants, arbitrarily far away from the leaves. This theorem is proved in Section \[sec:consistency-mle\], along with a related result involving the branching number.

\[thm:criterion-consistency\] Let $({\mathbb{T}}_k)_k$ be a sequence of trees satisfying Assumption \[assumption:unbounded\]. Let $({\hat{\mu}_{\mathrm{ML}}}^{(k)})_k$ be the corresponding sequence of MLEs of $\mu$ given $\alpha$. Denote by $\tilde{\pi}_t^{k}$ the cutset of ${\mathbb{T}}_k$ [*at time $t$ away from the leaves*]{} and let $T_k$ be the height of ${\mathbb{T}}_k$. Then $({\hat{\mu}_{\mathrm{ML}}}^{(k)})_k$ is consistent for $\mu$ if and only if for all $s \in (0,+\infty)$ $$\label{eq:consistency-condition}
\liminf_k \left|\tilde{\pi}^k_{s}\right| = +\infty.$$

We further obtain bounds on the variance of the MLE to characterize the rate of convergence of the MLE. In particular we give conditions for $\sqrt{n_k}$-consistency. We show that the latter undergoes a phase transition, generalizing Example \[ex:first-transition-continued\]. When the upper growth is above $2\alpha$, we show that the MLE of $\mu$ cannot be $\sqrt{n_k}$-consistent. If further the branching number is above $2\alpha$, we give tight bounds on the convergence rate of the MLE. Roughly we show that, in the latter case, the variance behaves like $n_k^{2\alpha/{\Lambda^\mathrm{g}}}$. Or perhaps a more accurate way to put it is that the “effective number of samples” ${n^{\mathrm{eff}}}_k$ is $e^{2\alpha T_k}$, in the sense that ${\mathrm{Var}}_{{\mathbb{T}}_k}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}] = \Theta(({n^{\mathrm{eff}}}_k)^{-1})$.

\[thm:rate-supercritical\] Let $({\mathbb{T}}_k)_k$ be a tree sequence. If ${\overline{\Lambda}^\mathrm{g}}> 2\alpha$, then for all $\epsilon > 0$ there is a subsequence $(k_j)_j$ along which $$\label{eq:non-sqrt-supercritical}
{\mathrm{Var}}_{{\mathbb{T}}_{k_j}}[{\hat{\mu}_{\mathrm{ML}}}^{(k_j)}] 
\geq \gamma n_{k_j}^{-2\alpha/({\overline{\Lambda}^\mathrm{g}}- \epsilon)}.$$ In particular $({\hat{\mu}_{\mathrm{ML}}}^{(k)})_k$ is not $\sqrt{n_k}$-consistent. If, further,

1.  ${\Lambda^\mathrm{b}}> 2\alpha$: then $${\mathrm{Var}}_{{\mathbb{T}}_k}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}] = \Theta\left(e^{-2\alpha T_k}\right).$$ Moreover in terms of $n_k$, for all $\epsilon > 0$, there are constants $0 < C', C < +\infty$ such that $$\label{eq:br-above-upper}
    C' n_k^{-2\alpha/({\underline{\Lambda}^\mathrm{g}}- \epsilon)}
    \leq 
    {\mathrm{Var}}_{{\mathbb{T}}_k}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}]
    \leq 
    C n_k^{-2\alpha/({\overline{\Lambda}^\mathrm{g}}+ \epsilon)},$$ and, in addition to , $$\exists\ \text{subsequence}\ (k'_j)_j,\ \text{s.t.}\ {\mathrm{Var}}_{{\mathbb{T}}_{k'_j}}[{\hat{\mu}_{\mathrm{ML}}}^{(k'_j)}]
    \leq \gamma n_{k'_j}^{-2\alpha/({\underline{\Lambda}^\mathrm{g}}+ \epsilon)}.$$

2.  ${\Lambda^\mathrm{b}}< 2\alpha$: then, for all $\epsilon > 0$, there are constants $0 < C', C < +\infty$ such that $$\label{eq:br-below-upper}
    C' n_k^{-2\alpha/({\underline{\Lambda}^\mathrm{g}}- \epsilon)}
    \leq 
    {\mathrm{Var}}_{{\mathbb{T}}_k}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}]
    \leq 
    C n_k^{-({\Lambda^\mathrm{b}}- \epsilon)/({\overline{\Lambda}^\mathrm{g}}+ \epsilon)},$$ where the lower bound in  above holds provided ${\underline{\Lambda}^\mathrm{g}}> 0$, and

    $$\exists\ \text{subsequence}\ (k'_j)_j,\ \text{s.t.}\ {\mathrm{Var}}_{{\mathbb{T}}_{k'_j}}[{\hat{\mu}_{\mathrm{ML}}}^{(k'_j)}]
    \leq \gamma n_{k'_j}^{-({\Lambda^\mathrm{b}}- \epsilon)/({\underline{\Lambda}^\mathrm{g}}+ \epsilon)}.$$

The following example shows that, when ${\Lambda^\mathrm{b}}< 2\alpha$, the upper bound in may not be achieved, but cannot be improved in general.

\[ex:twolevel\] Let $({\mathbb{T}}_k)_k$ be a tree sequence with two levels of nodes below the root: $D^{(k)}_0= e^{ \Lambda_0 \tau^{(k)}_0}$ nodes are attached to the root by edges of length $\tau^{(k)}_0 = \sigma T_k$, for some arbitrary choice of tree height $T_k\to\infty$ and $0<\sigma<1$. Each of these $D^{(k)}_0$ nodes has itself $D^{(k)}_1 = e^{ \Lambda_1 \tau^{(k)}_1}$ children along edges of length $\tau^{(k)}_1 = (1-\sigma) T_k$, and these form the leaves of ${\mathbb{T}}_k$.

\[prop:example2level\] For $0<\Lambda_0 < \Lambda_1$ and ${\mathcal{T}}=({\mathbb{T}}_k)_k$ described above, we have that ${\Lambda^\mathrm{b}}= \Lambda_0$, ${\Lambda^\mathrm{g}}= \sigma \Lambda_0 + (1-\sigma) \Lambda_1$, and $$\begin{aligned}
{\mathrm{Var}}_{{\mathbb{T}}_k}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}] &=& \gamma  e^{-2\alpha T_k}  + 
\gamma (1 - e^{-2\alpha \sigma T_k}) \, e^{-(\sigma\Lambda_0 + (1-\sigma)2\alpha) T_k} \nonumber\\
&&\qquad\qquad\qquad\qquad\qquad + \gamma (1 - e^{-2\alpha (1-\sigma) T_k}) \, e^{-{\Lambda^\mathrm{g}}T_k} .\label{eq:variance-2-level}\end{aligned}$$

This proposition is proved in Section \[sec:convergence-rate-mle\]. It implies that if $ 2\alpha \leq \Lambda_0 = {\Lambda^\mathrm{b}}$, the dominant term in the variance is $\gamma e^{-2\alpha T_k} = \gamma n_k^{-2\alpha /{\Lambda^\mathrm{g}}}$, as predicted by in Theorem \[thm:rate-supercritical\]. If instead $2\alpha \geq \Lambda_1 $, the dominant term in the variance is $\gamma e^{- {\Lambda^\mathrm{g}}T_k} = \gamma n_k^{-1},$ and we have $\sqrt{n_k}$-consistency. In the intermediate case when $\Lambda_0 < 2\alpha < \Lambda_1$, the dominant term in the variance is $ \gamma e^{- (\sigma \Lambda_0 + (1-\sigma) 2\alpha) T_k}
= \gamma n_k^{-\sigma ({\Lambda^\mathrm{b}}/{\Lambda^\mathrm{g}}) - (1-\sigma) (2\alpha/{\Lambda^\mathrm{g}})}.$ Therefore, depending on the value of $\sigma$, we can get the full range of exponent values between $-2\alpha/{\Lambda^\mathrm{g}}$ and $-{\Lambda^\mathrm{b}}/{\Lambda^\mathrm{g}}$, as given in .

In the other direction when ${\overline{\Lambda}^\mathrm{g}}< 2\alpha$, the picture is somewhat murkier. For example, by taking $\sigma$ close enough to $1$ in Example \[ex:twolevel\], it is possible to have ${\Lambda^\mathrm{g}}< 2\alpha$, yet not $\sqrt{n_k}$-consistency. The issue in Example \[ex:twolevel\] is the inhomogeneous growth rate. However, under extra regularity conditions, $\sqrt{n_k}$-consistency can be established. In words, the growth of the tree must be sufficiently homogeneous. In Theorem \[thm:subcritical\] below, we consider imposing the extra condition ${\Lambda^\mathrm{b}}= {\Lambda^\mathrm{g}}$, which does not hold in Example \[ex:twolevel\].

\[thm:subcritical\] Let $({\mathbb{T}}_k)_k$ be a tree sequence with ${\overline{\Lambda}^\mathrm{g}}< 2\alpha$. Then $${\mathrm{Var}}_{{\mathbb{T}}_k}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}] = \Omega\left(n_k^{-1}\right).$$ Further if:

1.  ${\Lambda^\mathrm{b}}= {\overline{\Lambda}^\mathrm{g}}> 0$ then, for all $\epsilon > 0$, $
    {\mathrm{Var}}_{{\mathbb{T}}_k}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}] = O\left(n_k^{-(1-\epsilon)}\right).
    $

2.  ${\Lambda^\mathrm{ug}}< 2\alpha$ then $
    {\mathrm{Var}}_{{\mathbb{T}}_k}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}] = O\left(n_k^{-1}\right).
    $

Theorems \[thm:rate-supercritical\] and \[thm:subcritical\] are proved in Section \[sec:convergence-rate-mle\]. All our results on the estimation of $\mu$ leverage a useful characterization of the variance of linear estimators in terms of electrical networks. An analogous characterization is used in ancestral state reconstruction [@peres1999climb]. Note that our results are not as clean as those obtained for ancestral state reconstruction. As Example \[ex:twolevel\] showed, estimation of $\mu$ is somewhat sensitive to the “homogeneity” of the growth. In Section \[sec:sensitivity-to-alpha\], we show that assuming $\alpha$ is known is inconsequential, provided a good estimate of $\alpha$ is available. Such an estimate is discussed next.

#### Results concerning the parameters $\alpha$ and $\gamma$

Our main result for $\alpha$ and $\gamma$ is a $\sqrt{n_k}$-consistent estimator under the following assumption: there are two separate “bands” of node ages, each containing a number of internal nodes growing linearly with the number of leaves.

\[assump:alpha\] Define $n_k(c,c')$ as the number of nodes in ${\mathbb{T}}_k$ of age (height from the leaves) in $(c,c')$. Assume that there are constants $\beta >0$ and $0 < c_1 < c_1' < c_2 < c_2' < \infty$ such that $n_k(c_i,c_i') \geq \beta n_k, i =1,2,$ for all $k$ large enough.

As shown in Corollary \[cor:yule\_alpha\_gamma\], this assumption holds for the Yule process, a speciation model frequently used in practice.

\[thm:main-alpha-rate\] Let $({\mathbb{T}}_k)_k$ be a sequence of ultrametric trees satisfying Assumptions \[assumption:unbounded\] and \[assump:alpha\]. Then there is an estimator $(\hat\alpha_k,\hat\gamma_k)_k$ of $(\alpha,\gamma)$ such that $|\hat\alpha_k -\alpha| = O_p(n_k^{-1/2})$ and $|\hat\gamma_k -\gamma| = O_p(n_k^{-1/2})$.

The proof, found in Section \[sec:alpha-rate\], is based on the common notion of contrasts. Assumption \[assump:alpha\] ensures the existence of an appropriate set of such contrasts. The key point is that this extra assumption can be satisfied no matter what the growth and branching number are, indicating that the estimation of $\alpha$ and $\gamma$ is unaffected by the growth of the tree unlike $\mu$. Intuitively, $\mu$ is a more “global” parameter.

Special cases {#sec:special-cases-mu}
-------------

We apply here the results stated in Section \[sec:statementresults\] to a number of scenarios. The tree of life naturally gives rise to two types of tree sequences. If one imagines sampling an increasing number of contemporary species, one obtains a nested sequence, defined as follows.

A sequence of trees ${\mathcal{T}}=({\mathbb{T}}_k)_{k}$ is [*nested*]{} if, for all $k$, $n_k = k$ and ${\mathbb{T}}_k$ restricted to the first $k-1$ species is identical to ${\mathbb{T}}_{k-1}$ as an ultrametric.

An example of nested trees is given by a caterpillar sequence.

\[ex:caterpillar\] Let $(t_k)_k$ be a sequence of nonnegative numbers such that $\limsup_k t_k = +\infty$. Let ${\mathbb{T}}_1$ be a one-leaf star with height $T_1 = t_1$. For $k > 1$, let ${\mathbb{T}}_k$ be the caterpillar-like tree obtained by adding a leaf edge with leaf $k$ to ${\mathbb{T}}_{k-1}$ at height $t_k$ on the path between leaf $1$ and the root of ${\mathbb{T}}_{k-1}$, if $t_k \leq T_{k-1}$. If instead $t_k > T_{k-1}$, create a new root at height $t_k$ with an edge attached to the root of ${\mathbb{T}}_{k-1}$ and an edge attached to $k$ (see Figure \[fig:caterpillar\]).

\[cor:consistency-nested\] Let ${\mathcal{T}}$ be a nested sequence such that the height $T_k$ goes to infinity. Then ${\mathcal{T}}$ satisfies Assumption \[assumption:unbounded\] and the MLE for $\mu$ is consistent on ${\mathcal{T}}$.

Let $k_j$ be the subsequence such that $T_{k_{j+1}} > T_{k_j}$ for every $j$ and $T_i = T_{k_j}$ for all $i = k_j + 1, \ldots, k_{j+1}-1$. Then, for all $s \in (0,+\infty)$, as $k$ goes to $+\infty$ $\pi^k_s$ eventually contains all leaves $k_j$ such that $T_{k_j} \geq s$. Since $T_k \to +\infty$, the result follows.

If one is modeling the growth of the tree of life in time, instead of modeling increased sampling of contemporary species, one obtains a growing sequence as follows. Let ${\mathbb{T}}_0$ be a rooted infinite tree of bounded degree, with branch lengths and no leaves. Think of the branches of ${\mathbb{T}}_0$ as a continuum of [*points*]{} whose distance from the endpoints grows linearly. Then, for $t \geq 0$, we define $\mathcal{B}_t({\mathbb{T}}_0)$ as the tree made of the set of points of ${\mathbb{T}}_0$ at distance at most $t$ from the root.

A sequence of trees $({\mathbb{T}}_k)_k$ is a [*growing sequence*]{} of trees if there is an infinite tree ${\mathbb{T}}_0$ as above and an increasing sequence of non-negative reals $(t_k)_k$ such that ${\mathbb{T}}_k$ is isomorphic to $\mathcal{B}_{t_k}({\mathbb{T}}_0)$ as an ultrametric.

Let $({\mathbb{T}}_k)_k$ be a growing sequence such that the height $T_k = t_k$ goes to infinity. Then ${\mathcal{T}}$ satisfies Assumption \[assumption:unbounded\] and the MLE for $\mu$ is consistent on $({\mathbb{T}}_k)_k$.

Fix $s \in (0,+\infty)$. For $L =1,2,\ldots$, let $k'_L$ be the smallest $k$ such that $n_k \geq L$ and let $k''_L$ be the smallest $k > k'_L$ such that $T_{k} \geq T_{k'_L} + s$. Then, for all $k \geq k''_L$, $|\pi^k_s|\geq L$. Letting $L$ go to $+\infty$ gives the result.

\[ex:yule\] Let ${\mathbb{T}}_0$ be a tree generated by a pure-birth (Yule) process with rate $\lambda > 0$: starting with one lineage, each current lineage splits independently after an exponential time with mean $\lambda^{-1}$ (see e.g. [@semple2003phylogenetics]). For any (possibly random) sequence of increasing non-negative reals $(t_k)_k$ with $t_k \to +\infty$, $\mathcal{B}_{t_k}({\mathbb{T}}_0)$ (that is, ${\mathbb{T}}_0$ run up to time $t_k$), forms a growing sequence.

The following result is proved in Section \[section:proofs-special\].

\[cor:mu-yule\] Let $({\mathbb{T}}_k)_k$ be a Yule sequence with rate $0 < \lambda < +\infty$. Then, with probability $1$ (on the generation of ${\mathbb{T}}_0$),

1.  $({\hat{\mu}_{\mathrm{ML}}}^{(k)})_k$ is consistent.

2.  If $\lambda < 2 \alpha$, $({\hat{\mu}_{\mathrm{ML}}}^{(k)})_k$ is $\sqrt{n_k}$-consistent.

3.  If $\lambda > 2 \alpha$, $({\hat{\mu}_{\mathrm{ML}}}^{(k)})_k$ is [*not*]{} $\sqrt{n_k}$-consistent and for all $\epsilon > 0$ there is $0 < C', C < +\infty$ such that $$C' n_k^{-2\alpha\lambda^{-1}-\epsilon}
    \leq
    {\mathrm{Var}}_{{\mathbb{T}}_k}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}]
    \leq
    C n_k^{-2\alpha\lambda^{-1}+\epsilon}.$$

\[cor:yule\]

We also apply the estimators $\hat\alpha$ and $\hat\gamma$ to the Yule model. For simplicity, we take the sequence of times at which new speciation events occur (although this assumption is not crucial). For $k \geq 1$, let $t_k$ be the first time at which ${\mathbb{T}}_0$ has $k+1$ lineages. Then $n_k = k$ for all $k$ and $t_k \to +\infty$ so that Assumption \[assumption:unbounded\] is satisfied. The following result is proved in Section \[section:proofs-special\].

\[cor:yule\_alpha\_gamma\] Let $({\mathbb{T}}_k)_k$ be a Yule sequence with $n_k=k$ as above. Then Assumption \[assump:alpha\] is satisfied asymptotically, and hence $|\hat\alpha_k -\alpha| = O_p(n_k^{-1/2})$ and $|\hat\gamma_k -\gamma| = O_p(n_k^{-1/2})$.

Application to experimental design for trait evolution studies {#sec:app}
==============================================================

Thanks to recent developments in technology, scientists have reconstructed several large phylogenetic trees with thousands of species such as trees containing 4507 mammal species [@binindaEmonds-etal07] and 9993 bird species [@jetz12]. However, researchers may not be able to collect trait data from all species, due to limited resources and funding. Thus, many studies are only based on a subset of species in the available tree. For example, to study the evolution of body size in mammals, Cooper and Purvis [@cooperPurvis10] used 3473 of the 4507 species in their tree, and Venditti et al. [@vendittiMeadePagel11] incorporated 3185 species in their analysis. When considering extra data collection, an important question arises: can additional species increase the precision of our estimates? Our theoretical results help answering this question for the OU tree model:

1.  If $\hat \lambda \ll 2 \hat \alpha$, additional species tend to be very informative for estimating $\mu$ (Corollary \[cor:mu-yule\]).

2.  If $\hat \lambda \gg 2 \hat \alpha$, additional species that do not increase tree height tend to be non-informative for estimating $\mu$ (Corollary \[cor:mu-yule\]).

3.  When $\hat \lambda$ is around $2 \hat \alpha$, it is not clear whether additional species are informative for estimating $\mu$.

4.  Additional species tend to be informative for estimating $\alpha$ and $\gamma$ (Corollary \[cor:yule\_alpha\_gamma\]).

#### Example:   {#example .unnumbered}

In [@vendittiMeadePagel11], body size evolution was studied using 3185 mammal species. Would it be worth the effort to collect data for the remaining 1322 species in the tree, to increase the precision of estimating $\mu$? To answer this question about sampling utility, we first need to estimate the speciation rate $\lambda$ and the selection strength $\alpha$. The 4507-species mammal tree was rescaled to have height $1$ and its speciation rate was estimated to be $11.83$ using maximum likelihood (`yule` function in the R package `ape` [@ape]). We also estimated $\hat\alpha=0.01$ using maximum likelihood (`phylolm` function in the R package `phylolm` [@HoAne14sb]). Note that the tree formed by the 3185 species has the same height as the full tree with all 4507 species. Since $\hat \lambda \approx 11.83 \gg 0.02 \approx 2 \hat \alpha$, additional species tend to be non-informative and our recommendation is to stop data collection. Our conclusion is consistent with simulations in [@HoAne13; @ho2014intrinsic], which showed that additional species are non-informative for estimating $\mu$ if they do not increase tree height, when $\alpha$ is low. Our recommendation here specifies the critical value of $\alpha$ below which additional sampling is of little utility.

To further demonstrate the relationship between sampling utility and $\alpha$ (or $\lambda$) at fixed tree height, we simulated data according to the OU model along the 4507-species mammal tree with $\mu = 0$, $\gamma = 1$, and several values of $\alpha$ ranging from $0.01$ to $300$. For every set of parameters, we simulated 2000 data sets using the `rTrait` function (R package `phylolm`). Then, ${\hat{\mu}_{\mathrm{ML}}}$ was computed for each data set using the `phylolm` function. The sample variance of ${\hat{\mu}_{\mathrm{ML}}}$ (Figure \[fig:fig3\]) was found to be about $e^{-2 \alpha}$ when $2 \alpha \ll \hat \lambda$, and about $1/n = 1/4507$ when $2 \alpha \gg \hat \lambda$.

![Sample variance of ${\hat{\mu}_{\mathrm{ML}}}$ (red points) as a function of $\alpha$, from the simulation on the mammal tree. When $\alpha$ is small, $\text{var}({\hat{\mu}_{\mathrm{ML}}}) \approx e^{-2 \alpha}$ (purple line). When $\alpha$ is large, $\text{var}({\hat{\mu}_{\mathrm{ML}}}) \approx 1/4507$ (blue line).[]{data-label="fig:fig3"}](change_alpha)

To illustrate the relationship between sampling utility and $\alpha$ (or $\lambda$) when the tree height varies, we simulated 400 trees under the Yule process using the `sim.bdtree` function (R package `geiger` [@geiger]). We used speciation rate $\lambda = 11.83$, which was the maximum likelihood estimate from the mammal tree. The tree height was varied from $0.05$ to $1$ and we simulated $20$ trees for each tree height. We calculated $\text{var}({\hat{\mu}_{\mathrm{ML}}})$ corresponding to three fixed values of $\alpha$ $(0.1,\lambda/2,30)$ using and the `three.point.compute` function (R package `phylolm`). The results showed that (Figure \[fig:fig4\]) when $2 \alpha \ll \lambda$, $e^{-2 \alpha T}$ approximates $\text{var}({\hat{\mu}_{\mathrm{ML}}})$ better than $1/n$. On the other hand, when $2 \alpha \gg \lambda$, $1/n$ is a better approximation.

![Variance of ${\hat{\mu}_{\mathrm{ML}}}$ on ramdom trees, simulated under the Yule process. Top: $\text{var}({\hat{\mu}_{\mathrm{ML}}})$ against number of leaves $n$. Bottom: $\text{var}({\hat{\mu}_{\mathrm{ML}}})$ against tree height $T$. The true value of $\alpha$ is either small (0.1, left), or $\lambda/2$ (5.9, middle), or large (30, right).[]{data-label="fig:fig4"}](yule)

Taken together, our results show that when $2 \alpha \ll \hat \lambda$, the variance of ${\hat{\mu}_{\mathrm{ML}}}$ depends on the tree height, not the sample size. So, additional sampling that does not increase tree height is not recommended. On the other hand, when $2 \alpha \gg \hat \lambda$ the variance of ${\hat{\mu}_{\mathrm{ML}}}$ is of order $1/n$, as if we had $n$ independent samples. In this case additional species are very informative, and additional sampling is recommended if affordable.

Proofs of results for estimating $\mu$ {#sec:proofs}
======================================

We develop here necessary tools (Section \[sec:mu-tools\]), then prove Theorem \[thm:criterion-consistency\] (Section \[sec:consistency-mle\]), Theorems \[thm:rate-supercritical\] and \[thm:subcritical\] (Section \[sec:convergence-rate-mle\]), which assume that $\alpha$ is known. Using arguments from the proofs, we also identify examples showing that the sample mean $\overline Y$ can perform significantly worse than ${\hat{\mu}_{\mathrm{ML}}}$, and we show that Assumption \[assumption:unbounded\] is not sufficient in Theorem \[thm:criterion-consistency\] for the consistency of ${\hat{\mu}_{\mathrm{ML}}}$. We prove an alternative sufficient condition based on the branching number (Proposition \[thm:consistency-branching\] below). In Section \[section:proofs-special\], we prove Corollaries \[cor:mu-yule\] and \[cor:yule\_alpha\_gamma\]. Finally, in Section \[sec:sensitivity-to-alpha\] we discuss the sensitivity of the MLE to estimation errors on $\alpha$.

Bounding the variance of the MLE {#sec:mu-tools}
--------------------------------

Fix an ultrametric species tree ${\mathbb{T}}$ with leaf set ${\mathscr{L}}$, number of tips $n = |{\mathscr{L}}|$, and root $\rho$. We also fix $\alpha > 0$.

#### A formula for the variance

Let ${{\boldsymbol\theta}}= (\theta_\ell)_{\ell\in{\mathscr{L}}}$, with ${{\boldsymbol\theta}}' {\mathbf{1}}= 1$ and $\theta_\ell \in [0,1]$ for all $\ell$, and recall that $
Y_{{\boldsymbol\theta}}= \sum_{\ell \in {\mathscr{L}}} \theta_\ell Y_\ell
$ is an unbiased estimator of $\mu$. By defining, for each branch $b$, $$\label{eq:flow-on-branch}
\theta_b = \sum_{\ell \in \mathscr{L}} {\mathbbm{1}}_{b \in p(\rho,\ell)}\theta_\ell,$$ where $p(\rho,\ell)$ is the path from $\rho$ to $\ell$, we naturally associate to the coefficients ${{\boldsymbol\theta}}$ a flow on the edges of ${\mathbb{T}}$, defined as follows.

A [*flow*]{} ${{\boldsymbol\eta}}$ is a mapping from the set of edges to the set of positive numbers such that, for every edge $b$, we have $
\eta_b = \sum_{b' \in O_b}{\eta_{b'}}
$ where $O_b$ is the set of outgoing edges stemming from $b$ (with the edges oriented away from the root). Define $\|{{\boldsymbol\eta}}\| = \sum_{b \in O_\rho}{\eta_b}$. We say that ${{\boldsymbol\eta}}$ is a [*unit flow*]{} if $\|{{\boldsymbol\eta}}\| = 1$. We extend ${{\boldsymbol\eta}}$ to vertices $v$ in ${\mathbb{T}}$ by defining $\eta_v$ as the flow on the edge entering $v$. Similarly, for a point $x$ in ${\mathbb{T}}$, we let $\eta_x$ be the flow on the corresponding edge or vertex.

For every edge $b$ of ${\mathbb{T}}$, we set $
R_b = (1 - e^{-2 \alpha |b|}) e^{2 \alpha \delta(\rho,b)}
$ where $|b|$ is the length of $b$ and $\delta(\rho,b)$ is the length of the path from the root to $b$ (inclusive).

\[prop:key\] Let $\mathcal{F}$ be the set of unit flows from $\rho$ to ${\mathscr{L}}$. Let $E$ be the set of edges and $T$ be the height of tree ${\mathbb{T}}$. For any ${{\boldsymbol\theta}}\in \mathcal{F}$, we have $$\label{eq:variance-formula}
{\mathrm{Var}}[Y_{{\boldsymbol\theta}}] = \gamma e^{- 2 \alpha T} \left( 1 + \sum_{b \in E}{R_b\theta^2_b} \right)$$ so that ${\mathrm{Var}}[{\hat{\mu}_{\mathrm{ML}}}] = 
\inf_{{{\boldsymbol\theta}}\in \mathcal{F}}\  
\gamma e^{- 2 \alpha T} ( 1 + 
\sum_{b \in E}{R_b\theta^2_b}).
$

As detailed in [@peres1999climb], a species tree can be interpreted as an electrical network with resistance $R_b$ on edge $b$. The minimum $\mathcal{R}_{\mathbb{T}}$ of $\sum_{b \in E}{R_b\theta^2_b}$ over unit flows (corresponding to the MLE) is known as the [*effective resistance*]{} of ${\mathbb{T}}$, which can be interpreted in terms of a random walk on the tree. See [@peres1999climb] for details.

The second part follows from the first, because ${\hat{\mu}_{\mathrm{ML}}}$ is the best unbiased linear estimator of $\mu$. The proof of follows from a computation in [@EvKePeSc:00 Lemma 5.1]. For every node $u$ of the tree, by a telescoping argument, $$\label{eq:telescope}
e^{2 \alpha \delta(\rho,u)} - 1 
= \sum_{b \in p(\rho,u)}{R_b}$$ where $\delta(\rho,u)$ is the distance from $\rho$ to $u$, and $p(\rho,u)$ is the path from $\rho$ to $u$. Denote by $v \wedge w$ the most recent common ancestor of $v$ and $w$. Then $$\begin{aligned}
{\mathrm{Var}}[Y_{{\boldsymbol\theta}}] &=& \gamma  \sum_{v,w \in \mathscr{L}}{\theta_v\theta_w \frac{e^{-2 \alpha T}}{e^{-2 \alpha \delta(\rho, v \wedge w)}}}
= \gamma e^{- 2 \alpha T} \sum_{v,w \in \mathscr{L}}{\theta_v\theta_w\bigg( 1 + \sum_{b \in p(\rho,v \wedge w)}{R_b} \bigg)}\\
&=& \gamma e^{- 2 \alpha T} \bigg( 1 + \sum_{b \in E}{R_b \sum_{v,w \in \mathscr{L}}{{\mathbbm{1}}_{b \in p(\rho,v \wedge w)}\theta_v\theta_w}} \bigg)\\
&=& \gamma e^{- 2 \alpha T} \left [ 1 + \sum_{b \in E}{R_b \Big( \sum_{v \in \mathscr{L}}{{\mathbbm{1}}_{b \in p(\rho,v)}\theta_v} \Big) \Big( \sum_{w \in \mathscr{L}}{{\mathbbm{1}}_{b \in p(\rho,w)}\theta_w} \Big)} \right ]\\
&=& \gamma e^{- 2 \alpha T} \bigg( 1 + \sum_{b \in E}{R_b\theta^2_b} \bigg),\end{aligned}$$ where the second equality follows from , the fourth equality follows from ${\mathbbm{1}}_{b \in p(\rho,v \wedge w)} 
= {\mathbbm{1}}_{b \in p(\rho,v)} {\mathbbm{1}}_{b \in p(\rho,w)}$, and the last equality follows from .

For $0 \leq t \leq T$, let $\pi_t$ be the set of points at distance $t$ from the root (that is, the cutset corresponding to time $t$ away from the root). Noting that $$R_b = 2\alpha \int_{\delta(\rho,b)-|b|}^{\delta(\rho,b)}
e^{2\alpha s} \mathrm{d} s,$$ we get the following convenient formula:

\[cor:integral-form\] For any unit flow ${{\boldsymbol\theta}}$ from $\rho$ to $\mathscr{L}$, we have $${\mathrm{Var}}[Y_{{\boldsymbol\theta}}] = \gamma e^{- 2 \alpha T} \left [ 1 + 
2\alpha \int_0^T e^{2\alpha s} 
\Big(\sum_{x\in \pi_s} \theta_x^2\Big)
\mathrm{d}s \right ].$$

As a first important application of Proposition \[prop:key\] and Corollary \[cor:integral-form\], we show that the variance of the MLE of $\mu$ can be controlled by the branching number. The result is characterized by a transition at ${\Lambda^\mathrm{b}}=2\alpha$, similarly to Example \[ex:first-transition-continued\].

\[prop:variance-branching\] Let ${\mathcal{T}}= ({\mathbb{T}}_k)_k$ be a tree sequence with branching number ${\Lambda^\mathrm{b}}> 0$. Then, for all $\Lambda < {\Lambda^\mathrm{b}}$, there is ${\mathcal{I}}_\Lambda$ such that $${\mathrm{Var}}_{{\mathbb{T}}_k}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}]
\leq 
\begin{cases}
\gamma\Big(1 + \frac{2\alpha}{{\mathcal{I}}_\Lambda(2\alpha-\Lambda)}\Big)
e^{- \Lambda T_k}, &\text{if $\Lambda < 2\alpha$,}\\
\gamma\;\Big(1 + \frac{2\alpha T_k}{{\mathcal{I}}_\Lambda}\Big)
e^{-2\alpha T_k}, &\text{if $\Lambda = 2\alpha$,}\\
\gamma\Big(1 + \frac{2\alpha}{{\mathcal{I}}_\Lambda(\Lambda-2\alpha)}\Big)
e^{-2\alpha T_k}, &\text{if $\Lambda > 2\alpha$.}
\end{cases}$$

For $\Lambda < {\Lambda^\mathrm{b}}$, let $
{\mathcal{I}}_\Lambda = \inf_{k, \pi \in \Pi^k} 
\sum_{x \in \pi}{ e^{- \Lambda \delta_k(\rho,x)}  > 0 }.
$ By the max-flow min-cut theorem (see e.g. [@lawler1976combinatorial]), there is a flow ${{\boldsymbol\eta}}^{(k)}$ on ${\mathbb{T}}_k$ with $$\label{eq:max-flow-value}
\|{{\boldsymbol\eta}}^{(k)}\| \geq {\mathcal{I}}_\Lambda$$ and $$\label{eq:max-flow-capacity}
\eta^{(k)}_x \leq e^{- \Lambda \delta_k(\rho,x)},$$ for all points $x$ in ${\mathbb{T}}_k$. Normalize ${{\boldsymbol\eta}}^{(k)}$ as ${{\boldsymbol\theta}}^{(k)} = {{\boldsymbol\eta}}^{(k)}/\|{{\boldsymbol\eta}}^{(k)}\|$. By Proposition \[prop:key\] and Corollary \[cor:integral-form\], for $\Lambda \neq 2\alpha$, $$\begin{aligned}
{\mathrm{Var}}_{{\mathbb{T}}_k}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}]
&\leq& 
\gamma e^{- 2 \alpha T_k} \left [ 1 + 
2\alpha \int_0^{T_k} e^{2\alpha s} 
\Big(\sum_{x\in \pi^k_s} (\theta^{(k)}_x)^2\Big)
\mathrm{d}s \right ]\\
&\leq& 
\gamma e^{- 2 \alpha T_k} \left [ 1 + 
2\alpha \int_0^{T_k} e^{2\alpha s} 
\Big(\sum_{x\in \pi^k_s} \theta^{(k)}_x \frac{e^{- \Lambda \delta_k(\rho,x)}}{{\mathcal{I}}_\Lambda}\Big)
\mathrm{d}s \right ]\\
&\leq& 
\gamma e^{- 2 \alpha T_k} \left [ 1 + 
\frac{2\alpha}{{\mathcal{I}}_\Lambda}
\int_0^{T_k} e^{( 2\alpha -\Lambda)s} 
\mathrm{d}s \right ]\\
&=& 
\gamma \left [ e^{- 2 \alpha T_k} + 
\frac{2\alpha}{{\mathcal{I}}_\Lambda (2\alpha-\Lambda)}\big(e^{- \Lambda T_k} - e^{- 2 \alpha T_k}\big) 
 \right ].\end{aligned}$$ where the second inequality follows from  and , and the third inequality follows from the fact that $\delta_k(\rho,x) = s$ for $x \in \pi^k_s$ by definition and that $\sum_{x\in \pi^k_s} \theta^{(k)}_x = 1$. Similarly if $\Lambda = 2\alpha$ $${\mathrm{Var}}_{{\mathbb{T}}_k}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}]
\leq \gamma \left[ e^{- 2 \alpha T_k} + 
\frac{2\alpha e^{- 2 \alpha T_k}  T_k}{{\mathcal{I}}_\Lambda}
\right].$$

#### Removing bottlenecks

Examining , one sees that a natural bound on ${\mathrm{Var}}[Y_{{\boldsymbol\theta}}]$ is obtained by “splitting an edge” in ${\mathbb{T}}$.

\[def:splitting\] Let ${\mathbb{T}}$ be an ultrametric tree with edge set $E$. Let $b_0 = (x_0,y_0)$ be a branch in ${\mathbb{T}}$ (where $x_0$ is closer to the root) and let $b_i = (y_0,y_i)$, $i = 1,\ldots, D$, be the outgoing edges at $y_0$. The operation of [*splitting branch $b_0$*]{} to obtain a new tree ${\mathbb{T}}'$ with edge set $E'$ is defined as follows: remove $b_0, b_1,\ldots, b_D$ from ${\mathbb{T}}$; add $D$ new edges $b'_i = (x_0,y_i)$ of length $|b_0| + |b_i|$, $i=1,\ldots,D$ (see Figure \[fig:fig2\]). We call [*merging*]{} the opposite operation of undoing the above splitting.

Note that the number of tips in ${\mathbb{T}}$ and ${\mathbb{T}}'$ above are the same, and therefore we can use the same estimator $Y_{{\boldsymbol\theta}}$ on both of them.

\[lemma:splitting\] Let ${\mathbb{T}}$ be an ultrametric tree, let $b_0$ be a branch in ${\mathbb{T}}$, and let ${\mathbb{T}}'$ be obtained from ${\mathbb{T}}$ by splitting $b_0$. Then for any nonnegative ${{\boldsymbol\theta}}= (\theta_\ell)_{\ell\in{\mathscr{L}}}$ $${\mathrm{Var}}_{{\mathbb{T}}'}[Y_{{\boldsymbol\theta}}] \leq {\mathrm{Var}}_{{\mathbb{T}}}[Y_{{\boldsymbol\theta}}].$$

We use the notation of Definition \[def:splitting\]. Denote by $(\theta_b)_{b\in E}$ and $(\theta'_{b})_{b\in E'}$ the flows associated to ${{\boldsymbol\theta}}$ by  on ${\mathbb{T}}$ and ${\mathbb{T}}'$ respectively. For any branch $b$, except $b_0,b_1,\ldots,b_D$ and $b'_1,\ldots,b'_D$, we have $\theta_b = \theta'_b$, as the descendant leaves of $b$ on ${\mathbb{T}}$ and ${\mathbb{T}}'$ are the same. Think of $b'_i = (x_0,y_i)$, $i=1,\ldots,D$, as being made of two consecutive edges $b''_i = (x_0,y_i')$ and $b'''_i = (y'_i,y_i)$ with $|b''_i| = |b_0|$ and $|b'''_i| = |b_i|$ (and note, for sanity check, that $R_{b'_i} = R_{b''_i}  + R_{b'''_i}$). Then, $\theta_{b_i} = \theta_{b'''_i}$ and $R_{b_i} = R_{b'''_i}$, and by  $$\frac{{\mathrm{Var}}_{{\mathbb{T}}}[Y_{{\boldsymbol\theta}}] - {\mathrm{Var}}_{{\mathbb{T}}'}[Y_{{\boldsymbol\theta}}]}{\gamma e^{-2 \alpha T}}
= R_{b_0} \theta_{b_0}^2 - \sum_{i=1}^D R_{b''_i} \theta_{b''_i}^2
= R_{b_0} \Big(\sum_{i=1}^D \theta_{b''_i} \Big)^2 - 
R_{b_0} \sum_{i=1}^D  \theta_{b''_i}^2 \geq 0,$$ where we used that $R_{b_0} = R_{b''_i}$ and the nonnegativity of the $\theta_{b''_i}$’s.

Comparing ${\mathbb{T}}$ to a star we then get:

\[prop:variance-lower\] Let ${\mathbb{T}}$ be an ultrametric tree with $n$ tips and height $T$. Then $${\mathrm{Var}}_{\mathbb{T}}[{\hat{\mu}_{\mathrm{ML}}}] \geq  \gamma
\left(e^{-2\alpha T} + \frac{1 - e^{-2\alpha T}}{n}\right).$$

Split all edges in ${\mathbb{T}}$ by repeatedly applying Lemma \[lemma:splitting\] until a star tree with $n$ leaves and height $T$ is obtained. The result then follows from .

The following example will be useful when proceeding in reverse, to find an upper bound on the variance of ${\hat{\mu}_{\mathrm{ML}}}$.

\[ex:spherical\] Let ${\mathbb{T}}$ be a spherically symmetric, ultrametric tree, that is, a tree such that all vertices at the same graph distance from the root have the same number of outgoing edges, all of the same length. Let $D_h$, $h=0,\ldots,H-1$, be the out-degree of vertices at graph distance $h$ (where $h=0$ and $h=H$ correspond to the root and leaves respectively) and let $\tau_h$ be the corresponding branch length. Notice that $\beta_1^2+\cdots+\beta_d^2$, subject to $\beta_1+\cdots+\beta_d = 1$, is minimized at $\beta_1=\cdots=\beta_d = 1/d$. Hence, since ${\hat{\mu}_{\mathrm{ML}}}$ is the best unbiased linear estimator and arguing inductively from the leaves in , we see that ${\hat{\mu}_{\mathrm{ML}}}=\overline{Y}$ in this case. The mean squared error is, by , $$\begin{aligned}
{\mathrm{Var}}[{\hat{\mu}_{\mathrm{ML}}}] 
&=& \gamma e^{-2\alpha T} \left[
1 + \sum_{h=0}^{H-1} \left(\prod_{h'=0}^{h} D_{h'}\right) (1 - e^{-2\alpha \tau_h})
e^{2\alpha \sum_{h'=0}^{h} \tau_{h'}} 
\prod_{h'=0}^{h}\frac{1}{D_{h'}^2} 
\right]\nonumber\\
&=& \gamma e^{-2\alpha T} \left[
1 + \sum_{h=0}^{H-1} (1 - e^{-2\alpha \tau_h})
 \prod_{h'=0}^{h}\frac{e^{2\alpha \tau_{h'}}}{D_{h'}} 
\right].\label{eq:spherical}\end{aligned}$$

\[prop:variance-upper\] Let ${\mathbb{T}}$ be an ultrametric tree with height $T$. Recall that $\pi_t$ be the set of points at distance $t$ from the root. Then $${\mathrm{Var}}_{\mathbb{T}}[{\hat{\mu}_{\mathrm{ML}}}]
\leq
\inf_{0 \leq t \leq T}
\gamma 
\left(
e^{-2\alpha (T-t)} 
+ \frac{1 - e^{-2\alpha(T-t)}}{|\pi_t|}
\right).$$

Let $0 \leq t \leq T$. For all points $x$ in $\pi_t$, choose one descendant leaf $\ell_x$ of $x$ and define ${{\boldsymbol\theta}}$ as $$\theta_\ell
=
\begin{cases}
1/|\pi_t| & \text{if $\ell = \ell_x$ for some $x$},\\
0         & \text{otherwise}.
\end{cases}$$ Divide all branches crossing $\pi_t$ into two branches meeting at $\pi_t$. Then merge all branches above $\pi_t$ (that is, closer to the root) by repeatedly applying Lemma \[lemma:splitting\]. By , removing all branches $b$ with $\theta_b = 0$ does not affect the variance, and from Example \[ex:spherical\] with $H = 2$, $D_0 = 1$, $D_1 = |\pi_t|$, $\tau_0 = t$, and $\tau_1 = T-t$, we get $$\begin{aligned}
{\mathrm{Var}}_{\mathbb{T}}[{\hat{\mu}_{\mathrm{ML}}}]
&\leq& 
\gamma e^{-2\alpha T}
\left[
1+ (1 - e^{-2\alpha t})e^{2\alpha t} 
+ (1 - e^{-2\alpha(T-t)}) e^{2\alpha t} \frac{e^{2\alpha (T-t)}}{|\pi_t|}
\right]\\
&\leq& 
\gamma 
\left[
e^{-2\alpha (T-t)} 
+ \frac{1 - e^{-2\alpha(T-t)}}{|\pi_t|}
\right].\end{aligned}$$

#### The two estimators ${\hat{\mu}_{\mathrm{ML}}}$ vs. $\overline{Y}$

As an application of the previous proposition, we provide an example where ${\hat{\mu}_{\mathrm{ML}}}$ performs significantly better than $\overline{Y}$. Roughly, the example shows that $\overline{Y}$ can perform poorly on asymmetric trees.

\[ex:caterpillar-bad\] Consider a caterpillar sequence $(\mathbb{T}_k)_{k}$, as defined in Example \[ex:caterpillar\], with $t_{2m+1} = m$ and $t_{2m} = 1$ for all $m$, as shown in Figure \[fig:fig1\]. Note that the tree height is $T_{2m+1} = T_{2m+2} = m$ and the cut sets $\pi_t^{k}$ of ${\mathbb{T}}_k$ at time $t$ satisfy $|\pi_{m-1}^{2m+1}| = |\pi_{m-1}^{2m+2}| = m$. Therefore, by Proposition \[prop:variance-upper\], $$\max\left\{{\mathrm{Var}}_{{\mathbb{T}}_{2m+1}}[{\hat{\mu}_{\mathrm{ML}}}],
{\mathrm{Var}}_{{\mathbb{T}}_{2m+2}}[{\hat{\mu}_{\mathrm{ML}}}]\right\}
\leq \gamma \left[
e^{-2\alpha (m-1)} + \frac{1}{m} 
\right] \to 0,$$ as $m \to +\infty$, and hence ${\hat{\mu}_{\mathrm{ML}}}$ is consistent. On the other hand, note that ${\mathrm{Cov}}[Y_i,Y_j] \geq 0$ for all pairs of leaves $i, j$ in ${\mathbb{T}}_k$. Therefore, $$\begin{aligned}
{\mathrm{Var}}_{{\mathbb{T}}_{2m}}\left[\overline Y\right] 
&=& \frac{1}{4m^2} {\mathrm{Var}}\left[
\sum_{\ell=1}^{2m}{Y_\ell} \right] \geq \frac{1}{4m^2} {\mathrm{Var}}\left[
\sum_{i=1}^m{Y_{2i}} \right] = \frac{1}{4m^2} \sum_{i,j=1}^m 
{\mathrm{Cov}}[Y_{2i},Y_{2j}]\\ 
&\geq& \frac{1}{4m^2} m^2 \gamma e^{-2 \alpha}
= \frac{\gamma e^{- 2 \alpha}}{4}.\end{aligned}$$ So, $\overline Y$ is not consistent.

Proof of Theorem \[thm:criterion-consistency\] and Sufficiency of Conditions {#sec:consistency-mle}
----------------------------------------------------------------------------

First assume . From Proposition \[prop:variance-upper\], for all $s$, $$\limsup_k {\mathrm{Var}}_{{\mathbb{T}}_k}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}]
\leq
\limsup_k \gamma 
\left[
e^{-2\alpha s} 
+ \frac{1 - e^{-2\alpha s}}{\left|\tilde{\pi}^k_s\right|}
\right]
\leq \gamma e^{-2\alpha s}.$$ Taking $s$ to $+\infty$ gives consistency. On the other hand, assume by contradiction that $({\hat{\mu}_{\mathrm{ML}}}^{(k)})_k$ is consistent but that $ \liminf_k \left|\tilde{\pi}^k_{s}\right| = L < +\infty $ for some $s \in (0, +\infty)$. Let $(k_j)_j$ be the corresponding subsequence. Divide all branches in ${\mathbb{T}}_{k_j}$ crossing $\tilde{\pi}^{k_j}_s$ into two branches meeting at $\tilde{\pi}^{k_j}_s$. Split edges in ${\mathbb{T}}_{k_j}$ above $\tilde{\pi}^{k_j}_s$ (closer to the root) repeatedly until the tree above $\tilde{\pi}^{k_j}_s$ forms a star. Let ${\mathbb{T}}'$ be the resulting tree, let $b'_1,\ldots,b'_{D}$ be the branches emanating from the root, where $D \leq L$ by assumption, and let $\tilde{\pi}'$ be the cutset at time $s$ from the leaves. For the unit flow ${{\boldsymbol\theta}}'$ corresponding to the MLE on ${\mathbb{T}}'$, by Lemma \[lemma:splitting\] and counting only those edges above $\tilde{\pi}'$ in ${\mathbb{T}}'$ in , we have $$\begin{aligned}
{\mathrm{Var}}_{{\mathbb{T}}_{k_j}}[{\hat{\mu}_{\mathrm{ML}}}^{(k_j)}]
&\geq& \gamma e^{-2\alpha T_{k_j}}
\left[
1 + \left(1 - e^{-2\alpha (T_{k_j} - s)}\right)e^{2\alpha (T_{k_j} - s)}
\sum_{i=1}^D (\theta'_{b'_i})^2
\right]\\
&\geq& \gamma 
\left[
e^{-2\alpha T_{k_j}} + \left(1 - e^{-2\alpha (T_{k_j} - s)}\right) \frac{e^{- 2\alpha s}}{L}
\right],\end{aligned}$$ where we used the fact that $\beta_1^2+\cdots+\beta_D^2$, subject to $\beta_1+\cdots+\beta_D = 1$, is minimized at $\beta_1=\cdots=\beta_D = 1/D$. Since $T_{k_j} \to +\infty$ under Assumption \[assumption:unbounded\], $$\limsup_k {\mathrm{Var}}_{{\mathbb{T}}_{k}}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}]
\geq \gamma \frac{e^{- 2\alpha s}}{L} > 0,$$ and we get a contradiction.

We note that, by Proposition \[prop:variance-branching\], the branching number provides a simple, sufficient condition for consistency.

\[thm:consistency-branching\] Let ${\mathcal{T}}= ({\mathbb{T}}_k)_k$ be a tree sequence satisfying Assumption \[assumption:unbounded\] with branching number ${\Lambda^\mathrm{b}}$. Then ${\Lambda^\mathrm{b}}> 0$ suffices for the consistency of the MLE of $\mu$.

Phase transition on the rate of convergence of the MLE {#sec:convergence-rate-mle}
------------------------------------------------------

Theorems \[thm:rate-supercritical\] and \[thm:subcritical\] show a phase transition for the $\sqrt{n_k}$-consistency of ${\hat{\mu}_{\mathrm{ML}}}$, which we prove now.

Assume ${\overline{\Lambda}^\mathrm{g}}> 2\alpha$. As remarked after Definition \[def:growth\], for all $\epsilon > 0$, eventually $$\label{eq:nk-bound}
\exp\left( 
({\underline{\Lambda}^\mathrm{g}}- \epsilon)
T_k\right)
\leq n_k \leq \exp\left( 
({\overline{\Lambda}^\mathrm{g}}+ \epsilon) 
T_k\right),$$ that is, $
n_k^{-2\alpha/({\underline{\Lambda}^\mathrm{g}}- \epsilon)}
\leq e^{-2\alpha T_k} 
\leq 
n_k^{-2\alpha/({\overline{\Lambda}^\mathrm{g}}+ \epsilon)}.
$ Moreover for all $\epsilon> 0$ there are subsequences $(k_j)_j$ and $(k'_j)_j$ such that $$\label{eq:subseq-supercritical}
n_{k_j} \geq \exp\left( 
({\overline{\Lambda}^\mathrm{g}}- \epsilon) 
T_{k_j}\right)
\mbox{ and }
n_{k'_j} \leq 
\exp\left( 
({\underline{\Lambda}^\mathrm{g}}+ \epsilon) 
T_{k'_j}\right).$$ By Proposition \[prop:variance-lower\], $$\label{eq:super-lower}
{\mathrm{Var}}_{{\mathbb{T}}_k}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}]
\geq \gamma\left[
e^{-2\alpha T_k}
+ \frac{1 - e^{-2\alpha T_k}}{n_k}
\right] \geq \gamma e^{-2\alpha T_k}.$$ Then  follows from  and . Hence $n_k {\mathrm{Var}}_{{\mathbb{T}}_k}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}] \to +\infty$ along a subsequence and $({\hat{\mu}_{\mathrm{ML}}}^{(k)})_k$ is not $\sqrt{n_k}$-consistent (using that ${\hat{\mu}_{\mathrm{ML}}}$ is unbiased and normally distributed).

Assume ${\Lambda^\mathrm{b}}> 2\alpha$. Let $2\alpha < \Lambda < {\Lambda^\mathrm{b}}$. By Proposition \[prop:variance-branching\] $$\label{eq:super-upper}
{\mathrm{Var}}_{{\mathbb{T}}_k}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}]
\leq \gamma\left[
1 + \frac{2\alpha}{{\mathcal{I}}_\Lambda (\Lambda - 2\alpha)}
\right]
e^{-2\alpha T_k}.$$ Note that ${\overline{\Lambda}^\mathrm{g}}\geq {\underline{\Lambda}^\mathrm{g}}\geq {\Lambda^\mathrm{b}}> 2\alpha$ and hence, by  and , $
{\mathrm{Var}}_{{\mathbb{T}}_k}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}] = \Theta\left(e^{-2\alpha T_k}\right).
$ Combining this with  gives the result in terms of $n_k$.

Assume instead that ${\Lambda^\mathrm{b}}< 2\alpha$. Let $\Lambda < {\Lambda^\mathrm{b}}$. By Proposition \[prop:variance-branching\] $${\mathrm{Var}}_{{\mathbb{T}}_k}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}]
\leq \gamma\left[
1 + \frac{2\alpha}{{\mathcal{I}}_\Lambda (2\alpha-\Lambda)}
\right]
e^{- \Lambda T_k}.$$ The rest of the argument is similar to the previous case.

Note that Example \[ex:twolevel\] considers a spherically symmetric tree. By , $${\mathrm{Var}}_{{\mathbb{T}}_k}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}] 
= \gamma e^{-2\alpha (\tau^{(k)}_0 + \tau^{(k)}_1)} \left[
1 + \sum_{h=0,1} (1 - e^{-2\alpha \tau^{(k)}_h})
 \prod_{h'=0}^{h}\frac{e^{2\alpha \tau^{(k)}_{h'}}}{D^{(k)}_{h'}} 
\right]$$ which then gives . Note that $$\frac{\log n_k}{ T_k}
= \frac{\Lambda_0\tau^{(k)}_0 + \Lambda_1 \tau^{(k)}_1}
{\tau^{(k)}_0 + \tau^{(k)}_1}
= \Lambda_0 \sigma + \Lambda_1 (1-\sigma)
= {\Lambda^\mathrm{g}}.$$ To compute the branching number, it suffices to consider cutsets with $m_0$ level-1 vertices and the $D^{(k)}_1(D^{(k)}_0 - m_0)$ tips below the rest of the level-1 vertices. Then $${\mathcal{J}}_k 
\equiv
\inf_{\pi \in \Pi^k}\ 
\sum_{x \in \pi} e^{- \Lambda \delta_k(\rho,x)}
=
\begin{cases}
D^{(k)}_0 e^{-\Lambda \tau^{(k)}_0},
& \text{if $D^{(k)}_1 > e^{\Lambda \tau^{(k)}_1}$}\\
n_k e^{- \Lambda T_k}, & \text{otherwise.}
\end{cases}$$ Hence if $\Lambda \geq \Lambda_1 > {\Lambda^\mathrm{g}}$ we are in the second case and $ n_k e^{- \Lambda T_k} = e^{- (\Lambda-{\Lambda^\mathrm{g}}) T_k} \to 0,$ as $k \to +\infty$. If $\Lambda < \Lambda_1$ we are in the first case and $ D^{(k)}_0 e^{-\Lambda \tau^{(k)}_0} = e^{- (\Lambda-\Lambda_0) \tau^{(k)}_0},$ so that ${\Lambda^\mathrm{b}}= \Lambda_0$.

One direction follows immediately from Proposition \[prop:variance-lower\] which implies $${\mathrm{Var}}_{{\mathbb{T}}_k}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}] \geq \gamma \frac{1 - e^{-2\alpha T_k}}{n_k}
= \Omega(n_k^{-1}).$$ We prove the other direction separately in each case. Assume first that $0 < {\Lambda^\mathrm{b}}= {\overline{\Lambda}^\mathrm{g}}< 2\alpha$. For $\epsilon > 0$ (small), choose $\Lambda$ such that $${\overline{\Lambda}^\mathrm{g}}- \epsilon < \Lambda < {\overline{\Lambda}^\mathrm{g}}= {\Lambda^\mathrm{b}}< 2\alpha.$$ By Proposition \[prop:variance-branching\], eventually $${\mathrm{Var}}_{{\mathbb{T}}_k}[{\hat{\mu}_{\mathrm{ML}}}^{(k)}]
\leq
\gamma\left[1 + \frac{2\alpha}{{\mathcal{I}}_\Lambda(2\alpha-\Lambda)}\right]
e^{- \Lambda T_k}
\leq 
\gamma\left[1 + \frac{2\alpha}{{\mathcal{I}}_\Lambda(2\alpha-\Lambda)}\right]
n_k^{-({\overline{\Lambda}^\mathrm{g}}-\epsilon)/({\overline{\Lambda}^\mathrm{g}}+ \epsilon)}.$$ Assume instead that ${\Lambda^\mathrm{ug}}< 2\alpha$. We show that $\overline{Y}$ (and hence the MLE by Proposition \[prop:key\]) achieves $\sqrt{n_k}$-consistency in this case. Let ${{\boldsymbol\theta}}$ be the corresponding flow on ${\mathbb{T}}_k$. By Corollary \[cor:integral-form\], letting ${\Lambda^\mathrm{ug}}< \Lambda < 2\alpha$, for $k$ large enough $$\begin{aligned}
{\mathrm{Var}}_{{\mathbb{T}}_k}[\overline{Y}] 
&=& 
\gamma e^{- 2 \alpha T_k} \left [ 1 + 
2\alpha \int_0^{T_k} e^{2\alpha s} 
\left(\sum_{x\in \pi^k_s} \left(\frac{n_k(x)}{n_k}\right)^2\right)
\mathrm{d}s \right ]\\
&\leq&
\gamma e^{- 2 \alpha T_k} \left [ 1 + 
2\alpha \int_0^{T_k} e^{2\alpha s} 
\left(\sum_{x\in \pi^k_s} \left(\frac{n_k(x)}{n_k}\right)
\frac{e^{\Lambda [(T_k-s)+M]}}{n_k}\right)
\mathrm{d}s \right ]\\
&\leq&
\gamma  \left [ e^{- 2 \alpha T_k} + 
e^{\Lambda M}
\frac{2\alpha}{n_k (2\alpha-\Lambda)} \,
(1 - e^{-(2\alpha -\Lambda) T_k})
\right ].\end{aligned}$$ The result follows from the fact that $e^{ \Lambda [T_k+M]} \geq n_k$.

Proofs for special cases {#section:proofs-special}
------------------------

By Theorems \[thm:criterion-consistency\], \[thm:rate-supercritical\], and \[thm:subcritical\], it suffices to prove that ${\Lambda^\mathrm{b}}= {\overline{\Lambda}^\mathrm{g}}= \lambda$ with probability $1$. A Galton-Watson (GW) branching process is a discrete-time non-negative integer-valued population process defined as follows: at each time step, each individual in the population has an independent number of offsprings, according to a distribution $F$, that form the population at the next time. In [@peres1999climb Chapter 3], it is shown that a GW tree where $F$ has mean $m$ has branching number and upper growth equal to $\log m$.

To compute the branching number of an infinite Yule tree ${\mathbb{T}}_0$, we use a comparison to a GW tree. Fix $\epsilon > 0$. Let $F$ be the distribution of the number of lineages in ${\mathbb{T}}_0$ at time $\epsilon$. By standard branching process results [@athreya2004branching Equation (4) on page 108], $m = e^{\lambda \epsilon}$. By the memoryless property of the exponential, the number of lineages $|\pi_{N\epsilon}|$ in the Yule tree at time $N\epsilon$ is identically distributed to the population size $Z_{N}$ of a GW tree with offspring distribution $F$ at time $N$. Then $$\frac{\log |\pi_s|}{s} 
    \leq \frac{\log Z_{\lceil s/\epsilon \rceil}}{s}
    = \frac{\lceil s/\epsilon \rceil}{ s} \cdot 
    \frac{\log Z_{\lceil s/\epsilon \rceil}}{\lceil s/\epsilon \rceil},$$ which implies that $
    {\overline{\Lambda}^\mathrm{g}}\leq \frac{1}{\epsilon} \cdot
    \log e^{\lambda \epsilon} = \lambda.
    $

Similarly, let $\pi$ be a cutset in ${\mathbb{T}}_0$ and let $\pi_\epsilon$ be the cutset obtained by rounding up the points in $\pi$ to the next $\epsilon$-multiple closer to the root (removing duplicates). Let $\delta_{\mathrm{GW}}(v)$ be the distance from the root to vertex $v$ in the GW tree. Then $$\sum_{x \in \pi} 
    e^{- \Lambda \delta_0(\rho,x)}
    \geq \sum_{y \in \pi_\epsilon} 
    e^{- \Lambda (\delta_{\mathrm{GW}}(y)+1)\epsilon}
    = e^{- \Lambda \epsilon} \sum_{y \in \pi_\epsilon} 
    e^{- (\epsilon \Lambda) \delta_{\mathrm{GW}}(y)}
    > 0$$ whenever $\epsilon \Lambda < \log e^{\lambda \epsilon}$, so that $
    {\Lambda^\mathrm{b}}\geq \lambda.
    $

Let $\tau_i = t_{i} - t_{i-1}$ be the amount of time during which ${\mathbb{T}}_0$ has $i$ lineages (with $t_0 = 0$). Then $(\tau_i)_i$ are independent exponential random variables with parameters $(1/(i\lambda))_i$. Let $T^j_i =\sum_{r=i+1}^j \tau_r$. Note that $$    {{{\mathrm{I}\!\mathrm{E}}}\!\left[T^j_i\right]} 
    = \sum_{r=i+1}^j {{{{\mathrm{I}\!\mathrm{E}}}\!\left[\tau_r\right]}}
    \equiv \lambda^{-1} \sum_{r=i+1}^j \frac{1}{r}
    \in \left(\lambda^{-1}\log\Big(\frac{j}{i+1}\Big), \lambda^{-1}\log\Big(\frac{j}{i}\Big)\right).$$ Similarly, $$\label{eq:yule-alpha-var}
    {\mathrm{Var}}[T^j_i] 
    = \sum_{r=i+1}^j {\mathrm{Var}}[\tau_r] 
    = \lambda^{-2}
    \sum_{r=i+1}^j \frac{1}{r^2} \leq \frac{1}{\lambda^2 i}.$$ By Chebyshev’s inequality, for all $0<\sigma<1$, $${{{\mathrm{I}\!\mathrm{P}}}\!\left[
        |T_{\lfloor \sigma k \rfloor}^k
        - {{{\mathrm{I}\!\mathrm{E}}}\!\left[T_{\lfloor \sigma k \rfloor}^k\right]}|
        \geq \epsilon
    \right]} = O(k^{-1}),$$ where we used . Let $0 < \sigma_2' < \sigma_2
    < \sigma_1' < \sigma_1 < 1$. From the previous equation, we get for $\iota=1,2$ $${{{\mathrm{I}\!\mathrm{P}}}\!\left[
        T_{\lfloor \sigma_\iota k \rfloor}^k
        \leq \lambda^{-1}\log\left(\frac{k}{\lfloor \sigma_\iota k\rfloor + 1}\right) - \epsilon
    \right]} = O(k^{-1}),$$ and similarly for the other direction. Take $$a_\iota
    = \lambda^{-1}
    \log\left(\frac{1}{\sigma_\iota}\right),
    \quad
    a_\iota'
    = \lambda^{-1}
    \log\left(\frac{1}{\sigma_\iota'}\right)
    \quad\mbox{and }
    \epsilon 
    < a_1     \land
    \frac{1}{2}
    \left[ a_2      -      a_1'     \right].$$ Then Assumption \[assump:alpha\] is satisfied asymptotically with $
    c_\iota=a_\iota-\epsilon$, $c_\iota' = a_\iota'+\epsilon$ and $\beta 
    = [\sigma_1 - \sigma_1']
    \land
    [\sigma_2 - \sigma_2'],
    $ because then $${{{\mathrm{I}\!\mathrm{P}}}\!\left[c_1 < T_{\lfloor \sigma_1 k \rfloor}^k < T_{\lfloor \sigma_1' k \rfloor}^k< 
        c'_1 <c_2 <
        T_{\lfloor \sigma_2 k \rfloor}^k < T_{\lfloor \sigma_2' k \rfloor}^k<
        c'_2
    \right]}
    \geq 1-O(k^{-1}).$$

Sensitivity to estimate of $\alpha$ {#sec:sensitivity-to-alpha}
-----------------------------------

So far in this section, we considered the MLE of $\mu$ [*given $\alpha$*]{}. Here we look at the sensitivity of the MLE to estimation errors on $\alpha$. Theorem \[thm:main-alpha-rate\] shows that there exists a $\sqrt{n_k}$-consistent estimator of $\alpha$ under Assumption \[assump:alpha\], which is unrelated to the growth or height of the species tree. Moreover the estimator of $\alpha$ we derive does not require the knowledge of $\mu$.

Hence suppose that we have a $\sqrt{n_k}$-consistent estimator $\hat\alpha_k$ of $\alpha$. Let ${\widehat{\mathrm{Var}}}_{{\mathbb{T}}_k}$ denote the variance under the parameter $\alpha = \hat\alpha_k$ (with $\mu$ and $\gamma$ unchanged) and let $\hat{{\boldsymbol\theta}}_k$ be the corresponding weights of the MLE of $\mu$, that is, the choice of weights assuming that $\alpha = \hat\alpha_k$ and minimizing ${\widehat{\mathrm{Var}}}_{{\mathbb{T}}_k}[Y_{{{\boldsymbol\theta}}}]$.

For all $k$ and under the true $\alpha$, $Y_{\hat{{\boldsymbol\theta}}_k}$ is an unbiased estimator of $\mu$. Moreover, because $\hat\alpha_k = \alpha + o(1)$ and so on, the bounds in Theorems \[thm:rate-supercritical\] and \[thm:subcritical\] apply to ${\widehat{\mathrm{Var}}}_{{\mathbb{T}}_k}[Y_{\hat{{\boldsymbol\theta}}_k}]$ as well (for $k$ large enough). The quantity of interest is ${\mathrm{Var}}_{{\mathbb{T}}_k}[Y_{\hat{{\boldsymbol\theta}}_k}]$. By , $$\begin{aligned}
&&{\mathrm{Var}}_{{\mathbb{T}}_k}[Y_{\hat{{\boldsymbol\theta}}_k}] \\
&&\quad =   
\gamma e^{- 2 \alpha T_k} 
+ \gamma \sum_{b \in E_k}(1 - e^{-2 \alpha |b|}) 
e^{2 \alpha (\delta_k(\rho,b)-T_k)}
(\hat\theta_k)^2_b \\
&&\quad = (1+O(T_k n_k^{-1/2}))  
\left [\gamma e^{- 2 \hat\alpha_k T_k} 
+ \gamma \sum_{b \in E_k}(1 - e^{-2 \hat\alpha_k |b|}) 
e^{2 \hat\alpha_k (\delta_k(\rho,b)-T_k)}
(\hat\theta_k)^2_b \right ]\\
&&\quad = (1+O(T_k n_k^{-1/2})) {\widehat{\mathrm{Var}}}_{{\mathbb{T}}_k}[Y_{\hat{{\boldsymbol\theta}}_k}],\end{aligned}$$ provided $T_k n_k^{-1/2} = o(1)$. Hence, for instance if ${\underline{\Lambda}^\mathrm{g}}> 0$, $T_k = O(\log n_k)$ and we get that ${\mathrm{Var}}_{{\mathbb{T}}_k}[Y_{\hat{{\boldsymbol\theta}}_k}]$ satisfies the bounds in Theorems \[thm:rate-supercritical\] and \[thm:subcritical\].

Convergence rate of a new estimator for $\alpha$ and $\gamma$ {#sec:alpha_gamma}
=============================================================

In this section, we provide a novel estimator for $(\alpha, \gamma)$. Under natural assumptions on the species tree, we show that this estimator is $\sqrt{n_k}$-consistent. Moreover this estimator does not require the knowledge of $\mu$. Interestingly, in contrast to what we showed for $\mu$, the conditions for $\sqrt{n_k}$-consistency in this case do not involve the growth, or even the height, of the species tree. This is in line with the results in [@HoAne13], who found that $\mu$ requires an unbounded tree height to be microergodic, whereas $\alpha$ and $\gamma$ do not.

Note, however, that the MLE of $\alpha$ and $\gamma$ are not simple linear estimators, which makes them harder to study here. In particular, unlike in the case of $\mu$, we do not provide lower bounds on their rate of convergence.

Contrast-based estimator {#sec:contrasts}
------------------------

We first describe the estimator. The proof of its convergence rate is in Section \[sec:alpha-rate\].

#### Contrasts

Our estimator relies on an appropriately chosen set of contrasts, that is, differences between pairs of leaf states (see e.g. [@felsenstein2004inferring]). More specifically, we choose contrasts associated with internal nodes, as follows. Let ${\mathbb{T}}$ be an ultrametric species tree with leaves ${\mathscr{L}}$ and internal vertices ${\mathscr{I}}$. For two leaves $\ell$ and $\ell'$, we let $\ell\land \ell'$ be their most recent common ancestor. Assume that all internal vertices of ${\mathbb{T}}$ have out-degree at least $2$. Let $i \in {\mathscr{I}}$ be an internal vertex of ${\mathbb{T}}$, and let $\ell^i_1\neq\ell^i_2$ be two leaves such that $\ell^i_1 \land \ell^i_2 = i$. Let $P_i$ be the path connecting $\ell^i_1$ and $\ell^i_2$. We define the corresponding contrast ${\mathcal{C}}_i = Y_{\ell^i_1} - Y_{\ell^i_2}$. Let $T(i)$ be the height of $i$ from the leaves. We say that $T(i)$ is the height of ${\mathcal{C}}_i$.

Let $i_1,\ldots,i_m$ be a collection of internal nodes of ${\mathbb{T}}$. Let ${\mathcal{C}}_{i_1}, \ldots, {\mathcal{C}}_{i_m}$ be an arbitrary set of associated contrasts. Assume that the corresponding paths $P_{i_1},\ldots,P_{i_m}$ are pairwise non-intersecting, that is, none of the pairs of paths share a vertex. Then ${\mathcal{C}}_{i_1}, \ldots, {\mathcal{C}}_{i_m}$ are mutually independent, multivariate normal with ${\mathcal{C}}_i \sim {\cal N}(0,2\gamma(1-e^{-2\alpha T(i)}))$.\[lem\_addmic06\]

Indeed, expanding the covariance, we get for $j\neq j'$ $$\begin{aligned}
\gamma^{-1}{\mathrm{Cov}}[{\mathcal{C}}_j, {\mathcal{C}}_{j'}]
&=& e^{-\alpha d_{\ell^{j}_1\ell^{j'}_1}}
- e^{-\alpha d_{\ell^{j}_1\ell^{j'}_2}}
- e^{-\alpha d_{\ell^{j}_2\ell^{j'}_1}}
+ e^{-\alpha d_{\ell^{j}_2\ell^{j'}_2}} = 0,\end{aligned}$$ since, by assumption, $\ell^{j}_\iota \land \ell^{j'}_{\iota'}$ is the same vertex for all $\iota, \iota' = 1,2$.

The following lemma will be useful in identifying an appropriate collection of contrasts.

Let $\mathbb{T}$ be an ultrametric tree and let ${\mathscr{I}}_{(a,b)}$ be the set of internal nodes of $\mathbb{T}$ whose height from the leaves lies in $(a,b)$. For every $a < b$, we can select a set of independent contrasts ${\mathscr{C}}$, associated with internal nodes in ${\mathscr{I}}(a,b)$, such that $$|{\mathscr{C}}| \geq n(a,b)/2,$$ where $n(a,b) = |{\mathscr{I}}(a,b)| $. In particular, the heights of the contrasts in ${\mathscr{C}}$ lie in $(a,b)$ and their corresponding paths are pairwise non-intersecting. \[lem:choose\]

Start with the lowest vertex $i$ in ${\mathscr{I}}_{(a,b)}$ and choose a pair of vertices $\ell^i_1$ and $\ell^i_2$ such that $\ell^i_1 \land \ell^i_2 = i$. Remove $i$ and its descendants as well as the edge immediately above $i$ (and fuse consecutive edges separated by degree-$2$ vertices). As a result, the number of internal vertices in $(a,b)$ decreases by at most $2$. Repeat until no vertex is left in ${\mathscr{I}}_{(a,b)}$.

#### The estimator

For a sequence of trees ${\mathcal{T}}= ({\mathbb{T}}_k)_k$, let ${\mathscr{L}}_k$ be the leaf set of ${\mathbb{T}}_k$; ${\mathscr{I}}_k$, the set of its internal vertices; $n_k = |{\mathscr{L}}_k|$ and $n_k(a,b) = |{\mathscr{I}}_k(a,b)|$; and $T_k(i)$, the height of $i$, for each $i\in{\mathscr{I}}_k$. The idea behind our estimator is to set up a system of equations that characterize $\alpha$ and $\gamma$ uniquely. Our construction relies on the following condition. We illustrate this condition on two special cases below.

We set up our equations as follows. Let $m_k = \lfloor \beta n_k/2 \rfloor$. Under Assumption \[assump:alpha\], by Lemma \[lem:choose\], for each $k$ we can choose [*two*]{} collections of independent contrasts $({\mathcal{C}}^k_{i_r})_{r=1}^{m_k}$ and $({\mathcal{C}}^k_{j_r})_{r=1}^{m_k}$ with corresponding heights $T_k(i_r) \in (c_1,c_1')$ and $T_k(j_r) \in (c_2,c_2')$ for every $r = 1, 2, \ldots, m_k$. (Note that the two collections are [*not*]{} independent.) For $r = 1,\ldots,m$, let $$\begin{aligned}
\hat a_{k} 
= \frac{1}{m_k} \sum_{r=1}^{m_k}{({\mathcal{C}}^{(k)}_{i_r})^2},
\qquad
\hat b_{k} 
= \frac{1}{m_k} \sum_{r=1}^{m_k}{({\mathcal{C}}^{(k)}_{j_r})^2},\end{aligned}$$ and note that $$\begin{aligned}
a_{k} &\equiv& {{{\mathrm{I}\!\mathrm{E}}}\!\left[\hat a_{k}\right]} = 2 \gamma \left (1 - \frac{1}{m_k} \sum_{r=1}^{m_k}{e^{- 2 \alpha T_k(i_r)}} \right )
\equiv 2\gamma h^1_{k}(\alpha),\\
b_{k} &\equiv& {{{\mathrm{I}\!\mathrm{E}}}\!\left[\hat b_{k}\right]} = 2 \gamma \left (1 - \frac{1}{m_k} \sum_{r=1}^{m_k}{e^{- 2 \alpha T_k(j_r)}} \right )
\equiv 2\gamma h^2_{k}(\alpha).\end{aligned}$$ Notice that, under Assumption \[assump:alpha\], $a_{k} \in [2\gamma(1-e^{-2 \alpha c_2}),2\gamma(1-e^{-2 \alpha c_1})] \equiv [\underline{a}_\alpha,\bar{a}_\alpha]$ and $b_{k} \in [2\gamma(1-e^{-2 \alpha c_4}),2\gamma(1-e^{-2 \alpha c_3})] \equiv [\underline{b}_\alpha,\bar{b}_\alpha]$. As shown below, $$H_k(\alpha) = \frac{a_k}{b_k} = \frac{h^1_{k}(\alpha)}
{h^2_{k}(\alpha)}$$ is invertible in $\alpha$ on $(0,+\infty)$. Hence a natural estimator of $(\alpha,\gamma)$ is obtained by setting $$\hat\alpha_k = H_k^{-1}\left(\frac{\hat a_k}{\hat b_k}\right) 
\quad\mbox{and}\quad
\hat\gamma_k = \frac{\hat a_k}{2 h_k^1(\hat\alpha_k)}.$$ We will show in the proof of invertibility below that $H_k$ is actually strictly increasing, and therefore relatively straightforward to invert numerically. It remains to prove invertibility.

\[lemma:invertibility\] Under Assumption \[assump:alpha\], $H_{k}(\alpha)$ is strictly positive, differentiable, and invertible on $(0,+\infty)$.

We have that $$\begin{aligned}
\frac{\partial \log H_{k}(\alpha)}{\partial \alpha} 
&=& \frac{\sum_{r=1}^{m_k} {2T_k(i_r)}e^{-2 \alpha T_k(i_r)}}{\sum_{r=1}^{m_k}{(1 - e^{-2 \alpha T_k(i_r)}})} - \frac{\sum_{r=1}^{m_k}{2T_k(j_r) e^{-2 \alpha T_k(j_r)}}}{\sum_{r=1}^{m_k}{(1 - e^{-2 \alpha T_k(j_r)})}}\nonumber\\
&=& \frac{\sum\sum_{r, r'=1}^{m_k} {2T_k(i_r)}e^{-2 \alpha T_k(i_r)}
{(1 - e^{-2 \alpha T_k(j_{r'})})}}{\sum_{r=1}^{m_k}{(1 - e^{-2 \alpha T_k(i_r)}}) \sum_{r=1}^{m_k}{(1 - e^{-2 \alpha T_k(j_r)})}}\nonumber\\
&& \qquad - \frac{\sum\sum_{r,r'=1}^{m_k}{2 T_k(j_{r'}) e^{-2 \alpha T_k(j_{r'})}}(1 - e^{-2 \alpha T_k(i_{r})})}{\sum_{r=1}^{m_k}{(1 - e^{-2 \alpha T_k(i_r)}})\sum_{r=1}^{m_k}{(1 - e^{-2 \alpha T_k(j_r)})}}
\label{eq:invertibility-1}\end{aligned}$$ Note that the function $\frac{xe^{-x}}{1 - e^{-x}}$ is strictly decreasing on $(0, \infty)$ because its derivative is $\frac{e^{-x}(1 - x - e^{-x})}{(1 - e^{-x})^2} < 0$ on $(0,+\infty)$. Therefore $$\frac{2T_k(i_r)e^{-2 \alpha T_k(i_r)}}{1 - e^{-2 \alpha T_k(i_r)}} \geq 
\frac{2 c_1' e^{-2 \alpha c_1'}}{1 - e^{-2 \alpha c_1'}} >
\frac{2 c_2 e^{-2 \alpha c_2 }}{1 - e^{-2 \alpha c_2}} \geq
\frac{2T_k(j_{r'}) e^{-2 \alpha T_k(j_{r'})}}{1 - e^{-2 \alpha T_k(j_{r'})}},$$ that is, $$\begin{aligned}
&&2T_k(i_r)e^{-2 \alpha T_k(i_r)} (1 - e^{-2 \alpha T_k(j_{r'})})\nonumber\\
&& \qquad - 2T_k(j_{r'}) e^{-2 \alpha T_k(j_{r'})} (1 - e^{-2 \alpha T_k(i_r)})
> 0,\label{eq:logh-difference}\end{aligned}$$ for every $r,r'$, so that each $(r,r')$-term in  is strictly positive. Hence, we can deduce that $\partial \log H_k(\alpha)/\partial \alpha > 0$, that is, $\log H_k$ (and hence $H_k$ itself) is strictly increasing on $(0,+\infty)$ and continuous, and therefore invertible.

Note that we cannot use the law of large numbers to derive consistency (despite the independence of the contrasts) because $a_k/b_k$ is a bounded, but not necessarily convergent, sequence and $H_k^{-1}$ is continuous, but depends on $k$. Instead we argue directly about $\sqrt{n_k}$-consistency below.

Proof of Theorem \[thm:main-alpha-rate\] {#sec:alpha-rate}
----------------------------------------

Note that ${{{\mathrm{I}\!\mathrm{E}}}\!\left[\hat a_k\right]} = a_k$ and $${\mathrm{Var}}[\hat a_k]
= \frac{8\gamma^2}{m_k^2}
\sum_{r=1}^{m_k} 
(1-e^{-2\alpha T_k(i_r)})^2
\leq \frac{8\gamma^2}{m_k}
(1-e^{-2\alpha c_1})^2 
= O(m_k^{-1})
= O(n_k^{-1}),$$ where we used that $([2\gamma(1-e^{-2\alpha T_k(i_r)})]^{-1/2} {\mathcal{C}}^k_{i_r})^2 $ is $\chi^2_1$-distributed and, therefore, has variance $2$. Hence $|\hat a_k - a_k| = O_p(n_k^{-1/2})$ by Chebyshev’s inequality. Similarly, $|\hat b_k - b_k| = O_p(n_k^{-1/2})$. Our claim that $|\hat\alpha_k - \alpha_k| = O_p(n_k^{-1/2})$ then follows from the following straightforward lemma.

If $0 < z_* \leq z \leq z^* < \infty$, $|z' - z| \leq \epsilon$ and $\epsilon < z_*/2$, then there is a constant $\Delta(z_*,z^*)$ depending on $c_1,c_1',c_2,c_2'$ such that for all $k$ $$\sup_{t \in [0,1]}{|(H_k^{-1})' (t z' + (1-t)z )|} \leq \Delta(z_*,z^*).$$ \[lem:monotone\]

We use the proof of Lemma \[lemma:invertibility\]. Let $\zeta_\alpha = \zeta_\alpha(c_1,c_1',c_2,c_2') > 0$ be the smallest possible difference in  for a fixed $\alpha$. Let $\alpha_*, \alpha^*$ be defined as $$\frac{1}{2}z_* = \frac{\bar{a}_{\alpha_*}}{\underline b_{\alpha_*}},
\qquad \frac{3}{2}z^* = \frac{\underline{a}_{\alpha^*}}{\bar b_{\alpha^*}}.$$ Then $[\alpha_*,\alpha^*]
\supseteq H_k^{-1}\left(\left[\frac{1}{2}z_*,\frac{3}{2}z^*\right]\right)$ for all $k$. Note that $$\begin{aligned}
\sup_{t \in [0,1]}{|(H_k^{-1})' (t z' + (1-t)z )|}
&\leq& \sup_{z\in\left[\frac{1}{2}z_*,\frac{3}{2}z^*\right]}\left|(H_k^{-1})'(z)\right|\\
&=& \sup_{z\in\left[\frac{1}{2}z_*,\frac{3}{2}z^*\right]}\left(\frac{\partial H_k}{\partial \alpha} (H_k^{-1}(z))\right)^{-1}\\
&=& \sup_{z\in\left[\frac{1}{2}z_*,\frac{3}{2}z^*\right]}\left(\left[H_k\frac{\partial \log H_k}{\partial \alpha}\right] (H_k^{-1}(z))\right)^{-1}\\
&\leq& \sup_{\alpha\in[\alpha_*,\alpha^*]}\frac{\underline{b}_\alpha}{\bar a_\alpha}
\cdot \frac{(1-e^{-2\alpha c_1'})(1 - e^{-2\alpha c_2'})}{\zeta_\alpha}\\
&\equiv& \Delta(z_*,z^*).\end{aligned}$$

We finish the proof of Theorem \[thm:main-alpha-rate\]. We use the following observation: for $0 < x_* \leq x \leq x^* < \infty $ and $0 < y_* \leq y \leq y^* < \infty$ such that $|x - x'| \leq \epsilon$ and $|y - y'| \leq \epsilon$ with $\epsilon < y_*/2$, we have $$\left | \frac{x'}{y'} - \frac{x}{y} \right |
= \left | \frac{y(x' - x) + x(y - y')}{y y'} \right | \leq \frac{y^*|x' - x|}{y_*(y_*/2)} + \frac{x^*|y - y'|}{y_*(y_*/2)}
 < \frac{4(x^* + y^*)}{y_*^2}\epsilon.$$ Fix $\delta > 0$ (small) and pick $M_\delta$ such that ${{{\mathrm{I}\!\mathrm{P}}}\!\left[|\hat a_k - a_k|\geq M_\delta n_k^{-1/2}\right]} <\delta/2$ and similarly for $\hat b_k$. Then, by Assumption \[assumption:unbounded\], for $k$ large enough [$$\begin{aligned}
&&{{{\mathrm{I}\!\mathrm{P}}}\!\left[\left | \frac{\hat a_k}{\hat b_k} - \frac{a_k}{b_k} \right | 
\geq \frac{4(\bar a_\alpha + \bar b_\alpha)}{\underline{b}_\alpha^2} M_\delta n_k^{-1/2}\right]}\\ 
&&\qquad \leq {{{\mathrm{I}\!\mathrm{P}}}\!\left[ \left |\frac{\hat a_k}{\hat b_k} - \frac{a_k}{b_k} \right | \geq \frac{4(\bar a_\alpha + \bar b_\alpha)}{\underline{b}_\alpha^2} M_\delta n_k^{-1/2}, |\hat a_k - a_k| \leq M_\delta n_k^{-1/2}, |\hat b_k - b_k| \leq M_\delta n_k^{-1/2}\right]}\\
&& \qquad\qquad  + {{{\mathrm{I}\!\mathrm{P}}}\!\left[|\hat a_k - a_k| \geq M_\delta n_k^{-1/2}\right]} 
+ {{{\mathrm{I}\!\mathrm{P}}}\!\left[|\hat b_k - b_k| \geq M_\delta n_k^{-1/2}\right]}\\
&&\qquad \leq 0 + \frac{\delta}{2} + \frac{\delta}{2} = \delta,\end{aligned}$$]{} so that $
\left | \frac{\hat a_k}{\hat b_k} - \frac{a_k}{b_k} \right |
= O_p(n_k^{-1/2}).
$

Secondly, using Rolle’s theorem, we have $$|\hat \alpha_k - \alpha| 
\leq \sup_{t \in [0,1]}{\left |(H_k^{-1})' \left (t\frac{\hat a_k}{\hat b_k} + (1-t)\frac{a_k}{b_k} \right ) \right|. \left |\frac{\hat a_k}{\hat b_k} - \frac{a_k}{b_k} \right |}.$$ Let $M_\delta$ be such that $${{{\mathrm{I}\!\mathrm{P}}}\!\left[\left|\frac{\hat a_k}{\hat b_k} - \frac{a_k}{b_k}\right|\geq M_\delta n_k^{-1/2}\right]} <\delta.$$ Fix $\epsilon' > 0$ and let $$z_* = \frac{\underline{a}_{\alpha - \epsilon'}}{\bar{b}_{\alpha - \epsilon'}}, \qquad 
z^* = \frac{\bar{a}_{\alpha + \epsilon'}}{\underline{b}_{\alpha + \epsilon'}}.$$ Then, by Lemma \[lem:monotone\], letting $$\mathcal{H}_k
= \left\{
\sup_{t \in [0,1]}{\left |(H_k^{-1})' \left (t\frac{\hat a_k}{\hat b_k} + (1-t)\frac{a_k}{b_k} \right ) \right|. 
\left |\frac{\hat a_k}{\hat b_m} - \frac{a_m}{b_m} \right |} \geq \Delta^{-1}(z_*,z^*) M_\delta n_k^{-1/2}
\right\},$$ we have for $k$ large enough $$\begin{aligned}
&&{{{\mathrm{I}\!\mathrm{P}}}\!\left[|\hat \alpha_k - \alpha| \geq \Delta^{-1}(z_*,z^*) M_\delta n_k^{-1/2}\right]}\\
&&\quad\leq {{{\mathrm{I}\!\mathrm{P}}}\!\left[\mathcal{H}_k\right]}\\
&&\quad\leq {{{\mathrm{I}\!\mathrm{P}}}\!\left[\mathcal{H}_k, \left |\frac{\hat a_k}{\hat b_k} - \frac{a_k}{b_k} \right | < M_\delta n_k^{-1/2}\right]}
+ {{{\mathrm{I}\!\mathrm{P}}}\!\left[\left |\frac{\hat a_k}{\hat b_k} - \frac{a_k}{b_k} \right | \geq M_\delta n_k^{-1/2}\right]}\\
&&\quad\leq 0 + \delta 
= \delta.\end{aligned}$$ That implies $ |\hat\alpha_k - \alpha|= O_p(n_k^{-1/2}). $ The argument for $\hat\gamma_k$ is similar.

[10]{} \[1\][[\#1]{}]{} urlstyle \[1\][DOI \#1]{}

Adamczak, R., Mi[ł]{}o[ś]{}, P.: [CLT]{} for [O]{}rnstein-[U]{}hlenbeck branching particle system. Electronic Journal of Probability **20**(42), 1–35 (2015)

Adamczak, R., Mi[ł]{}o[ś]{}, P.: [U]{}-[S]{}tatistics of [O]{}rnstein-[U]{}hlenbeck branching particle system. Journal of Theoretical Probability **27**(4), 1071–1111 (2014)

Anderson, T.W.: An introduction to multivariate statistical analysis, 2nd edn. Wiley, Chichester (1984)

Athreya, K., Ney, P.: Branching Processes. Dover Books on Mathematics Series. Dover Publications (2004)

Bartoszek, K., Pienaar, J., Mostad, P., Andersson, S., Hansen, T.F.: A phylogenetic comparative method for studying multivariate adaptation. Journal of Theoretical Biology **314**, 204–215 (2012)

Bartoszek, K., Sagitov, S.: Phylogenetic confidence intervals for the optimal trait value. Journal of Applied Probability **52**(4), 1115–1132 (2015).

Bininda-Emonds, O., Cardillo, M., Jones, K.E., MacPhee, R.D.E., Beck, R.M.D., Grenyer, R., Price, S.A., Vos, R.A., Gittleman, J.L., Purvis, A.: The delayed rise of present-day mammals. Nature **446**(7135), 507–512 (2007)

Brawand, D., Soumillon, M., Necsulea, A., Julien, P., Csardi, G., Harrigan, P., Weier, M., Liechti, A., Aximu-Petri, A., Kircher, M., Albert, F.W., Zeller, U., Khaitovich, P., Grutzner, F., Bergmann, S., Nielsen, R., Pääbo, S., Kaessmann, H.: The evolution of gene expression levels in mammalian organs. Nature **478**(7369), 343–348 (2011)

Butler, M.A., King, A.A.: Phylogenetic comparative analysis: a modeling approach for adaptive evolution. The American Naturalist **164**(6), 683–695 (2004)

Cooper, N., Purvis, A.: Body size evolution in mammals: Complexity in tempo and mode. The American Naturalist **175**(6), 727–738 (2010)

Crawford, F.W., Suchard, M.A.: Diversity, disparity, and evolutionary rate estimation for unresolved [Y]{}ule trees. Systematic Biology **62**(3), 439–455 (2013)

Evans, W.S., Kenyon, C., Peres, Y., Schulman, L.J.: Broadcasting on trees and the [I]{}sing model. Ann. Appl. Probab. **10**(2), 410–433 (2000)

Felsenstein, J.: Phylogenies and the comparative method. American Naturalist **125**(1), 1–15 (1985)

Felsenstein, J.: Inferring Phylogenies. Sinauer Associates (2004)

Hansen, T.F.: Stabilizing selection and the comparative analysis of adaptation. Evolution **51**(5), 1341–1351 (1997)

Harmon, L., Weir, J., Brock, C., Glor, R., Challenger, W.: G[EIGER]{}: investigating evolutionary radiations. Bioinformatics **24**, 129–131 (2008)

Harmon, L.J., Losos, J.B., Jonathan Davies, T., Gillespie, R.G., Gittleman, J.L., Bryan Jennings, W., Kozak, K.H., McPeek, M.A., Moreno-Roark, F., Near, T.J., Purvis, A., Ricklefs, R.E., Schluter, D., Schulte II, J.A., Seehausen, O., Sidlauskas, B.L., Torres-Carvajal, O., Weir, J.T., Mooers, A.Ø.: Early bursts of body size and shape evolution are rare in comparative data. Evolution **64**(8), 2385–2396 (2010)

Ho, L.S.T., Ané, C.: Asymptotic theory with hierarchical autocorrelation: [Ornstein-Uhlenbeck]{} tree models. Annals of Statistics **41**, 957–981 (2013)

Ho, L.S.T., An[é]{}, C.: Intrinsic inference difficulties for trait evolution with [Ornstein-Uhlenbeck]{} models. Methods in Ecology and Evolution **5**(11), 1133–1146 (2014)

Ho, L.S.T., Ané, C.: A linear-time algorithm for [G]{}aussian and non-[G]{}aussian trait evolution models. Systematic Biology **63**(3), 397–408 (2014)

Jetz, W., Thomas, G., Joy, J., Hartmann, K., Mooers, A.: The global diversity of birds in space and time. Nature **491**(7424), 444–448 (2012)

Lawler, E.: Combinatorial Optimization: Networks and Matroids. Holt, Rinehart and Winston (1976)

Mossel, E., Roch, S., Sly, A.: Robust estimation of latent tree graphical models: Inferring hidden states with inexact parameters. IEEE transactions on information theory **59**(7), 4357–4373 (2013)

Mossel, E., Steel, M.: Majority rule has transition ratio 4 on yule trees under a 2-state symmetric model. Journal of Theoretical Biology **360**(7), 315–318 (2014).

Paradis, E., Claude, J., Strimmer, K.: A[PE]{}: analyses of phylogenetics and evolution in [R]{} language. Bioinformatics **20**, 289–290 (2004)

Peres, Y.: Probability on trees: An introductory climb. In: P. Bernard (ed.) Lectures on Probability Theory and Statistics, *Lecture Notes in Mathematics*, vol. 1717, 193–280. Springer Berlin Heidelberg (1999)

Rohlfs, R.V., Harrigan, P., Nielsen, R.: Modeling gene expression evolution with an extended [O]{}rnstein-[U]{}hlenbeck process accounting for within-species variation. Molecular Biology and Evolution **31**(1), 201–211 (2014)

Semple, C., Steel, A.: Phylogenetics. Oxford lecture series in mathematics and its applications. Oxford University Press (2003)

Shao, J.: Mathematical Statistics. Springer (2003)

Venditti, C., Meade, A., Pagel, M.: Multiple routes to mammalian diversity. Nature **479**(7373), 393–396 (2011)

Yule, G.U.: A mathematical theory of evolution, based on the conclusions of [Dr. JC Willis, FRS]{}. Philosophical Transactions of the Royal Society of London. Series B **213**, 21–87 (1925)

[^1]: Departments of Statistics and of Botany, University of Wisconsin-Madison. Work supported by NSF grants DMS-1106483.

[^2]: Departments of Statistics, University of Wisconsin-Madison.

[^3]: Departments of Mathematics and Statistics (by courtesy), University of Wisconsin-Madison. Work supported by NSF grants DMS-1007144 and DMS-1149312 (CAREER), and an Alfred P. Sloan Research Fellowship.
