Additional Analytical Support for a New Method to Compute the Likelihood of Diversification Models

Molecular phylogenies have been increasingly recognized as an important source of information on species diversification. For many models of macroevolution, analytical likelihood formulas have been derived to infer macroevolutionary parameters from phylogenies. A few years ago, a general framework to numerically compute such likelihood formulas was proposed, which accommodates models that allow speciation and/or extinction rates to depend on diversity. This framework calculates the likelihood as the probability of the diversification process being consistent with the phylogeny from the root to the tips. However, while some readers found the framework presented in Etienne et al. (Proc R Soc Lond B Biol Sci 279(1732):1300–1309, 2012) convincing, others still questioned it (personal communication), despite numerical evidence that for special cases the framework yields the same (i.e., within double precision) numerical value for the likelihood as analytical formulas do that were independently derived for these special cases. Here we prove analytically that the likelihoods calculated in the new framework are correct for all special cases with known analytical likelihood formula. Our results thus add substantial mathematical support for the overall coherence of the general framework.


Introduction
One of the major challenges in the field of macroevolution is understanding the mechanisms underlying patterns of diversity and diversification. A very fruitful approach has been to model macroevolution as a birth-death process which reduces the problem to the specification of macroevolutionary events (i.e., speciation and extinction). However, providing likelihood expressions for these models given empirical data on speciation and extinction events is quite challenging, for the following reason. While such a likelihood is very easy to derive when full information is available for all events, typically the data involve phylogenetic trees constructed with molecular data collected from extant species alone. Hence, no extinction events and speciation events leading to extinct species are recorded in such phylogenetic trees. For a variety of models, this problem can be overcome by considering a reconstructed process, whereby the phylogeny of extant species can be regarded as a pure-birth process with time-dependent speciation rate (Nee et al. 1994). But this approach is not generally valid.
Thus, the methods employed to derive likelihood expressions are usually applicable to a limited set of models. They do not apply to models that assume that speciation and extinction rates depend on the number of species in the system. Hence, potential feedback of diversity itself on diversification rates, due to interspecific competition or niche filling, is completely ignored. The first to incorporate such feedbacks were Rabosky and Lovette (2008), who made rates dependent on the number of species present at every given moment in time, analogously to logistic growth models used in population biology. However, their model had to assume that there is no extinction for mathematical tractability, which stands in stark contrast to the empirical data: the fossil record provides us with many examples of extinct species. Etienne et al. (2012) presented a framework to compute the likelihood of phylogenetic branching times under a diversity-dependent diversification process that explicitly accounts for the influence of species that are not in the phylogeny, because they have become extinct. We note that diversity dependence as implemented in the approach of Etienne et al. (2012) does not need to start at the crown of a branching process: It can already act earlier. This feature has already been used in applications to island biogeography (Valente et al. 2015). Some of our colleagues have doubts that this framework contains a formal argument that the solution of the set of ordinary differential equations that together constitute the framework gives the likelihood of the model for a given phylogenetic tree. Instead, only numerical evidence for a small set of parameter combinations has been provided that the method yields, in the appropriate limit, the known likelihood for the standard diversity-independent (i.e., using constant rates) birth-death model. This likelihood was first provided by Nee et al. (1994), using a breaking-the-tree approach. Later, Lambert and Stadler (2013) used coalescent point process theory to provide an approach to obtain likelihood formulas for a wider set of models. These models did not include diversity dependence. For example, Lambert et al. (2015) applied their framework to the protracted birth-death model (Etienne et al. 2014), which is a generalization of the diversity-independent model where speciation is no longer an instantaneous event (Etienne and Rosindell 2011). For this model, they provided an explicit likelihood expression.
Here we provide an analytical proof that the likelihood of Etienne et al. (2012) reduces to the likelihood of Lambert et al. (2015)-and hence to that of Nee et al. (1994)-for the case of diversity-independent diversification.
The extant species belonging to a clade are often not all available for sequencing, because some species are difficult to obtain tissue from (either because they are difficult to find/catch, or because they are endangered, or because they have recently become extinct due to anthropogenic rather than natural causes) or because it is difficult to extract their DNA. This means that our data consist of a phylogenetic tree of an incomplete sample of species, and thus of an incomplete set of speciation events, even for those that lead to the species that we observe today. This incomplete sampling has been described by two different random models. The first model assumes that a fixed number of extant species are not represented in the phylogeny. This model might be appropriate for well-described taxonomic groups, such as birds, where we have a good idea of the species that are evolutionarily related, but we are simply missing some data points for the reasons mentioned above. This sampling model is called n-sampling (Lambert et al. 2015). The second model assumes that extant species are represented in the phylogeny with a fixed probability ρ. This sampling scheme is called ρ-sampling (Lambert et al. 2015), but is also referred to as f -sampling (Nee et al. 1994). The framework of Etienne et al. (2012) assumes n-sampling, but in this paper, we show that it can also be extended to incorporate ρ-sampling.
In the next section, we summarize the framework of Etienne et al. (2012), and we provide the likelihood formula analytically derived by Lambert et al. (2015) for the special case of diversity-independent but time-dependent diversification with nsampling. Then we proceed by showing that the probability generating functions of these two likelihoods are identical. We end with a discussion where we point out how the framework of Etienne et al. (2012) can be extended to include ρ-sampling and how it relates to the likelihood formula of Rabosky and Lovette (2008) for the diversity-dependent birth-death model without extinction.

