DETECTION OF ADAPTIVE SHIFTS ON PHYLOGENIES USING SHIFTED STOCHASTIC PROCESSES ON A TREE

TL;DR

This paper introduces a method to detect adaptive shifts in species traits on phylogenetic trees using statistical models and algorithms.

Contribution

A novel algorithm and statistical framework for identifying adaptive shifts in evolutionary traits while addressing model identifiability issues.

Findings

Models with adaptive shifts are generally not identifiable without constraints.

A recursive algorithm can enumerate equivalent evolutionary scenarios.

A maximum likelihood estimation procedure and model selection method are developed.

Abstract

Comparative and evolutive ecologists are interested in the distribution of quantitative traits among related species. The classical framework for these distributions consists of a random process running along the branches of a phylogenetic tree relating the species. We consider shifts in the process parameters, which reveal fast adaptation to changes of ecological niches. We show that models with shifts are not identifiable in general. Constraining the models to be parsimonious in the number of shifts partially alleviates the problem but several evolutionary scenarios can still provide the same joint distribution for the extant species. We provide a recursive algorithm to enumerate all the equivalent scenarios and to count the effectively different scenarios. We introduce an incomplete-data framework and develop a maximum likelihood estimation procedure based on the EM algorithm. Finally, we propose a model selection procedure, based on the cardinal of effective scenarios, to estimate the number of shifts and prove an oracle inequality.

Figures

Trait evolution under a Brownian Motion. The ancestral value of the trait is 0, and the observed values (time 800) range from −4 to 11 for extant species. One shift occurs on the parent branch of (D,E), changing the trajectory of their ancestral trait value from the grey one to the colored one. The shift increases the observed dispersion.

A rooted and time calibrated phylogenetic tree with the notations used to parametrize the tree (l,t,d,b) and the observed (Y) and non-observed (Z) variables.

Equivalent allocations in the BM case. Mean tip values are represented by colors and equal for all allocations. The two allocations on the right are parsimonious.

Five equivalent shift allocations that produce colorings that are adapted to an OUsun, with α = 3 and γ2 = 0.1. The box at the root represents the ancestral optimum β1, and the boxes on the branches represent the positions and values of the shifts on the optimal value. While accounting for very different evolutionary scenarios, all allocations produce the same trait distribution at the tips.

Partitions in 2 groups.

Left: Simulated configuration (with t1/2 = 0.75, γ2 = 0.5 and K = 5). The shifts positions and values are marked on the tree. The value of the character generated (positive or negative) is represented on the right. The colors of the branches correspond to the true regimes, black being the ancestral state. Right: One of the three equivalent allocations of shifts for the model inferred from the data, with corresponding vector of mean tip values. Shifts not recovered are located on pendant edges, and have low influence on the data. The two other equivalent allocations can be easily deducted from this one.

Box plots over the 200 repetitions of each set of parameters, for the phylogenetic halflife (top) and root variance (bottom) with K estimated, and α fixed to its true value (blue) or estimated (red), on a tree with 128 taxa. For better legibility, the y-axis of these two rows were re-scaled, omitting some outliers (respectively, for t1/2 and γ2, 0.82% and 0.46% of points are omitted). The whisker of the first box for t1/2 goes up to 7.5.

Same for β1(top), the number of shifts (middle) and ARI (bottom), with a estimated, and K fixed to its true value (blue) or estimated (red). As previously, the y-axis of the first row (β1) was re-scaled, omitting some outliers (1.38% of points are omitted).

Phylogenetic tree of the Chelonians. Log-transformed traits values are represented on the right. Branch colors represent the habitats. The shifts found by our EM algorithm are shown as circles, with a color indicating the value of the shift, from blue (negative) to red (positive). Boxes highlight the groups induced by the shifts. The x-scale is in million years.

Tables

Equations

Full text

Introduction

1

Motivations: Environmental Shifts

1.1

An important goal of comparative and evolutionary biology is to decipher the past evolutionary mechanisms that shaped present day diversity, and more specifically to detect the dramatic changes that occurred in the past (see for instance Losos, 1990; Mahler et al., 2013; Davis et al., 2007; Jaffe et al., 2011). It is well established that related organisms do not evolve independently (Felsenstein, 1985): their shared evolutionary history is well represented by a phylogenetic tree. In order to explain the pattern of traits measured on a set of related species, one needs to take these correlations into account. Indeed, a given species will be more likely to have a similar trait value to her “sister” (a closely related species) than to her “cousin” (a distantly related species), just because of the structure of the tree. On top of that structure, when considering a functional trait (i.e. a trait directly linked to the fitness of its bearer), such as shell size for turtles (Jaffe et al., 2011), one needs to take into account the effect of the species environment on its traits. Indeed, a change in the environment for a subset of species, like a move to the Galàpagos Islands for turtles, will affect the observed trait distribution, here with a shift towards giant shell sizes compared to mainland turtles. The observed present-day trait distribution hence contains the footprint of adaptive events, and should allow us to detect unobserved past events, like the migration of one ancestral species to a new environment. Our goal here is to devise a statistical method based on a rigorous maximum likelihood framework to automatically detect the past environmental shifts that shaped the present day trait distribution.

optimum. The value will then be passed down to its descendant, potentially leading to unexpectedly large differences between sub-families of species. The distribution of trait values across extant species hence contains the footprint of adaptive events and should in principle allow us to detect unobserved past events. Our goal here is to devise a statistical framework to automatically detect the past environmental shifts that shaped the present day trait distribution.

Stochastic Process on a tree

1.2

We model the evolution of a quantitative adaptive trait using the framework of stochastic processes on a tree. Specifically, given a rooted phylogenetic tree, we assume that the trait evolves according to a given stochastic process on each branch of the tree. At each speciation event, or equivalently node of the tree, one independent copy with the same initial conditions and parameters is created for each daughter species, or outgoing branches.

Tree Structure. This model is our null model: it accounts for the tree-induced distribution of trait values in the absence of shifts. Depending on the phenomenon studied, several stochastic processes can be used to capture the dynamic of the trait evolution. In the following, we will use the Brownian Motion (BM) and the Ornstein-Uhlenbeck (OU) processes.

Trait evolution under a Brownian Motion. The ancestral value of the trait is 0, and the observed values (time 800) range from −4 to 11 for extant species. One shift occurs on the parent branch of (D,E), changing the trajectory of their ancestral trait value from the grey one to the colored one. The shift increases the observed dispersion.

Brownian Motion. Since the seminal article of Felsenstein (1985), the BM has been used as a neutral model of trait evolution. If (B_t_; t ≥ 0) is the Brownian motion, a character (W_t_;t ≥ 0) evolves on a lineage according the the stochastic differential equation: dW_t_ = σdB_t_, σ^2^ being a variance parameter. If μ is the ancestral value at the root of the tree (t = 0), then [eqn]. The variance σ^2^t of the trait is proportional to the time of evolution and the covariance σ^2^t_ij_ between two species i and j is proportional to their time of shared evolution.

