Phylogenetic confidence intervals for the optimal trait value

TL;DR

This paper introduces a new method to calculate confidence intervals for optimal trait values in evolutionary models involving related species.

Contribution

The novelty lies in deriving confidence interval formulae for the optimal trait value using limit theorems from a stochastic evolutionary model.

Findings

Three limit theorems for the sample mean are established based on adaptation rate domains.

For fast adaptation, the normalized sample mean is approximately normally distributed for large n.

New confidence interval formulae for the optimal trait value are developed using these theorems.

Abstract

We consider a stochastic evolutionary model for a phenotype developing amongst n related species with unknown phylogeny. The unknown tree is modelled by a Yule process conditioned on n contemporary nodes. The trait value is assumed to evolve along lineages as an Ornstein-Uhlenbeck process. As a result, the trait values of the n species form a sample with dependent observations. We establish three limit theorems for the sample mean corresponding to three domains for the adaptation rate. In the case of fast adaptation, we show that for large n the normalized sample mean is approximately normally distributed. Using these limit theorems, we develop novel confidence interval formulae for the optimal trait value.

Figures

On the left: a branching Ornstein-Uhlenbeck process simulated on a realization of the Yule n-tree with n = 5 tips using the TreeSim [38, 39] and mvS-LOUCH [7] R [32] packages. Parameters used are α = 1, σ = 1, X0 −θ = 2, after the tree height Un was scaled to 1. On the right: the species tree disregarding the trait values supplied with the notation for the inter-speciation times. For the pair of tips (2,3) the time τ(n) to their most recent common ancestor is marked on the time axis (starting at present and going back to the time of origin).

Equations

Full text

Introduction

1

Phylogenetic comparative methods deal with multi-species trait value data. This is an established and rapidly expanding area of research concerning evolution of phenotypes in groups of related species living under various environmental conditions. An important feature of such data is the branching structure of evolution causing dependence among the observed trait values. For this reason the usual starting point for phylogenetic comparative studies is an inferred phylogeny describing the evolutionary relationships. The likelihood can be computed by assuming a model for trait evolution along the branches of this fixed tree, such as the ornstein-Uhlenbeck process.

The one-dimensional Ornstein-Uhlenbeck model is characterized by four parameters: the optimal value θ, the adaptation rate α > 0, the ancestral value X0, and the noise size σ. The classical Brownian motion model [14] can be viewed as a special case with α = 0 and θ being irrelevant. As with any statistical procedure, it is important to be able to compute confidence intervals for these parameters. However, confidence intervals are often not mentioned in phylogenetic comparative studies [8].

There are a number of possible numerical ways of calculating such confidence intervals when the underlying phylogenetic tree is known. Using a regression framework one can apply standard regression theory methods to compute confidence intervals for (θ, X0) conditionally on (α, σ^2^) [15, 20, 28, 33]. Notably in [16] the authors derive analytical formulae for confidence intervals for X0 under the Brownian motion model. In more complicated situations a parametric bootstrap is a (computationally very demanding) way out [8, 11, 27]. Another approach is to report a support surface [20, 21], or consider the curvature of the likelihood surface [7].

All of the above methods have in common that they assume that the phylogeny describing the evolutionary relationships is fully resolved. Possible errors in the topology can cause problems-the closer to the tips they occur, the more problematic they can be [43]. On the other hand, the regression estimators will remain unbiased even with a misspecified tree [34] and also seem to be robust with respect to errors in the phylogeny at least for the Brownian motion model [42]. There are only few papers addressing the issue of phylogenetic uncertainty. An MCMC procedure to jointly estimate the phylogeny and parameters of the Brownian model of trait evolution was suggested in [26, 25]. Recently, [36] develops an Approximate Bayesian Computation framework to estimate Brownian motion parameters in the case of an incomplete tree.

Our paper studies a situation when nothing is known about the phylogeny. The simplest stochastic model addressing this case is a combination of a Yule tree and the Brownian motion on top of it: already in the 1970s, a joint maximum likelihood estimation procedure of a Yule tree and Brownian motion on top of it was proposed in [13]. This basic evolutionary model allows for far reaching analytical analysis [6, 12, 35]. A more realistic stochastic model of this kind combines the Brownian motion with a birth-death tree allowing for extinction of species [10]. For the latter model [35] explicitly compute the so-called interspecies correlation coefficient. Such “tree-free” models are appropriate for working with fossil data when there may be available rich fossilized phenotypic information but the molecular material might have degraded so much that it is impossible to infer evolutionary relationships. In [12] the usefulness of the tree-free approach for contemporary species is demonstrated in an Carnivora order case study and in [31] the distribution over the space of Yule trees of the interspecies correlation coefficient is calculated.