The Diversity-Dependent Diversification Model
Diversification models are birth-death processes in which "birth" and "death" correspond to speciation and extinction events, respectively. In the simplest case, the per-species speciation rate λ and the per-species extinction rate μ are constants. Here we consider diversification models in which the per-species speciation and extinction rates depend on the number of species n present at time t, i.e., diversity-dependent, which we denote by λ n and μ n . We also allow the speciation and extinction rates to depend on time t, i.e., λ n (t) and μ n (t), although the latter dependence is often not explicit in our notation.
We assume that the diversification process starts at time t c from a crown, i.e., from two ancestor species. Assuming that at a later time t > t c , the process has n species, the transition probabilities in the infinitesimal time interval [t, t + dt] are from n to n + 1 species with probability nλ n (t) dt from n to n − 1 species with probability nμ n (t) dt number of species n unchanged with probability 1 − nλ n (t) dt − nμ n (t) dt.
The diversification process runs until the present time t p . We denote by P n (t) the probability that the process has n species at time t. This probability satisfies the following ordinary differential equation [ODE, called master equation or forward Kolmogorov equation (Bailey 1990 where we omit in the notation the time dependence of the speciation and extinction rates.

Sampling Models
At the present time t p , a subset of the n extant species are observed and sampled. This sampling process can been modeled in two different ways (see introduction). The first model assumes that a fixed number of species is unsampled, which corresponds to the n-sampling scheme of Ref. Lambert et al. (2015). That is, the number of extant species at t p that are not sampled, a number we denote by m p , is a model parameter.
The second model assumes that each extant species at the present time is sampled with a given probability, which has been called f -sampling (Nee et al. 1994) or ρsampling (Lambert et al. 2015). In this case, the number of unsampled species m p is a random variable, and the probability with which extant species are sampled is a model parameter, which we denote by f p .

Reconstructed Tree
A realization of the diversification process from t c to t p can be represented graphically as a tree; see Fig. 1. The complete tree shows all the species that have originated in the process (Fig. 1a). However, in practice, we have only access to the reconstructed tree, i.e., the complete tree from which we remove all the species that became extinct before the present or that were not sampled (Fig. 1b). While it would be straightforward to infer information about the diversification process based on the complete tree, this task is much more challenging when only the reconstructed tree is available. This paper deals with likelihood formulas for a reconstructed phylogenetic tree. The number of tips equals the number of sampled extant species k p . We assume that also the number of unsampled extant species is known, a number we denote by m p . The information contained in a phylogenetic tree consists of a topology and a set of branching times. For a large set of diversification models, including the diversitydependent one, all trees having the same branching times but different topologies are equally probable (Lambert and Stadler 2013). Hence, rather than computing the (a) (b) (c) Fig. 1 a Full tree where missing species are plotted as red dashed lines: the ones ending in a cross become extinct before the present, whereas the ones ending with a red dot are unsampled species at the present; b Corresponding reconstructed tree in which only extant species are present. This is the type of tree we usually work with because actual phylogenetic trees are usually obtained from molecular data taken from extant species; c Lineages-through-time plot: The green line represents the number of lineages leading to extant species (k), the red line represents lineages leading to extinct or unsampled species (m), and the blue line represents the total number of lineages (n = k + m) likelihood of a specific topology, we present formulas for the likelihood of the vector of branching times. We denote the vector of branching times by t = (t 2 , t 3 , . . . , t k p −1 ), where t k is the branching time at which the phylogenetic tree changes from k to k + 1 branches. It will be convenient to set t 0 = t 1 = t c and t k p = t p .

Likelihood Conditioning
It is common practice to condition the likelihood on the survival of both ancestor lineages to the present time (Nee et al. 1994). Indeed, we would only do an analysis on trees that have actually survived to the present. To incorporate this fact, we need to divide the unconditioned likelihood by the probability for each of the ancestor species at the crown age to have sampled extant descendants. This probability would necessarily depend on the way extant species were sampled, i.e., using either the nsampling or the f -sampling model. However, for the sake of simplicity, here we apply the same conditioning as presented in the original paper (Etienne et al. 2012), where it is required that the descendants survive to the present, but not that they are sampled. In this way, the conditioning becomes independent of the choice of the sampling scheme.
3 The Framework of Etienne et al. Etienne et al. (2012) presented an approach to compute the likelihood of a phylogeny for the diversity-dependent model. It is based on a new variable, Q k m (t), which they described as "the probability that a realization of the diversification process is consistent with the phylogeny up to time t, and has n = m +k species at time t" (Ref. Etienne et al. 2012, Box 1), where k lineages are represented in the phylogenetic tree (because they are ancestral to one of the k p species extant and sampled at present) and m additional species are present but unobserved (Fig. 1c). These species might not be in the phylogenetic tree because they became extinct before the present or because they are either not discovered or not sampled yet (see introduction). From hereon, we will refer to these species denoted by m as missing species. We cannot ignore these missing species, because in a diversity-dependent speciation process, they can influence the speciation and extinction rates.
We start by describing the computation of the variable Q k p m p (t p ), which proceeds from the crown age t c to the present time t p . It is convenient to arrange the values Q k m (t), with m = 0, 1, 2, . . ., into the vector Q k (t). The initial vector Q k=2 (t c ) is transformed into the vector Q k (t) at a later time t as follows (Ref. Etienne et al. 2012, Appendix S1, Eq. (S1)): The operators A k and B k are infinite-dimensional matrices that operate along the tree, on branches, and nodes, respectively (Fig. 2). Continuing this computation until the present time t p , we get Note that Eq.
(2) generalizes Eq. (S1) of Ref. Etienne et al. (2012) to the case in which the rates are time-dependent.

Fig. 2
An example of how to build a likelihood for a tree with k p = 4 tips. We start with a vector Q 2 (t c ) at the crown age. We use A k (t k , t k−1 ) and B k (t k ) to evolve the vector across the entire tree (on branches and nodes, respectively) up to the present time t p according to . At the present time, the likelihood accounting for m p missing species will be proportional to the m p th component of the vector We specify the different terms appearing in the right-hand side of Eq. (2): -For the initial vector Q k=2 (t c ), we assume that there are no missing species at crown age, that is, Q k=2 , during which the phylogenetic tree has k branches. Etienne et al. (2012) argued that these dynamics are given by the following ODE system (Ref. Etienne et al. 2012, Box 1, Eq. (B2)): The quantity Q k p m (t p ) is proportional to the likelihood of the tree at the present time with m unsampled extant species (see Claim 3.1 for the precise statement, including the constant of proportionality). We can collect the coefficients of Q k m (t) on the right-hand side of the ODE system in a matrix V k (t). If we do so, the system can be rewritten as which has solution -The matrix B k transforms the solution of the ODE system ending at t k into the initial condition of the ODE system starting at t k . It is a diagonal matrix with components kλ k+m dt, so that The multiplication by λ k+m dt corresponds to the probability that a speciation occurs in the time interval [t k , t k + dt]. We multiply by a factor k because we do not specify which lineage speciates (recall that we compute the likelihood of a vector of branching times rather than of a specific topology). In the likelihood expressions, we will omit the differential (a choice that is widely adopted across the vast majority of this kind of models in the literature) as it is actually not essential in parameter estimation. Therefore, we will work with a likelihood density, but for simplicity, we will refer to it as a likelihood.
Claim 3.1 Consider the diversity-dependent diversification model, given by speciation rates λ n (t) and extinction rates μ n (t). The diversification process starts at crown age t c with two ancestor species and ends at the present time t p , at which a fixed number of species m p are not sampled. A phylogenetic tree is constructed for the sampled species. Then, the likelihood that the phylogenetic tree has k p tips and vector of branching times t = (t 1 , t 2 , . . . , t k p −1 ), conditional on the event that both crown lineages survive until the present, is equal to The term Q k p m p (t p ) in the numerator of this expression is obtained from Eq. (2). The term P c (t c , t p ) in the denominator, where the subscript c stands for conditioning, is equal to where Q k=2 n (t c ) = δ n,0 .
The structure of the likelihood expression (6) can be understood intuitively. It is proportional to Q k p m p (t p ), which in Etienne et al.'s interpretation is the probability that the diversification process generates the phylogenetic tree with k p tips and m p missing species at present time t p . The combinatorial factor k p +m p m p accounts for the number of ways to select m p missing species out of a total pool of k p + m p species. The factor P c (t c , t p ) is the probability that both ancestor species at crown age t c have descendant species at the present time t p . Hence, this factor applies the likelihood conditioning. Etienne et al. (2012) provided numerical evidence that Claim 3.1 is in agreement with the likelihood provided by Nee et al. (1994) under the hypothesis of diversityindependent speciation and extinction rates and no missing species at the present. However, a rigorous analytical proof, even for this specific case, has not yet been given. In this paper, we show that Claim 3.1 holds (1) for the diversity-independent (but possibly time-dependent) case and (2) for the diversity-dependent case without extinction (i.e., extinction rate μ = 0).

The Likelihood for the Diversity-Independent Case
Claim 3.1 proposes a likelihood expression for the case with a known number of unsampled species at the present, i.e., it accounts for n-sampling. For the diversityindependent case, i.e., λ n (t) = λ(t) and μ n (t) = μ(t), the likelihood is contained in a more general result established by Lambert et al. (2015). In the following proposition, we derive an explicit likelihood expression by restricting the result of Lambert et al. to the diversity-independent case. Proposition 4.1 Consider the diversity-independent diversification model, given by speciation rates λ(t) and extinction rate μ(t). The diversification process starts at crown age t c with two ancestor species, and ends at the present time t p , at which a fixed number of species m p are not sampled. A phylogenetic tree is constructed for the sampled species. Then, the likelihood that the phylogenetic tree has k p tips and vector of branching times t = (t 1 , t 2 , . . . , t k p −1 ), conditional on the event that both crown lineages survive until the present, is equal to where we used the convention t 0 = t 1 = t c . The components m j (with j = 0, 1, . . . , k p − 1) of the vectors m, in the sum on the second line, are nonnegative integers satisfying k p −1 j=0 m j = m p . The functions ξ(t, t p ) and η(t, t p ) are given by with The functions ξ(t, t p ) and η(t, t p ) are those appearing in Kendall's solution of the birth-death model [see Ref. Kendall 1948, Eqs. (10-12)], and are useful to describe the process when time-dependent rates are involved. Given the probability P n (t, t p ) of realizing a process starting with 1 species at time t and ending with n species at time t p , we have that ξ(t, t p ) = P 0 (t, t p ) and η(t, t p ) = P n +1 (t,t p ) P n (t,t p ) for any n > 0.
Proof The likelihood for n-sampling was originally provided by Ref. Lambert et al. (2015), Eq. (7), but we start from the explicit version provided in Ref. Etienne et al. (2014), Eq. (1); see corrigendum in Ref. Etienne (2017), Etienne et al. (2014), Lambert et al. (2015) specify the functions f (t, t p ) and g(t, t p ) as the solution of a system of ODEs for the case of protracted speciation, a model where speciation does not take place instantaneously but is initiated and needs time to complete. The standard diversification model is then obtained by taking the limit in which the speciation-completion rate tends to infinity. In this limit, the four-dimensional system of Etienne et al. (2014), Eq. (2), reduces to a two-dimensional system of ODEs, Note that in this paper time t runs from past to present rather than from present to past as in Etienne et al. (2014). The conditions at the present time t p are given by g(t p , t p ) = 1 and q(t p , t p ) = 0. The solution of this system of ODEs can be expressed in terms of η(t, t p ) and ξ(t, t p ), which can be checked using the derivatives of the expressions 10 and 9 Substituting the functions f (t, t p ) and g(t, t p ) into the likelihood expression (11) concludes the proof.
The functions ξ(t, t p ) and η(t, t p ) are directly related to the functions used by Nee et al. (1994). In particular, the functions they denoted by P(t, t p ) and u t correspond in our notation to 1 − ξ(t, t p ) and η(t, t p ), respectively.
This correspondence allows us to get an intuitive understanding of the likelihood expression (8). First consider the case without missing species. Setting m p = 0, we get which is identical to the breaking-the-tree likelihood of Nee et al. (1994, Eq. (20)). In the latter approach, the phylogenetic tree is broken into single branches: two for the interval [t c , t p ] and one for each interval [t i , t p ] with i = 2, 3, . . . , k p − 1. Each branch contributes a factor (1 − ξ(t i , t p ))(1 − η(t i , t p )), equal to the probability that the branch starting at t i has a single descendant species at t p . For the two branches originating at t c , the factor (1 − ξ(t i , t p )), equal to the probability of having (one or more) descendant species, drops due to the conditioning. For the other branches, there is an additional factor λ(t i ) for the speciation events.
Next consider the case with missing species. Each of the branches resulting from breaking the tree can contribute species to the pool of m p missing species. For the branch over the interval [t j , t p ], there are m j such species, each contributing a factor η(t j , t p ) to the likelihood. Indeed, (1 − ξ(t j , t p ))(1 − η(t j , t p ))(η(t j , t p )) m j is equal to probability of having exactly m j + 1 descendant species at the present time. One of these species is represented in the phylogenetic tree, justifying the combinatorial factor (m j + 1) in the second line of Eq. (8).
Finally, we recall the expressions for the functions ξ(t, t p ) and η(t, t p ) in the case of constant rates, λ(t) = λ and μ(t) = μ,

Equivalence for the Diversity-Independent Case
Likelihood formula (8) allows speciation and extinction rates to be arbitrary functions of time, λ(t) and μ(t). Here we show that, for the diversity-independent case, we find the same likelihood formula with the approach of Etienne et al. (2012). From now on, we will use the short-hand notation ∂ x for the partial derivative with respect to the generic variable x.
Theorem 5.1 Claim 3.1 holds for the diversity-independent case.
Proof The proof relies heavily on generating functions. First, we introduce the generating function for the variables Q k m (t), The set of ODEs satisfied by Q k m (t), Eq. (3), transforms into a partial differential equation (PDE) for the generating function F k (z, t), with Note that the number of branches k changes at each branching time, so that the PDE for F k (z, t) is valid only for t k−1 ≤ t ≤ t k (corresponding to the operator A k ). At branching time t k , the solution F k (z, t k ) has to be transformed to provide the initial condition for the PDE for F k+1 (z, t) at time t k (corresponding to the operator B k ). Using Eq. (5) and dropping the differential, we get The initial condition at crown age is F 2 (z, t c ) = 1 because Q k=2 m (t c ) = δ m,0 . Next, we define P n (s, t) as the probability that the birth-death process that started with one species at time s has n species at time t. The corresponding generating function is defined as, The set of ODEs satisfied by P n (s, t), Eq. (1), transforms into a PDE, Its solution was given by Kendall (1948, Eq. (9)), where ξ(s, t) and η(s, t) are given in Eqs. (9) and (10). The generating function F k (z, t) can be expressed in terms of the generating function G(z, s, t), as shown in the following lemma.

Lemma 5.1 The generating function F k (z, t) of the variables Q k m (t) is given by
To prove the lemma, let us suppose that the solution of Eq. (13) is of the form, where we used the convention t 0 = t 1 = t c and C k (t) is a constant depending on the branching times. This expression can be rewritten as: The partial derivatives of F k can now be computed, 1 Substituting into Eq. (7) and using Eq. (21), we get Finally, substituting Eqs. (22) and (23) into the likelihood formula (6) of Claim 3.1, which is identical to likelihood formula (8). This concludes the proof of the theorem. Nee et al. (1994) noted that one way to model the sampling of extant species is equivalent to a mass extinction just before the present. This sampling model corresponds to sampling each extant species with a given probability f p , which has also been called ρ-sampling (Lambert et al. 2015). We use the link with mass extinction to extend the previous formula for n-sampling to the case of ρ-sampling. First, we formulate the ρ-sampling version of Claim 3.1.

A Note on Sampling a Fraction of Extant Species
Claim 6.1 Consider the diversity-dependent diversification model, given by speciation rates λ n (t) and extinction rates μ n (t). The diversification process starts at crown age t c with two ancestor species and ends at the present time t p , at which extant species are sampled with probability f p . Then, the likelihood of a phylogenetic tree with k p tips and branching times t, conditional on the event that both crown lineages survive until the present, is equal to The term P s (t c , t, t p , f p ) in the numerator, where the subscript s stands for sampling, is equal to where Q k p m (t p ) is obtained from Eq. (2). The term P c (t c , t p ) in the denominator, where the subscript c stands for conditioning, is equal to where Q k=2 m (t p ) is again obtained from Eq. (2).
Next, we establish as a reference the likelihood formula for ρ-sampling in the diversity-independent case. Proposition 6.1 Consider the diversity-independent diversification model, given by speciation rates λ(t) and extinction rates μ(t). The diversification process starts at crown age t c with two ancestor species, and ends at the present time t p , at which extant species are sampled with probability f p . Then, the likelihood of a phylogenetic tree with k p tips and branching times t, conditional on the event that both crown lineages survive until the present, is equal to where the operator W( f p ) is an infinite-dimensional matrix with components otherwise.
Hence, the operator C( f p ), which is also an infinite-dimensional matrix, is equal to We need the row m = 0 to evaluate the likelihood, which is equal to We are then ready to evaluate likelihood formula (6) with the modified extinction rate. Setting m p = 0, we get .
Recall that the conditioning probability P c (t c , t p ) is not affected by the process of sampling extant species. We get which is identical to Eq. (28). This ends the proof.
Note that Eq. 26 is equal to f k p p F k p (1 − f p , t p ) which provides an alternative route to prove Claim 6.1 (Manceau et al. 2019).
Finally, we give the expressions for the functions ξ(t, t p ) and η(t, t p ) in the case of constant rates, λ(t) = λ and μ(t) = μ, which are identical to Eqs. (4) and (5) in the paper by Stadler (2013). Rabosky and Lovette (2008) derived the likelihood for a particular instance of the diversity-dependent diversification model, namely, when there is no extinction. This is the only case for which a diversity-dependent likelihood formula is available. Here we show that this case is dealt with correctly in the approach of Etienne et al. (2012). We start by reformulating the result of Rabosky and Lovette (2008) in our notation.

The Diversity-Dependent Case Without Extinction
Proposition 7.1 Consider the diversity-dependent model without extinction, given by speciation rates λ n (t). The diversification process starts at crown age t c with two ancestor species, and ends at the present time t p , at which all extant species are sampled. Then, the likelihood of a phylogenetic tree with k p tips and branching times t is equal to where we used the convention t 1 = t c and t k p = t p .
Proof Equation (33) follows from Eqs. (2.4) and (2.5) in Ref. Rabosky and Lovette (2008), by noting that ξ i in their notation corresponds to exp − k p j=i t j t j−1 λ j (s) ds in our notation.
Note that in the case without extinction likelihood conditioning has no effect.
where t belongs to [t k−1 , t k ]. Note that in this time interval there are exactly k species. Given the initial condition Q k 0 (t k−1 ) at t k−1 , the solution is Using the initial condition at crown age t c , Q k=2 Substituting into Eq. (6) yields the desired result.

Concluding Remarks
We have shown here that for the diversity-independent, but time-dependent birthdeath model with n-sampling, the framework of Etienne et al. (2012) yields the same likelihood derived by Lambert et al. (2015) (also presented in a more explicit form in Etienne 2017 andEtienne et al. 2014). This provides strong support for the correctness of this framework, but does not prove that it is also correct for the case of diversity dependence. We have thus far not been able to provide alternative evidence for this framework, apart from the fact that parameter estimations on simulations of this model provide reasonable, although sometimes biased, estimates (Etienne et al. 2012). We hope that our analysis here will suggest directions for a further substantiating of the framework. The approach taken by Manceau et al. (2019) may be promising, as it also provides numerical evidence for the correctness of the framework in the diversitydependent case.
Most existing macroevolutionary models rely on the hypothesis that the subcomponents of trees do not interact (and one can thus apply a breaking-the-tree approach, as in Nee et al. 1994, p. 308), therefore letting the likelihood be a factorization of terms that comes independently from the tree's edges and nodes. However, such a hypothesis is not always valid. The diversification process likely also depends on properties of other lineages than the lineage under consideration. The analytical treatment of Etienne et al. (2012) arguments presented in this work suggests a direction toward deriving the likelihood for much more complicated models with "interacting branches," with the arguably simplest case being diversity dependence, i.e., dependence only on the total number of lineages present at any time. Our work, showing analytically that Etienne et al.'s model agrees with existing formulas for likelihoods of simple diversification models, suggests that future models that aim to deal with interacting branches should consider such a structure as a reference point, in the same fashion as models dealing with "breakable" trees often refer to Nee et al. (1994) paradigm.
In this article, we have proved that the framework to compute a likelihood for diversity-dependent processes by Etienne et al. (2012) agrees with analytical results obtained for diversity-independent diversification models. This suggests that the framework is valid for more general models that take into account the effect of diversity of speciation and extinction rates while still being able to deal with unsampled species in the phylogeny, when this number is known. Our results can thus improve the understanding of the general architecture of macroevolutionary diversification models providing useful tools for the development of new models.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.