Ornstein-Uhlenbeck Process. An unbounded variance is quite unrealistic for adaptive traits (Butler and King, 2004). For that reason, the OU process, that models stabilizing selection around an adaptive optimum (Hansen, 1997) is usually preferred to the BM. It is defined by the stochastic differential equation dW_t_ = −α(W_t_ − β)dt + σdB_t_, and has stationary distribution [eqn]. In this equation, W_t_ is the secondary optimum of a species, a trade-off between all selective constraints - e.g. ecological - on the trait and can be approached by the population mean of that species. The term −α(W_t_ − β)dt of the equation represents the effects of stabilizing selection towards a primary optimum β, that depends only on the ecological niche of the species. The selection strength is controlled by the call-back parameter α. For interpretation purpose, we will use the phylogenetic half-life t1/2 = ln(2)/α, defined as the time needed for the expected trait value to move half the distance from the ancestral state to the primary optimum (Hansen, 1997). The term σdB_t_ represents the random effects of uncontrolled factors, ranging from genetic drift to environmental fluctuations. We refer to Hansen (1997); Hansen et al. (2008) for further discussion and deeper biological interpretations of the hypothesis underlying this model of evolution. The aim of our work is to detect environmental shifts.

Environmental Shifts. In addition to the previous mechanisms, we assume that abrupt environmental changes affected the ecological niche of some ancestral species. We model these changes as instantaneous shifts in the parameters of the stochastic process. Shifted parameters are inherited along time and thus naturally create clusters of extant species that share the same parameters trajectories. In the BM process, shifts affect the mean value of the trait and are thus instantaneously passed on to the trait itself (see Figure 1) whereas in the OU process, shifts affect the primary optimum β. In this case, the trait converges to its new stationary value with an exponential decay of half-life t_1/2_ inducing a lag that makes recent shifts harder to detect (Hansen and Bartoszek, 2012). In the remainder, we assume that all other parameters (σ^2^ for the BM and σ^2^, α for the OU) are fixed and constant (but see Beaulieu et al., 2012; Rabosky, 2014, for partial relaxations of this hypothesis).

Scope of this article

1.3

State of the Art. Phylogenetics Comparative Methods (PCM) are an active field that has seen many fruitful developments in the last few years (see Pennell and Harmon, 2013, for an extensive review). Several methods have been specifically developed to study adaptive evolution, starting with the work of Butler and King (2004). Butler and King (2004) only consider shifts in the optimal value β whereas Beaulieu et al. (2012) also allow for shifts in the selection strength α and the variance σ^2^. Both have in common that shift locations are assumed to be known. Several extensions of the model without or with known shifts have also been proposed: Hansen et al. (2008) extended the original work of Hansen (1997) on OU processes to a two-tiered model where β(t) is itself a stochastic process (either BM or OU). Bartoszek et al. (2012) extended it further to multivariate traits whereas Hansen and Bartoszek (2012) introduced errors in the observations. Expanding upon the BM, (Landis et al., 2013) replaced fixed shifts, known or unknown, by random jump processes using Levy processes. Non-Gaussian models of trait evolution were also recently considered by Hiscott et al. (2015), who adapted Felsenstein’s pruning algorithm for the likelihood computation of these models, using efficient integration techniques. Finally, Ho and Ané (2013) derived consistency results for estimation of the parameters of an OU on a tree and Bartoszek and Sagitov (2015); Sagitov and Bartoszek (2012) computed confidence intervals of the same parameters by assuming an unknown random tree topology and averaging over it.

The first steps toward automatic detection of shifts, which is the problem of interest in this paper, have been done in a Bayesian framework, for both the BM (Eastman et al., 2013) and the OU (Uyeda and Harmon, 2014). Using RJ-MCMC, they provide the user with the posterior distribution of the number and location of shifts on the tree. Convergence is however severely hampered by the size of the search space. The growing use of PCM in fields where large trees are the norm makes maximum likelihood based point estimates of the shift locations more practical. A stepwise selection procedure for the shifts has been proposed in Ingram and Mahler (2013). The procedure adds shifts one at the time and is therefore rather efficient but the selection criterion is heuristic and has no theoretical grounding for that problem, where observations are correlated through the tree structure. These limitations have been pointed out in Ho and Ané (2014b). In this article, the authors describe several identifiability problems that arise when trying to infer the position of the shifts on a tree, and propose a different stepwise algorithm based on a more stringent selection criterion, heuristically inspired by segmentation algorithms.

To rigorously tackle this issue, we introduced a framework where a univariate trait evolves according to an OU process with stationary root state (S) on an Ultrametric tree (U). Furthermore, as the exact position of a shift on a branch is not identifiable for an ultrametric tree, we assume that shifts are concomitant to speciation events and only occur at Nodes (N) of the tree. We refer to this model as OUsun hereafter.