Conditioned birth-death processes as stochastic models for species trees, have received significant attention in the last decade [3, 18, 29, 38, 39, 40]. In this work the unknown tree is modeled by the Yule process conditioned on n extant species while the evolution of a trait along a lineage is viewed as the Ornstein-Uhlenbeck process, see Fig. 1. We study the properties of the sample mean and sample variance computed from the vector of n trait values. Our main results are three asymptotic confidence interval formulae for the optimal trait value θ. These three formulae represent three asymptotic regimes for different values of the adaptation rate α.

In the discussion in [12] it is pointed out that “as evolutionary biologists further refine our knowledge of the tree of life, the number of clades whose phy-logeny is truly unknown may diminish, along with interest in tree-free estimation methods.” In our opinion the main contribution of such methods is that they indicate statistical and asymptotic properties of phylogenetic samples under given evolutionary models. These properties can then be verified for other models of tree growth or real phylogenies [4, 5, 17, 22, 23, 24, 30]. We believe furthermore that the easy-to-compute tree-free predictions will always play an important role of a sanity check to see whether the conclusions based on the inferred phylogeny deviate much from those from a “typical” phylogeny. Moreover, results like those presented here can also be used as a method of testing software for phylogenetic comparative models.

A detailed description of the evolutionary model along with our main results are presented in Section 2. Section 3 contains new formulae for the Laplace transforms of important characteristics of the conditioned Yule species tree: the time to origin U_n_ and the time τ^(n)^ to the most recent common ancestor for a pair of two species chosen at random out of n extant species. In Section 4 we calculate the interspecies correlation coefficient for the Yule-Ornstein-Uhlenbeck model and Section 5 contains the proof of our limit theorems. In Section 6 we establish the consistency of the stationary variance estimator, which is needed for our confidence interval formulae, cf [20] where the residual sum of squares was suggested to estimate the stationary variance. In Appendix A we calculate all the joint moments of U_n_ and τ^(n)^.

Our main result, Theorem 2.1, should be compared with the limit theorems obtained in [1, 2]. They also revealed three asymptotic regimes in a related, though different setting, dealing with a branching Ornstein-Uhlenbeck process. In their case the time of observation is deterministic and the number of the tree tips is random, while in our case the observation time is random and the number of the tips is deterministic. Although it is possible (with some effort) to deduce our results from [1, 2], our proof provides a much more elementary derivation. We believe that our approach will be useful in addressing other biologically relevant issues like the formulae for the higher moments given in Appendix A. Another similar limit theorem, but one conditional on the sequence of species trees generated by different mechanisms, is derived in [5].

On the left: a branching Ornstein-Uhlenbeck process simulated on a realization of the Yule n-tree with n = 5 tips using the TreeSim [38, 39] and mvS-LOUCH [7] R [32] packages. Parameters used are α = 1, σ = 1, X0 −θ = 2, after the tree height Un was scaled to 1. On the right: the species tree disregarding the trait values supplied with the notation for the inter-speciation times. For the pair of tips (2,3) the time τ(n) to their most recent common ancestor is marked on the time axis (starting at present and going back to the time of origin).

The model and main results

2

This work deals with what we call the Yule-Ornstein-Uhlenbeck model which is characterized by four parameters (X_0_, α, σ, θ) and consists of two ingredients

1. the species tree connecting n extant species is modeled by the pure birth Yule process [44] with a unit speciation rate λ = 1 and conditioned on having n tips [18],
2. the observed trait values [eqn] on the tips of the tree evolved from the ancestral state X0 following the Ornstein-Uhlenbeck process with parameters (α, σ, θ).

Definition 2.1 Let (T1,…, Tn) be independent exponential random variables with parameters (1,…, n). We define the Yule n-tree as a random tree with n tips which is constructed using a bottom-up algorithm based on the following two simple rules.