Our contribution. In this work, we make several major contributions to the problem at hand. First, we derive a statistical method to find a maximum likelihood estimate of the parameters of the model. When the number of shifts is fixed, we work out an Expectation Maximization (EM) algorithm that takes advantage of the tree structure of the data to efficiently maximize the likelihood. Second, we show that, given the model used and the kind of data available, some evolutionary scenarios remain indistinguishable. Formally, we exhibit some identifiability problems in the location of the shifts, even when their number is fixed, and subsequently give a precise characterization of the space of models that can be inferred from the data on extant species. Third, we provide a rigorous model selection criterion to choose the number of shifts needed to best explain the data. Thanks to our knowledge of the structure of the spaces of models, acquired through our identifiability study, we are able to mathematically derive a penalization term, together with an oracle inequality on the estimator found. Fourth and finally, we implement the method on the statistical software R (R Core Team, 2014), and show that it correctly recovers the structure of the model on simulated datasets. When applied to a biological example, it gives results that are easily interpretable, and coherent with previously developed methods. All the code used in this article is publicly available on GitHub (https://github.com/pbastide/PhylogeneticEM).

Outline. In Section 2, we present the model, using two different mathematical point of views, that are both useful in different aspects of the inference. In Section 3, we tackle the identifiability problems associated with this model, and describe efficient algorithms to enumerate, first, all equivalent models within a class, and, second, the number of truly different models for a given number of shifts. These two sections form the foundation of Section 4, in which we describe our fully integrated maximum likelihood inference procedure. Finally, in Sections 5 and 6, we conduct some numerical experiments on simulated and biological datasets.

Statistical Modeling

2

Probabilistic Model

2.1

A rooted and time calibrated phylogenetic tree with the notations used to parametrize the tree (l,t,d,b) and the observed (Y) and non-observed (Z) variables.

Tree Parametrization. As shown in Figure 2, we consider a rooted tree [eqn] with n tips and m internal nodes (m = n − 1 for binary trees). The internal nodes are numbered from 1 (the root) to m, and the tips from m + 1 to m + n. Let i be an integer, [eqn]. Then pa(i) is the unique parent of node i. The branch leading to i from pa(i) is noted b_i_ and has length ℓi = t_i_ − tpa(i) where t_i_ is the time elapsed between the root and node i. By convention, we set t1 = 0 and tpa(1) = ∞ for the root. The last convention ensures that the trait follows the stationary distribution (if any) of the process at the root. We denote Anc(i) = {pa^r^(i): r ≥ 0} the set composed of node i and of all its ancestors up to the root. For a couple of integers (i,j), [eqn], nodes i and j are at phylogenetic distance d_ij_ and the time of their most recent common ancestor (mrca) is t_ij_. We consider ultrametric trees, for which t_m+1_ = …= t_m*+n_* =: h and note h the tree height. In the following, the tree is fixed and assumed to be known.

Trait Values. We denote by X the vector of size m + n of the trait values at the nodes of the tree. We split this vector between non-observed values Z (size m) at the internal nodes, and observed values Y (size n) at the tips, so that X^T^ = (Z^T^, Y^T^). According to our model of trait evolution, the random variable [eqn], is the result of a stochastic process stopped at time t_i_. In the following, we assume that the inference in the BM case is done conditionally to a fixed root value X1 = μ. In the OUsun case, we assume that the root trait value is randomly drawn from the stationary distribution: [eqn], where β1 is the ancestral optimal value.

Shifts. We assume that K shifts occur on the tree, [eqn]. The k^th^ shift, [eqn], occurs at the beginning of branch τ_k_, [eqn], and has intensity δ_k_, δ_k_ ∈ ℝ. The interpretation of this intensity depends on the process. In the following, we use the vector Δ of shifts on the branches, of size m + n, with K + 1 non-zero entries, and defined as follows (see example 2.1):

$$[eqn]$$

Note that no proper shift occurs on the root branch, but that the root trait value or mean, μ, is formalized as an initial fictive shift on this fictive branch.

Parameters. The parameters needed to describe an OUsun (respectively, a BM) are θ = (γ, α,Δ) (resp. θ = (σ, Δ)). Note that, as σ^2^ = 2αγ^2^, only the two parameters α and γ are needed to describe the OUsun. We denote by OUsun(θ) (resp. BM(6)) the OUsun (resp. BM) process running on the tree with parameters θ.

Incomplete Data Model Point of View

2.2

If the trait values were observed at all nodes of the tree, including ancestral ones, shifts would be characterized by unexpectedly large differences between a node and its parent. A way to mimic this favorable case is to use an incomplete data model, as described below. This representation of the model will be useful for the parametric inference using an EM algorithm (Section 4.1).

Brownian Motion. As the shifts occur directly in the mean of the process, we get: [eqn]

The trait value at node i, [eqn], is centered on the value of its parent node Xpa(i), with a variance proportional to the evolution time ℓ_i_ between i and pa(i). The effect of a non-zero shift Δ_i_ on branch b_i_ is simply to translate the trait value by Δ_i_.

Ornstein-Uhlenbeck. The shifts occur on the primary optimum β, which is piecewise constant. As the shifts are assumed to occur at nodes, the primary optimum is entirely defined by its initial value β_1_ and its values β2,…, βn+m on branches of the tree, where β_i_ is the value on branch b_i_ leading to node i.

$$[eqn]$$

Assuming that the root node is in the stationary state, we get:

$$[eqn]$$

The trait value at node i depends on both the trait value at the father node Xpa(i) and the value β_i_ of the primary optimum on branch b_i_. Contrary to the BM case, the shifts only appear indirectly in the distributions of X_i_s, through the values of β, and with a shrinkage of 1 − e^−αd^ for shifts of age d, which makes recent shifts (d small compared to 1/α) harder to detect.

Linear Regression Model Point of View

2.3

A more compact and direct representation of the model is to use the tree incidence matrix to link linearly the observed values (at the tips) with the shift values, as explained below. We will use this linear regression framework for the Lasso (Tibshirani, 1996) initialization of the EM (Section 4.1) and the model selection procedure (Section 4.2). It will also help us to explore identifiability issues raised in the next section.

Matrix of a tree. It follows from the recursive definition of X that it is a Gaussian vector. In order to express its mean vector given the shifts, we introduce the tree squared matrix U, of size (m + n), defined by its general term: [eqn]. In other words, the j^th^ column of this matrix, [eqn], is the indicator vector of the descendants of node j. To express the mean vector of the observed values Y, we also need the sub-matrix T, of size n × (m + n), composed of the bottom n rows of matrix U, corresponding to the tips (see example 2.1 below). Likewise, the i^th^ row of T, [eqn], is the indicator vector of the ancestors of leaf m + i.

Brownian Motion. From the tree structure, we get:

$$[eqn]$$

Here, [eqn] is a Gaussian error vector with co-variances [∑XX]ij = σ^2^t_ij_ for any 1 ≤ i,j ≤ m + n, and EY is the vector made of the last n coordinates of E_x_.

Ornstein-Uhlenbeck. For the OUsun, shifts occur on the primary optimum, and there is a lag term, so that: [eqn] where [eqn] and [eqn] are diagonal matrices of size m + n. As previously, [eqn], but [eqn]. As the tree is ultrametric, this expression simplifies to the following one when considering only observed values: [eqn] where Ey is the Gaussian vector made of the last n coordinates of EX, and [eqn] is a diagonal matrix of size m + n. Note that if α is positive, then α(h − tpa(i) > 0 for any [eqn], and W(α) is invertible.

Example 2.1. The tree presented in Figure 2 has five tips and one shift on branch 4 + 3 = 7, so:

$$[eqn]$$

And, respectively, for a BM or an OUsun (μ =β1): [eqn] [eqn]

Space of Expectations. Expressions 2.4 and 2.6 allow us to link the parameter θ to the probability distribution of observations Y and to explore identifiability issues. In this linear formulation, detecting shifts boils down to identifying the non-zero components of Δ. The following lemma highlights the parallels between solutions of the BM and OUsun processes:

Lemma 2.1 (Similar Solutions). Let m_Y_ ∈ ℝ^n^ be a vector, [eqn] an ultrametric tree, α a positive real number, and α, γ non-negative real numbers. Then there exists at least one vector Δ^BM^, Δ^BM^ ∈ ℝ^m+n^ (respectively, Δ^OU^ ∈ ℝ^m+n^), such that the vector of expectations at the tips of a BM(σ, Δ^BM^) (respectively, an OUsun(γ,α, Δ^OU^)) running on the tree [eqn] is exactly mY. Furthermore, Δ^BM^ is a solution to this problem for the BM if and only if Δ^OU^ = W(α)^−1^ Δ^BM^ is a solution for the OUsun, and Δ^BM^ and W(α)^−1^Δ^BM^ have the same support. These two vectors are said to be similar.

Proof. The first part of this lemma follows directly from formulas 2.4 (BM) and 2.6 (OU). Indeed, the maps Δ ↦ TΔ and Δ ↦ TW(α)Δ both span ℝ^n^. The second part of the lemma is a consequence of W(α) being diagonal and invertible (for α > 0).

Remark 2.1. Lemma 2.1 shows that the OUsun and BM processes that induce a given mY use shifts located on the same branches, although they may differ on other parameters.

Identifiability and Complexity of a Model

3

Identifiability Issues

3.1

As we only have access to Y, and not X, we only have partial information about the shifts occurrence on the tree. In fact, several different allocations of the shifts can produce the same trait distribution at the tips, and hence are not identifiable. In other words, there exists parameters θ ≠ θ′ with the same likelihood function: p_θ_(·) = p_θ′_(·). Note that the notion of identifiability is intrinsic to the model and affects all estimation methods. Restricting ourselves to the parsimonious allocations of shifts only partially alleviates this issue, and, using a “random cluster model” representation of the problem, we are able to enumerate, first, all the equivalent solutions to a given problem, and, second, all the equivalence classes for a given number of shifts.

Equivalent allocations in the BM case. Mean tip values are represented by colors and equal for all allocations. The two allocations on the right are parsimonious.

No Homoplasy Assumption. We assume in the following that there is no convergent evolution. This means that each shift creates a new (and unique) mean trait value for extant species that are below it. This assumption is reasonable considering that shifts are real valued and makes the model similar to “infinite alleles” models in population genetics. This assumption confines but does not eliminate the identifiability issue, as seen in Figure 3.

Definition of the problem

3.1.1

Figure 3 shows a simple example where the model is not identifiable in the BM case. Here, four distinct allocations give the same mean values (μ + δ_1_,μ + δ2) at the tips. The lack of identifiability is due to the non-invertibility of the tree matrix T.

Proposition 3.1 (Kernel of the Tree Matrix T). Let i be an internal node, [eqn], with L_i_ children nodes [eqn]. Then the vector K^i^ defined as follow: [eqn] is in the kernel of T. In addition, the m vectors constructed this way form a basis of the kernel space of T.

These kernel vectors effectively “cancel out” a shift on a branch by balancing it with the opposite shift on all immediate child branches. Note that the root mean value is treated as a shift. The following lemma describes the relationships that exist between these kernel vectors and the tree matrix U defined in Section 2.3. The technical proofs of Proposition 3.1 and Lemma 3.1 are postponed to Appendix A.

Lemma 3.1. Let b be the canonical basis of ℝ^m+n^, and S a supplementary space of ker(T). Then b′ = (K^1^, …, K^m^, bm+1, …, bm+n) is a basis adapted to the decomposition ker(T) ⊕ S, and the matrix U (as defined in Section 2.3) is the change of basis matrix between b and b′. As a consequence, U is invertible.

“Random Cluster Model” Representation. When inferring the shifts, we have to keep in mind this problem of non-identifiability, and be able to choose, if necessary, one or several possible allocations among all the equivalent ones. In order to study the properties of the allocations, we use a random cluster model, as defined in Mossel and Steel (2004). The following definition states the problem as a node coloring problem.

Definition 3.1 (Node Coloring). Let [eqn] be a set of K arbitrary “colors”, K ∈ ℕ*. For a given shift allocation, the color of each node is given by the application [eqn] recursively defined in the following way:

• Choose a color [eqn] for the root: B(1) = c.
• For a node i, [eqn], set B(i) to B(pa(i)) if there is no shift on branch i, otherwise choose another color c, [eqn], and set B(i) to c.

Hereafter, we identify [eqn] with [eqn] and refer to a node coloring indifferently as an application or a vector.

As the shifts only affect [eqn] and we only have access to [eqn], we identify colors with the distinct values of [eqn]:

Definition 3.2 (Adapted Node Coloring). A node coloring is said to be adapted to a shifted random process on a tree if two tips have the same color if and only if they have the same mean value under that process.

Proposition 3.2 (Adapted Coloring for BM and OUsun). Let σ and γ be two non-negative real numbers, and α a positive real number. Then:

• In the BM case, if [eqn] is the set of possible mean values taken by the nodes of the tree, then the knowledge of the node colors is equivalent to the knowledge of Δ. Furthermore, the associated node coloring is adapted to the original BM.
• In the OUsun case, from lemma 2.1, we can find a similar BM process, i.e. with shifts on the same branches. Then the knowledge of the node coloring associated to this similar BM process is equivalent to the knowledge of the vector of shifts of the OUsun, and the node coloring obtained is adapted to the original OUsun.

Proof of Proposition 3.2. The proof of (i) relies on expression 2.4, that states that [eqn]. Defining [eqn] as the set off all distinct values of [eqn], we can identify [eqn] with the node coloring application that maps any node i with [eqn]. Since U is invertible (see lemma 3.1 above), we can go from one formalism to the other.

For (ii), we use lemma 2.1 to find a similar BM, and then use (i).

From now on, we will study the problem of shifts allocation as a discrete-state coloring problem.

Parsimony

3.1.2

As we saw on Figure 3 there are multiple colorings of the internal nodes that lead to a given tips coloring. Among all these solutions, we choose to study only the parsimonious ones. This property can be seen as an optimality condition, as defined below:

Definition 3.3 (Parsimonious Allocation). Given a vector of mean values at the tips produced by a given shifted stochastic process running on the tree, an adapted node coloring is said to be parsimonious if it has a minimum number of color changes. We denote by [eqn] the set of parsimonious allocations of K shifts on the (m + n − 1) branches of the tree (not counting the root branch). As K shifts cannot produce more than K + 1 colors, we can define an application [eqn] that maps a parsimonious allocation of shifts to its associated tip partition.

Definition 3.4 (Equivalence). Two allocations are said to be equivalent (noted ~) if they produce the same partition of the tips and are both parsimonious. Mathematically: [eqn] In other words, two allocations are equivalent if they produce the same tip coloring up to a permutation of the colors. Given [eqn] a coloring of the tips of [eqn] with K + 1 colors, ϕ^−1^(d) is the set of equivalent parsimonious node coloring that coincide with d (up to a permutation of the colors) on the tree leaves. Several dynamic programming algorithms already exist to compute the minimal number of shifts required to produce a given tips coloring, and to find one associated parsimonious solution (see Fitch, 1971; Sankoff, 1975; Felsenstein, 2004). Here, we need to be a little more precise, as we want to both count and enumerate all possible equivalent node colorings associated with a tip coloring. For the sake of brevity, we only present the algorithm that counts |ϕ^−1^(d)|, for [eqn]. This algorithm can be seen as a corollary of the enumeration algorithm (presented and proved in Appendix B) and an extension of Fitch algorithm where we keep track of both the cost of an optimal coloring and the number of such colorings. It has O(K^2^Ln) time complexity where L is the maximal number of children of the nodes of the tree.

Proposition 3.3 (Size of an equivalence class). Let d be a coloring of the tips,[eqn], and let i be a node of tree [eqn] with L_i_ daughter nodes (i1, …, i*L_i), Li*_ ≥ 2. Denote by [eqn] the sub-tree rooted at node i.

For [eqn], Si(k) is the cost of starting from node i with color k, i.e. the minimal number of shifts needed to get the coloring of the tips of [eqn] defined by d, when starting with node i in color k. Denote by T_i_(k) the number of allocations on [eqn] that achieve cost S_i_(k).

If i is a tip (m + 1 ≤ i ≤ m + n), then, [eqn]

Otherwise, if i is a node, for 1 ≤ l ≤ L_i_, define the set of admissible colors for daughter i_l_: [eqn]

As these sets are not empty, let [eqn]. Then: [eqn]

At the root, if [eqn], then [eqn].

OU Practical Case. We can illustrate this notion on a simple example. We consider an OUsun on a random tree of unit height (total height h = 1). We put three shifts on the tree, producing a given trait distribution. Then, using proposition 3.2 and our enumeration algorithm, we can reconstruct the 5 possible allocations of shifts that produce the exact same distribution at the tips. These solutions are shown in Figure 4. Note that the colors are not defined by the values of the optimal regime β, but by the mean values 𝔼 [Y] of the process at the tips. As a result, the groups shown in blue and red in the first solution have the same optimal value in this configuration, but not in any other. The second solution shown illustrates the fact that all the shifts values are inter-dependent, as changing the position of only one of them can have repercussions on all the others. Finally, the third solution shows that the timing of shifts matters: to have the same impact as an old shift, a recent one must have a much higher intensity (under constant selection strength such as in the OUsun). preferring cases with lower shifts values. According to this rule, solutions 2 or 4 would be preferred.

Five equivalent shift allocations that produce colorings that are adapted to an OUsun, with α = 3 and γ2 = 0.1. The box at the root represents the ancestral optimum β1, and the boxes on the branches represent the positions and values of the shifts on the optimal value. While accounting for very different evolutionary scenarios, all allocations produce the same trait distribution at the tips.

Possible relaxation of the No-homoplasy assumption. Note that the algorithms used for counting and enumerating the configurations of an equivalence class are valid even without the no-homoplasy hypothesis. The no-homoplasy hypothesis is however crucial in the next Section to establish a link between the number of shifts and the number of distinct tips colors.

Complexity of a Collection of Models

3.2

Number of different tips colors. As we make the inference on the parameters with a fixed number of shifts K (see Section 4.1), we need a model selection procedure to choose K. This procedure depends on the complexity of the collection of models that use K shifts, defined as the number of distinct models. To do that, we count the number of tree-compatible partitions of the tips into K + 1 groups, as defined in the next proposition:

Proposition 3.4. Under the no homoplasy assumption, an allocation of K shifts on a tree is parsimonious if and only if it creates exactly K + 1 tip colors. The tip partition into K + 1 groups associated with this coloring is said to be tree-compatible. The set [eqn] of such partitions is the image of [eqn] by the map ϕ defined in the previous section.

Proof of Proposition 3.4. First, note that K shifts create at most K + 1 colors. If each shift produces a new tip mean value (no homoplasy), the only way to create K or less colors is to “forget” one of the shifts, i.e. to put shifts on every descendant of the branch where it happens. Such an allocation is not parsimonious, as we could just add the value of the forgotten shift to all its descendant to get the same coloring of the tips with one less shift. So a parsimonious allocation cannot create less than K + 1 colors, and hence creates exactly K + 1 colors.

Reciprocally, if an allocation with K shifts that produces p groups is not parsimonious, then we can find another parsimonious one that produces the same p groups with p − 1 shifts, with p − 1 < K, i.e. p < K + 1. So, by contraposition, if the allocation produces K + 1 groups, then it is parsimonious.

Using the equivalence relation defined in Definition 3.4, we can formally take the quotient set of [eqn] by the relation ~ to get the set of parsimonious allocations of K shifts on the m + n − 1 branches of the tree that are identifiable: [eqn]. In other words, the set [eqn] is constituted of one representative of each equivalence class. Under the no homoplasy assumption, there is thus a bijection between identifiable parsimonious allocations of K shifts and tree-compatible partitions of the tips in K + 1 groups: [eqn].

The number [eqn] is the complexity of the class of models with K shifts defined as the number of distinct identifiable parsimonious possible configurations one can get with K shifts on the tree. To compute [eqn], we will need [eqn] the number of marked tree-compatible partitions in K groups. These are composed of all the tree-compatible partitions where one group, among those that could be in the same state as the root, is distinguished with a mark (see example 3.1 below).

Example 3.1 (Difference between [eqn] and [eqn].

• If we consider only unmarked partitions, then colorings 1, 2 and 3 induce the same partitions as, respectively, colorings 4, 5 and 6, and [eqn].
• For marked partitions, fix the root state to an arbitrary color, for instance white, and consider the white group as marked. Then colorings 5 and 6 are not tree-compatible (they require two shifts). And although they induce the same partition, colorings 1 and 4 correspond to different marked partitions: each marks a different group of leaves. Therefore [eqn].

Partitions in 2 groups.

Proposition 3.5 (Computation of the Number of Equivalent Classes). Let i be a node of tree [eqn], and K ∈ ℕ*.

If i is a tip, then [eqn].

Else, if i is a node with Li daughter nodes (i1, …, i*L_i), *L_i ≥ 2, then: [eqn]

In the binary case, this relation becomes, if i has two daughters iℓ and ir: [eqn]

Proof. We will prove this proposition in the binary case, the general case being a natural extension of it. If [eqn] is a binary tree with [eqn] and [eqn] as left and right sub-trees, one faces two situations when partitioning the tips in K groups:

• The left and right sub-trees have no group in common. Then, the number of groups in [eqn] is equal to the number of groups in its two sub-trees, and there are [eqn] such partitions. This is the first term of the equation on [eqn] in 3.2.
• The left and right sub-trees have at least one group in common. Then, from the no homoplasy assumption, they have exactly one group in common: the ancestral state of the root. Suppose that this ancestral state is marked. Then it must be present in the two sub-trees, and there are [eqn] such partitions. This ends the proof of the formula on [eqn].

To get the formula on [eqn], we use the same kind of arguments. The second part of the formula is the same as the one for [eqn], and the first part corresponds to trees for which the marked partition is present in only one of the two sub-trees.

The complexity of the algorithm described above is O(2^L^(K + L)^L^Ln). Note that [eqn] depends on the topology of the tree [eqn] in general. However, if the tree is binary, a closed form solution of the recurrence relation 3.1, which does not depend on the topology, exists.

Corollary 3.1 (Closed Formula Binary Trees). For a rooted binary tree with n tips, we have: [eqn]

The demonstration of this formula is not straightforward, and is based on a Vandermonde-like equality, detailed in Appendix C. The formula is then obtained using a strong induction on the number of tips of the tree.

Remark 3.1. Note that, when K is large compared to [eqn], the average number of configurations per equivalence class goes to infinity. This can be checked by comparing the total number of configurations [eqn] with the total number of classes [eqn]. As a consequence, we only consider models for which [eqn] in the remainder.

Remark 3.2. This formula was already obtained in a different context in Steel (1992) (Proposition 1) and, with a slightly different formulation, in Semple and Steel (2003, Proposition 4.1.4). In these works, the authors are interested in counting the “r-states convex characters on a binary tree”. Under the no-homoplasy assumption, this number can be shown to be equal to [eqn].

Another Characterization of Parsimony

3.3

The following proposition gives an alternative definition of parsimony under the no-homoplasy hypothesis using the linear formulation of the problem. It will be used for model selection in Section 4.2. Its technical proof is postponed to Appendix A.

Proposition 3.6 (Equivalence between parsimony and independence). Let mY be a given mean vector, mY ∈ ℝ^n^, and Δ a vector of shifts such that TΔ = mY, with T the tree matrix defined in Section 2.3. Under the no homoplasy assumption, the vector of shifts Δ is parsimonious if and only if the corresponding column-vectors of the tree matrix [eqn] are linearly independent.

Statistical Inference

4

Expectation Maximization

4.1

Principle. As shown in Section 2.2, both BM and OUsun models can be seen as incomplete data models. The Expectation Maximization algorithm (EM, Dempster et al., 1977) is a widely used algorithm for likelihood maximization of these kinds of models. It is based on the decomposition: [eqn]. Given an estimate θ^(h)^ of the parameters, we need to compute some moments of p_θ^(h)^_(Z | Y) (E step), and then find a new estimate [eqn]. The parameters are given for the BM and OUsun in subsection 2.1. We assume here that the number of shifts K is fixed. We only provide the main steps of the EM. Additional details can be found in Appendix D.

E step. As X is Gaussian, the law of the hidden variables Z knowing the observed variables Y is entirely defined by its expectation and variance-covariance matrix, and can be computed using classical formulas for Gaussian conditioning. The needed moments of Z | Y can also be computed using a procedure that is linear in the number of tips (called “Upward-downward”) that takes advantage of the tree structure and bypasses inversion of the variance-covariance matrix (see Lartillot, 2014, for a similar algorithm).

Complete Likelihood Computation. Using the model described in Section 2.2, we can use the following decomposition of the complete likelihood: [eqn]

Each term of this product is then known, and we easily get [eqn].

M step. The difficulty comes here from the discrete variables (location of shifts on the branches). The maximization is exact for the BM but we only raise the objective function for the OUsun, hence computing a Generalized EM (GEM, seeDempster et al., 1977). This stems from the independent increment nature of the BM: shifts only affect pθ (X_j_ | Xpa(j))) on the branches where they occur and the maximization reduces to finding the K highest components of a vector, which has complexity O(n + K log(n)). By contrast, OUsun has autocorrelated increments: shifts affect pθ (X_j_ | Xpa(j)) on the branches where they occur and on all subsequent branches. Maximization is therefore akin to segmentation on a tree, which has complexity O(n^K^).

Initialization. Initialization is always a crucial step when using an EM algorithm. Here, we use the linear formulation 2.4 or 2.6, and initialize the vector of shifts using a Lasso regression. The selection strength a is initialized using pairs of tips likely to be in the same group.

Model Selection

4.2

Model Selection in the iid Case with Unknown Variance. Model selection in a linear regression setting has received a lot of attention over the last few years. In Baraud et al. (2009), the authors developed a non-asymptotic method for model selection in the case where the errors are independent and identically distributed (iid), with an unknown variance. In the following, we first recall their main results, and then adapt it to our setting of non-independent errors.

We assume that we have the following model of independent observations: [eqn] and we define a collection [eqn] of linear subspaces of ℝ^n^ that we call models, and that are indexed by a finite or countable set ℳ. For each η ∈ ℳ, we denote by [eqn] the orthogonal projection of Y′ on [eqn], that is a least-square estimator of s′, and [eqn] the projection of s′.

We extract from Baraud et al. (2009) the following theorem, that bounds the risk of the selected estimator, and provides us with a non-asymptotic guarantee. It relies on a penalty depending on the EDkhi function, as defined below:

Definition 4.1 (Baraud et al. (2009), Section 4, definitions 2 and 3). Let D, N be two positive integers, and X_D_, X_N_ be two independent χ^2^ random variables with degrees of freedom D and N respectively. For x ≥ 0, define [eqn] And define EDkhi[D,N, q] as the unique solution of the equation Dkhi[D,N, EDkhi[D,N, q]] = q (for 0 < q ≤ 1).

Theorem 4.1 (Baraud et al. (2009), Section 4, theorem 2 and corollary 1). In the setting defined above, let Dη be the dimension of [eqn], and assume that Nη = n − Dη ≤ 2 for all η ∈ ℳ. Let [eqn] be some family of positive numbers such that [eqn], and assume that, for A > 1, [eqn]

Take [eqn] as the minimizer of the criterion: [eqn]. Then, assuming that Nη ≥ 7 and max(Lη, Dη) ≤ Kn for any η ∈ ℳ, with K < 1, the following non-asymptotic bound holds: [eqn]

The penalty used here ensures an oracle inequality: in expectation, the risk of the selected estimator is bounded by the risk of the best possible estimator of the collection of models, up to a multiplicative constant, and a residual term that depends on the dimension of the oracle model. Note that if the collection of models is poor, such an inequality has low value. We refer to Baraud et al. (2009) for a more detailed discussion of this result.

Adaptation to the Tree-Structured Framework. We use the linear formulation described in 2.3, and assume that we are in the OUsun model (this procedure would also work for a BM with a deterministic root). Then, if V is a matrix of size n, with [eqn], we have: [eqn]

We assume that a is fixed, so that the design matrix TW(α) and the structure matrix V are known and fixed. A model is defined here by the position of the shifts on the branches of the tree, i. e. by the non-zero components of Δ (with the constraint that the first component, the root, is always included in the model). We denote by [eqn] the set of allowed (parsimonious) allocations of shifts on branches (see Section 3.2), p being the maximum allowed dimension of a model. From proposition 3.6, for η ∈ ℳ, the columns vectors Tη are linearly independent, and the model Sη = Span(Ti, i ∈ η) is a linear sub-space of ℝ^n^ of dimension D_η_ = |η| = Kη + 1, Kη being the number of shifts in model η. Note that as W(α) is diagonal invertible, it does not affect the definition of the linear subspaces. The set of models is then [eqn].

We define the Mahalanobis norm associated to V^−1^ by: ‖R‖v^−1^ = R^T^V ^−1^R, ∀R ∈ ℝ^n^. The projection on S_n_ according to the metric defined by V^−1^ is then: [eqn]

For a given number of shifts K, we define the best model with K shifts as the one maximizing the likelihood, or, equivalently, minimizing the least-square criterion for models with K shifts: [eqn]

The idea is then to slice the collection of models by the number of shifts K they employ. Thanks to the EM algorithm above, we are able to select the best model in such a set. The problem is then to select a reasonable number of shifts. To compensate the increase in the likelihood due to over-fitting, using the model selection procedure described above, we select K using the following penalized criterion: [eqn]

As noted in Baraud et al. (2009), the previous criterion can equivalently be re-written in term of likelihood, as: [eqn] with [eqn]. As we use maximum-likelihood estimators, we chose this formulation for the implementation. The following proposition then holds:

Proposition 4.1 (Form of the Penalty and guaranties (α known)). Let [eqn], with [eqn], the maximum dimension of a model, with k < 1, and: [eqn]

Let A > 1 and assume that [eqn]

Suppose that [eqn] is a minimizer of 4.1 or 4.2 with this penalty. Then: [eqn] with C(A, k) a constant depending on A and k only.

The proof of this proposition can be found in Appendix E. It relies on theorem 4.1, adapting it to our tree-structured observations.

Remark 4.1. With this oracle inequality, we can see that we are missing the oracle by a log(n) term. This term is known to be unavoidable, see Baraud et al. (2009) for further explanations.

Remark 4.2. Note that the chosen penalty may depend on the topology of the tree through the term [eqn] (see Section 3.2).

Remark 4.3. The penalty involves a constant A > 1, that needs to be chosen by the user. Following Baraud et al. (2009) who tested a series of values, we fixed this constant to A = 1.1.

Simulations Studies

5

Simulations Scheme

5.1

We tested our algorithm on data simulated according to an OUsun, with varying parameters. The simulation scheme is inspired by the work of Uyeda and Harmon (2014). We first generated three distinct trees with, respectively, 64,128 and 256 tips, using a pure birth process with birth rate λ = 0.1. The tree heights were scaled to one, and their topology and branch lengths were fixed for the rest of the simulations. We then used a star-like simulation study scheme, fixing a base scenario, and exploring the space of parameters one direction at the time. The base scenario was taken to be relatively “easy”, with β1 = 0 (this parameter was fixed for the rest of the simulations), α_b_ = 3 (i.e t1/2,b = 23%), [eqn] and K_b_ = 5. The parameters then varied in the following ranges: the phylogenetic half life t_1/2_ = ln(2)/α took 11 values in [0.01, 10]; the root variance [eqn] took 9 values in [0.05, 25]; the number of shifts K took 9 values in [0, 16] (see Figures 6–7 for the exact values taken). The problem was all the more difficult that γ^2^, t1/2 or K were large.

For each simulation, the K shifts were generated in the following way. First, their values were drawn according to a mixture of two Gaussian distributions, 𝒩(4,1) and 𝒩(−4,1), in equal proportions. The mixture was chosen to avoid too many shifts of small amplitude. Then, their positions were chosen to be balanced: we first divided the tree in K segments of equal heights, and then randomly drew in each segment an edge where to place a shift. We only kept parsimonious allocations.

Each of these configurations was repeated 200 times, leading to 16200 simulated data sets. An instance of a tree with the generated data is plotted in Figure 5.

Left: Simulated configuration (with t1/2 = 0.75, γ2 = 0.5 and K = 5). The shifts positions and values are marked on the tree. The value of the character generated (positive or negative) is represented on the right. The colors of the branches correspond to the true regimes, black being the ancestral state. Right: One of the three equivalent allocations of shifts for the model inferred from the data, with corresponding vector of mean tip values. Shifts not recovered are located on pendant edges, and have low influence on the data. The two other equivalent allocations can be easily deducted from this one.

Inference Procedures

5.2

For each generated dataset, we ran our EM procedure with fixed values of [eqn], n being the number of tips of the tree. Remark that for n = 64, [eqn], and we have no hope of detecting true values of K above 8 (see remark 3.1 for an explanation of the bound in [eqn]). The number of shifts K_s_ was chosen thanks to our penalized criterion, and we kept inferences corresponding to both Ks and the true number K_t_.

We ran two sets of estimations for a either known or estimated. The computations took respectively 66 and 570 (cumulated) days of CPU time. This amounts to a mean computational time of around 6 minutes (367 seconds) for one estimation when a is fixed, and 52 minutes (3137 seconds) when a is estimated, with large differences between easy and difficult scenarios.

Scores Used

5.3

The convergence of the EM algorithm was assessed through the comparison of the likelihood of the true and estimated parameters, and the comparison of mean number of EM steps needed when a is fixed or estimated. The quality of the estimates of β1, t1/2 and γ^2^ was assessed using the coefficient of variation. The model selection procedure was evaluated by comparing the true number of shifts with the estimated one, which should be lower. We do not expect to find the exact number as some shifts, which are too small or too close of the tips, cannot be detected. To evaluate the quality of the clustering of the tips, the only quantity we can observe, we used the Adjusted Rand Index (ARI, Hubert and Arabie, 1985) between the true clustering of the tips, and the one induced by the estimated shifts. The ARI is proportional to the number of concordant pairs in two clusterings and has maximum value of 1 (for identical clusterings) and expected value of 0 (for random clusterings). Note that this score is conservative as shifts of small intensity, which are left aside by our model selection procedure, produce “artificial” groups that cannot be reconstructed.

Results

5.4

The selection strength is notoriously difficult to estimate, with large ranges of values giving similar behaviors (see Thomas et al., 2014). We hence first analyse the impact of estimating a on our estimations, showing that the main behavior of the algorithm stays the same. Then, we study the shifts reconstruction procedure.

Convergence and Likelihood. For a known, all estimations converged in less than 49 iterations, with a median number of 13 iterations. For a estimated, the number of iterations increased greatly, with a median of 69, and a fraction of estimations (around 3.3%) that reached the maximum allowed number (fixed at 1000 iterations) without converging. Unsurprisingly, the more difficult the problem, the more iterations were needed. The log-likelihoods of the estimated parameters are close to the true ones, even when a is estimated (see supplementary Figure 13 in Appendix F, first row).

Estimation of continuous parameters. Figure 6 (first row) shows that we tend to slightly over-estimate a in general. The estimation is particularly bad for large values of a (with a high variance on the result, see first box of the row), and low values of a. In this regime, the model is “over-confident”, as it finds a higher selection strength than the real one and therefore a smaller variance (second row of Figure 6). For smaller and bigger trees, the estimators behave in the same way, but with degraded or improved values, as expected. We also note that taking the true number of shifts instead of the estimated one slightly degrades our estimation of these parameters (see supplementary Figure 13 in Appendix F). The estimation of β1 is not affected by the knowledge of α or K (see Figure 7, first row), and only has an increased variance for more difficult configurations. In the remainder, we only show results obtained for estimated a as estimating α does not impact ARI, [eqn] and [eqn] (see supplementary Figure 14 in Appendix F). the algorithm with the true number of shifts, γ^2^ is highly under-estimated. We will see afterwards that this is due to “undetectable” shifts.

Estimation of the number of shifts. The way shifts were drawn ensures that they are not too small in average, and that they are located all along the tree. Still, some shifts have a very small influence on the data, and are hence hard to detect (see Figure 5). The selection model procedure almost always under-estimates the number of shifts, except in very favorable cases (Figure 7, second row). This behavior is nonetheless expected, as allowing more shifts does not guarantee that the right shifts will be found (see supplementary Figures 11 and 12 in Appendix F).

Box plots over the 200 repetitions of each set of parameters, for the phylogenetic halflife (top) and root variance (bottom) with K estimated, and α fixed to its true value (blue) or estimated (red), on a tree with 128 taxa. For better legibility, the y-axis of these two rows were re-scaled, omitting some outliers (respectively, for t1/2 and γ2, 0.82% and 0.46% of points are omitted). The whisker of the first box for t1/2 goes up to 7.5.

Same for β1(top), the number of shifts (middle) and ARI (bottom), with a estimated, and K fixed to its true value (blue) or estimated (red). As previously, the y-axis of the first row (β1) was re-scaled, omitting some outliers (1.38% of points are omitted).

Clustering of the tips. The ARI tends to be degraded for small values of a or high variance, but remains positive (Figure 7, third row). When only one shift occurs, the ARI is very unstable, but for any other value of K, it stays quite high. Finally, knowing the number of shifts does not improve the ARI.

Equivalent Solutions. When α and K are both estimated, only 5.8% of the configurations have 2 or more equivalent solutions. One inferred configuration with three equivalent solutions is presented Figure 5.

Comparison with bayou. As mentioned above, our simulation scheme, although not completely equivalent to the scheme used in Uyeda and Harmon (2014), is very similar, so that we can compare our results with theirs. The main differences lies in the facts that we took a grid on γ^2^ = σ^2^/2α instead of σ^2^, and that we took shifts with higher intensities, making the detection of shifts easier. We can see that we get the same qualitative behaviors for our estimators, with the selection strength α over or under estimated, respectively, in small or large values regions. The main difference lies in the estimation of the number of shifts. Maybe because of the priors they used (K ~ Conditional Poisson(λ= 9,K_max_ = n/2)), they tend to estimate similar numbers of shifts (centered on 9) for any set of parameters. In particular, while our method seems to be quite good at detecting situations where there are no shifts at all, theirs seems unable to catch these kind of configurations, despite the fact that their shifts have low intensity, leading to a possible over-fitting of the data.

Overall, the behavior of the algorithm is quite satisfying. Our model selection procedure avoids over-fitting, while recovering the correct clustering structure of the data. It furthermore allows for a reasonable estimation of the continuous parameters, except for a which is notoriously difficult to estimate.

Case Study: Chelonian Carapace Lenxsgth Evolution

6

Phylogenetic tree of the Chelonians. Log-transformed traits values are represented on the right. Branch colors represent the habitats. The shifts found by our EM algorithm are shown as circles, with a color indicating the value of the shift, from blue (negative) to red (positive). Boxes highlight the groups induced by the shifts. The x-scale is in million years.

Description of The Dataset

6.1

Extant species of the order Testudines, or Chelonii, are turtles and tortoises, living all across the globe, and exhibiting a wide variation in body size, from the small desert speckled tortoise (Homopus signatus, 10 cm), to the large marine leatherback sea turtle (Dermochelys coriacea, 244 cm). In order to test the hypothesis of island and marine gigantism, that could explain the extreme variations observed, Jaffe et al. (2011) compiled a dataset containing a measure of the carapace length for 226 species, along with a phylogenetic tree of theses species, spanning 210 million years (my) (see Figure 8). They assigned each species to one of four habitats: mainland-terrestrial, freshwater, marine and island-terrestrial. Then, testing several fixed regimes allocations on the branches of the tree using the method described in Butler and King (2004), they found the best support in favor of a “OU2” model that assigned one regime to each habitat. Following Uyeda and Harmon (2014), we will refer to this model as the “OU_habitat_” model. Note that this model is ambiguously defined, as it requires to assign a habitat to each ancestral species. Using proposition 3.3, we found that there were 48 equivalent parsimonious ways of doing so that respect the habitats observed at the tips of the tree. One of these habitat reconstruction is presented Figure 8.

Method

6.2

We used the version of the dataset embedded in the package geiger (Harmon et al., 2008), that contains a phylogenetic tree and a vector of log-carapaces lengths. The corresponding habitats are reported in the Appendix of Jaffe et al. (2011).

We ran our algorithm with a number of shifts going from 0 to 20. Rather than estimating α directly within the EM as we did for the simulations, we took α varying on a grid, taking 6 values regularly spaced between 0.01 and 0.1, but fixed for each estimation. We found that this approach, although computationally more intensive, gave better results. These 6 × 20 = 120 estimations took around 2 hours of CPU time. For each number of shifts, going again from 0 to 20, we kept the solution with the maximal likelihood, and we applied the model selection criterion to them. This method gave a solution with 5 shifts, and a selection strength of 0.06 (i.e. 5.5 % of the total height of the tree). Using a finer grid for α gives highly similar results, allocating shifts to the same edges. These last estimations are given below.

Results

6.3

Our method selected a solution with 5 shifts, a rather strong selection strength (t1/2 = 5.4% of the tree height), and a rather low root variance (γ^2^ = 0.22, see table 1, first column). Two of those shifts are closely related to the habitats defined in Jaffe et al. (2011) (see Figure 8). The ancestral optimal value, that applies here to two clades of freshwater turtles, is estimated to be around 38 cm. A small decrease in size for a large number of mainland and freshwater turtles is found (optimal value 24 cm). Marine turtles (super-family Chelonioidea) are found to have an increased carapace length (with an optimal value of 130 cm), as well as a clade containing soft-shell tortes (family Trionychidae, optimal size 110 cm), and a clade containing almost all island tortoises, including several sub-species of Galápagos tortoises (Geochelone nigra). Only the Ryukyu black-breasted leaf turtle (Geoemyda japonica), endemic to the Ryukyu Islands in Japan, and distant on the phylogenetic tree, is not included in this group. Note that the group also contains some mainland tortoises of the genus Geochelone, that are closely related to Galaapagos tortoises. This is typical of our method: it constructs groups that are both phenotypically and phylogenetically coherent. Finally, one species is found to have its own group, the black-knobbed map turtle (Graptemys nigrinoda), with a very low optimal value of 1.4 · 10^−20^ cm, for a measured trait of 15 cm. The fact that the shift has a very high negative value (−49 in log scale) is just an artifact due to the actualization factor on a very small branch (0.18 my, for an inferred phylogenetic half-life of 11 my). This is a rather unexpected choice of shift location. When considering the linear model as transformed by the cholesky matrix of the variance to get independent errors (as in the proof of proposition 4.1), we find a leverage of 0. 94, indicating that this species trait behaves in the transformed space as an outsider.

Comparison with other Methods

6.4

In order to compare our results to previously published ones, we reproduced some of the analysis already conducted on this dataset. We hence ran the methods described in Jaffe et al. (2011) (using the R package OUwie, with fixed positions for the shifts), Uyeda and Harmon (2014) (implemented in package bayou), Ingram and Mahler (2013) (package SURFACE), and Ho and Ané (2014b) (function OUshifts in package phylolm). See Section F.3 in Appendix F for more details on these methods and the parameters we used. hours of CPU time.

The shifts allocated on the tree by methods bayou, SURFACE and OUshifts are presented on Supplementary Figure 15 (Appendix F). We can see that 3 among the most strongly supported shifts in the posterior distribution given by bayou, as well as some among the oldest shifts found by SURFACE and OUshifts are similar to the ones found by our method. The bayou method finds equal support for many shifts, all over the tree, and the median of the posterior distribution is 17 shifts, which is pretty close to the mode of the prior put on the number of shifts (15). The SURFACE and OUshifts methods select respectively 33 and 8 shifts, including many on pendant edges, that are not easily interpretable. The backward step of SURFACE allowed to merge the regimes found for marine turtles and soft-shell tortoises that our method found to have very similar optimal values. The results of the five methods are summarized Table 1. Note that these models are not nested, due to the status assigned to the root, and to the possible convergences.

Compared to step-wise heuristics, our integrated maximum likelihood based approach allows us to have a more “global” view of the tree, and hence to select a solution that accounts better for the global structure of the trait distribution. Thanks to its rigorous model selection procedure, our model seems to report significant shifts only, that are more easily interpretable than the solutions found by other methods, and that do not rely on any chosen prior.

Table 1:: Summary of the results obtained with several methods for the Chelonian Dataset. For bayou, the median of the posterior distributions is given.