1. During the time period T_k_ the tree has k branches.
2. For k ∈ [2, n] the reduction from k to k − 1 branches occurs as two randomly chosen branches coalesce into one branch.

The height the Yule n-tree is now U_n_ = T1 + … + Tn. As shown in [18], this definition corresponds to the standard Yule tree conditioned on having n tips at the moment of observation, assuming that the time to the origin has the improper uniform prior. Following [11, 20], we model trait evolution along a lineage using the Ornstein-Uhlenbeck process X(t) given by the stochastic differential equation [eqn] Here α > 0 is the adaptation rate, θ is the optimal trait value, σ^2^ is the noise variance, and B(t) is the standard Wiener process. The distribution of X(t) is normal with, [eqn] implying that X(t) looses the effect of the ancestral state X0 at an exponential rate. In the long run the Ornstein-Uhlenbeck process acquires a stationary normal distribution with mean θ and variance σ^2^/2α. We propose asymptotic confidence interval formulae for the optimal value θ which take into account phylogenetic uncertainty. To this end we study properties of the sample mean and sample variance [eqn] Using the properties of the Yule-Ornstein-Uhlenbeck model we find explicit expressions for [eqn], study the asymptotics of [eqn], and prove the following limit theorem revealing three different asymptotic regimes.

Theorem 2.1 Let [eqn] be a normalized difference between the ancestral and optimal values. Consider the normalized sample mean [eqn] of the Yule-Ornstein-Uhlenbeck process with [eqn]. As n → ∞ the process [eqn] has the following limit behavior.

• If α > 0.5, then [eqn] is asymptotically normally distributed with zero mean and variance [eqn].
• If α = 0.5, then [eqn] is asymptotically normally distributed with zero mean and variance 2.
• If α < 0.5, then [eqn] converges a.s. and in L^2^ to a random variable Yα,δ with [eqn].

Let z_x_ be the x-quantile of the standard normal distribution, and q_x_ be the x-quantile of the limit Y_α,δ_. Denote by S_n_ the sample standard deviation defined as the square root of [eqn]. As it will be shown in Section 6, the sample variance [eqn] is a consistent estimator of [eqn]. This fact together with Theorem 2.1 allows us to state the following three approximate (1 − x)-level confidence intervals for θ assuming that we know the value of α: [eqn] Notably, the first of these confidence intervals differs from the classical confidence interval for the mean [eqn] just by a factor K_α_. The latter is larger than 1, as it should, in view of a positive correlation among the sample observations. The correction factor K_α_ becomes negligible in the case of a very strong adaptation, α ⪢ 1, when the dependence due to common ancestry can be neglected.

Remark 2.1 Observe that our standing assumption λ = 1, see Definition 2.1, of having one speciation event per unit of time causes no loss of generality. To incorporate an arbitrary speciation rate λ one has to replace in our formulae parameters α and σ^2^ by α/λ and σ^2^/λ. This transformation corresponds to the time scaling by factor λ in Eq. (1), it changes neither the optimal value θ nor the stationary variance σ^2^/(2α).

Sampling m leaves from the Yule n-tree

3

Here we consider the Yule n-tree, see Definition 2.1 and study some properties of its subtree joining m randomly (without replacement) chosen tips, where m ∈ [2, n]. In particular, we compute the joint Laplace transform of the height of the Yule n-tree U_n_ = T1 + … + T_n_ and τ^(n)^, the height of the most recent common ancestor for two randomly sampled tips, see Fig. 1. For other results concerning the distribution of τ^(n)^ and U_n_ see also [18, 19, 29, 35, 37, 38, 40, 41].

Lemma 3.1 Consider a random m-subtree of the conditioned Yule n-tree. It has m − 1 bifurcating events. Let [eqn] be the consecutive numbers of the bifurcation events in the Yule n-tree (counted from the root toward the leaves) corresponding to the m − 1 bifurcating events of the m-subtree. Put [eqn] and [eqn]. The sequence [eqn] forms a time inhomogeneous Markov chain with transition probabilities [eqn] where [eqn] for all i ≥ 1, and [eqn] PROOF Tracing the lineages of m randomly sampled tips of the Yule n-tree towards the root, the first coalescent event can be viewed as the success in a sequence of independent Bernoulli trials. This argument leads to the expression cf [38] [eqn] confirming the formula stated for the transition probabilities [eqn]. The transition probabilities [eqn] for j = 2,…, m − 1 are obtained similarly. Notice, as a check, that [eqn].

Lemma 3.2 Consider the inter-bifurcation times for the m-subtree of the Yule n-tree [eqn] so that [eqn] for any m ≤ n, and [eqn] Then for x_j_ > − 1 we have [eqn] where k_m_ = n, k0 = 0, and [eqn] PROOF The Laplace transform of the sum of independent exponentials: [eqn] together with Lemma 3.1 implies the stated equality [eqn]

Lemma 3.3 The joint Laplace transform of the height of the Yule n-tree and the height of the most recent common ancestor for two randomly sampled tips is given by [eqn] In particular, [eqn] [eqn] and, denoting the harmonic number h_n_: = 1 + 1/2 + … + 1/n, [eqn] PROOF Turning to Lemma 3.1 with m = 2 we get [eqn] and according to Lemma 3.2 [eqn] This implies the main formula claimed by Lemma 3.3 giving E[e^−U_n_^]= b_n,x_ after putting y = 0. With x = 0, [eqn] When y = 1 this directly becomes [eqn] In the case of y ≠ 1 we use the following relation (easily verified by induction when z ≠ y) [eqn] to derive [eqn]

Lemma 3.4 As n → ∞ for positive x and y we have the following asymptotic results [eqn] where [eqn] PROOF The stated results are obtained from Lemma 3.3 using the first of the following three asymptotic properties of the function b_n,x_ [eqn] These three relations will often be used tacitly in what follows.

Interspecies correlation

4

Denote by [eqn] the σ-algebra containing all information on the Yule n-tree. The scaledtrait values [eqn] with in view of Eq. (2), are conditionally normal with [eqn] which together with the results from Section 3 entails [eqn]

Lemma 4.1 In the framework of the Yule-Ornstein-Uhlenbeck model, for an arbitrary pair of traits we have [eqn] where [eqn] is the backward time to the most recent common ancestor of the tips (i, j). PROOF Denote by [eqn] the normalized trait value of the most recent common ancestor of the tips (i, j). Let [eqn] stand for the σ-algebra generated by [eqn], then using Eq. (2) we get [eqn] implying the statement of this lemma [eqn]

Lemma 4.2 Consider the interspecies correlation coefficient, the unconditioned correlation between two randomly sampled trait values [eqn] If α ≠ 0.5, then [eqn] and in the case of α = 0.5 [eqn] PROOF According to Lemma 4.1 we have, [eqn] leading to [eqn] Applying the results of Section 3 we arrive at the asserted relations for ρ_n_. Observe that asymptotically as n → ∞ the interspecies correlation coefficient decays to 0 as [eqn]

Lemma 4.3 Consider the sample mean [eqn] and the sample variance [eqn] of the scaled trait values. For all α > 0 we have [eqn]. For α ≠ 0.5 [eqn] and in the singular case α = 0.5 [eqn] PROOF Obviously, [eqn]. To prove the other assertions we turn to [35], where the concept of interspecies correlation was originally introduced. It was shown there that the variance of the sample average and the expectation of the sample variance can be compactly expressed as [eqn] Since [eqn], [eqn], and [eqn] it remains to combine Lemma 4.2 with the known expression for [eqn]. A more direct proof of Lemma 4.3 can be obtained using the following result on conditional expectations.

Lemma 4.4 We have [eqn] PROOF The main assertion follows from [eqn]

Proof of Theorem 2.1

5

Lemma 5.1 Put [eqn] with [eqn]. For any x > −1 the sequence [eqn] forms a martingale converging a.s. and in L^2^. Moreover, (U_n_ − log^n^) converges in distribution to a random variable having the standard Gumbel distribution. PROOF The martingale property is obvious [eqn] Since the second moments [eqn] are uniformly bounded over n, we may conclude that [eqn] a.s. and in L^2^ with E[V^(x)^] = 1. It follows that [eqn], and therefore, E[e^−x(U_n_ −log^n^)^] → Γ(x + 1). The latter is a convergence of Laplace transforms confirming the stated convergence in distribution. Observe that the Gumbel limit for U_n_ −logn can be obtained using the classical extreme value theory, in view of the representation [eqn] in terms of independent exponentials with parameter 1. Notice also that Un+1/2 has the same distribution as the total branch length of Kingman’s n-coalescent.

Lemma 5.2 Denote by ℱ_n_ the σ—algebra containing information on the Yule n-tree realization as well as the corresponding information on the evolution of trait values. Set [eqn] The sequence [eqn] forms a martingale with E[Hn] = H0 = δ. PROOF Notice that, [eqn] Hence [eqn]

Lemma 5.3 For all positive α we have [eqn] as n → ∞. PROOF For a given realization of the Yule n-tree we denote by [eqn] and [eqn] two independent versions of τ^(n)^ corresponding to two independent choices of pairs of tips out of n available. We have, [eqn] Writing [eqn] and using the ideas of Section 3 we obtain [eqn] On the other hand, [eqn] Taking the difference between the last two expressions we find [eqn] Using the simple equality [eqn] we see that it suffices to prove that, [eqn] where [eqn] To verify these two asymptotic relations observe that [eqn] Since [eqn], it follows [eqn] and [eqn] PROOF OF THEOREM 2.1 (I) AND (II). Let α > 0.5. To establish the stated normal approximation it is enough to prove the convergence in probability of the first two conditional moments [eqn] since then, due to the conditional normality of [eqn], we will get the following convergence of characteristic functions [eqn] Now, due to Lemma 4.4 we can write [eqn] Using relations from Section 4 we see that [eqn] It remains to observe that on one hand, according to Lemma 5.3 [eqn] and on the other hand, [eqn] implying that [eqn]. This together with μn → 0 holding in L^2^ and therefore in probability, entails [eqn], finishing the proof of part (i). Part (ii) is proven similarly. PROOF OF THEOREM 2.1 (III). Let 0 < α < 0:5. Turning to Lemma 5.2 observe that the martingale [eqn] has uniformly bounded second moments. Indeed, due to Lemma 4.4 [eqn] Thus, according to Lemma 3.4 we have [eqn]. Referring to the martingale L^2^-convergence theorem we conclude that H_n_ → H_∞_ almost surely and in L^2^. Due to Lemma 5.1 it follows that [eqn] Finally, as n → ∞ [eqn]

Consistency of the Sample Variance

6

Recall that [eqn], and according to Lemma 4.3 we have [eqn]. The aim of this section is to show that [eqn] as n → ∞ which is equivalent to [eqn]

To this end we will need the following formula, see Eq. (13) in [9] valid for any normally distributed vector (Z_1_,Z_2_,Z_3_,Z_4_) with means (m1,m2,m3,m4) and covariances Cov [Z_i, Zj] = *c_ij*: [eqn]

In the special case with m_i_ = m it follows [eqn]

Writing Yi instead of [eqn] we use the representation [eqn] to find out that [eqn]

Denoting by (W1, W2, W3, W4) a random sample without replacement of four trait values out of n available, so that [eqn] we derive [eqn]

We compute the five fourth-order moments in the last expression using the conditional normality of the random quadruple (W1, W2, W3, W4) with conditional moments given by [eqn] where [eqn] is the time to the most recent ancestor for the pair of tips (i, j) among m randomly chosen tips of the Yule n-tree. Clearly, all [eqn] have the same distribution as τ^(n)^, and for [eqn] we can find the asymptotics of [eqn] using Lemma 3.4. Notice also that [eqn]

This follows from Lemma 5.3 and Lemma 3.4 as [eqn]

• With Z_1_ = Z_2_ = Z_3_ = Z_4_ = W1 in Eq. (9), we obtain [eqn] and therefore [eqn] as n → ∞.
• Using Eq. (9) with Z_1_ = Z_2_ = Z_3_ = Z_4_ = W1 and Z_4_ = W2 we obtain [eqn] resulting in [eqn] as n → ∞.
• Using Eq. (9) with Z_1_ = Z_2_ = W1 and Z_3_ = Z_4_ = W2 gives [eqn] so that [eqn] as n → ∞.
• Using a consequence of Eq. (9), [eqn] we get [eqn] Using the Cauchy-Schwarz inequality [eqn] we obtain [eqn] as n → ∞.
• According to Eq. (9) we have [eqn] implying [eqn] Using an estimate for [eqn] similar to that we used in (iv), we find E [W1W2W3W4] → 0 as n → ∞.

Finally, putting the above results (i)-(v) into Eq. (10) we arrive at Eq. (8).