1 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 \(\rho \). This sampling scheme is called \(\rho \)-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 \(\rho \)-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 n-sampling. 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 \(\rho \)-sampling and how it relates to the likelihood formula of Rabosky and Lovette (2008) for the diversity-dependent birth–death model without extinction.

2 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 \(\lambda \) and the per-species extinction rate \(\mu \) 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 \(\lambda _{n}\) and \(\mu _{n}\). We also allow the speciation and extinction rates to depend on time t, i.e., \(\lambda _{n}(t)\) and \(\mu _{n}(t)\), although the latter dependence is often not explicit in our notation.

We assume that the diversification process starts at time \(t_\mathrm{c}\) from a crown, i.e., from two ancestor species. Assuming that at a later time \(t>t_\mathrm{c}\), the process has n species, the transition probabilities in the infinitesimal time interval \([t,t+\mathrm{d}t]\) are

The diversification process runs until the present time \(t_\mathrm{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)],

$$\begin{aligned} \frac{\mathrm{d}P_{n}(t)}{\mathrm{d}t}=\mu _{n+1}\;(n+1)P_{n+1}(t)+\lambda _{n-1}\; (n-1)P_{n-1}(t)-(\lambda _{n}+\mu _{n})\;nP_{n}(t),\nonumber \\ \end{aligned}$$
(1)

where we omit in the notation the time dependence of the speciation and extinction rates.

2.1 Sampling Models

At the present time \(t_\mathrm{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_\mathrm{p}\) that are not sampled, a number we denote by \(m_\mathrm{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 \(\rho \)-sampling (Lambert et al. 2015). In this case, the number of unsampled species \(m_\mathrm{p}\) is a random variable, and the probability with which extant species are sampled is a model parameter, which we denote by \(f_\mathrm{p}\).

Fig. 1
figure 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\))

Full size image

2.2 Reconstructed Tree

A realization of the diversification process from \(t_\mathrm{c}\) to \(t_\mathrm{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_\mathrm{p}\). We assume that also the number of unsampled extant species is known, a number we denote by \(m_\mathrm{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 diversity-dependent one, all trees having the same branching times but different topologies are equally probable (Lambert and Stadler 2013). Hence, rather than computing the 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 \({\mathbf {t}}=(t_{2},t_{3},\ldots ,t_{k_\mathrm{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_\mathrm{c}\) and \(t_{k_\mathrm{p}}=t_\mathrm{p}\).

2.3 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 n-sampling 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_{m}^{k}(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_\mathrm{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_{m_\mathrm{p}}^{k_\mathrm{p}}(t_\mathrm{p})\), which proceeds from the crown age \(t_\mathrm{c}\) to the present time \(t_\mathrm{p}\). It is convenient to arrange the values \(Q_{m}^{k}(t)\), with \(m=0,1,2,\ldots \), into the vector \({\mathbf {Q}}^{k}(t)\). The initial vector \({\mathbf {Q}}^{k=2}(t_\mathrm{c})\) is transformed into the vector \({\mathbf {Q}}^{k}(t)\) at a later time t as follows (Ref. Etienne et al. 2012, Appendix S1, Eq. (S1)):

$$\begin{aligned} {\mathbf {Q}}^{k}(t)&=\mathrm {A}_{k}(t_{k-1},t)\, \mathrm {B}_{k-1}(t_{k-1})\,\mathrm {A}_{k-1}(t_{k-2},t_{k-1})\\&\qquad \ldots \mathrm {A}_{3}(t_{2},t_{3})\,\mathrm {B}_{2}(t_{2}) \,\mathrm {A}_{2}(t_\mathrm{c},t_{2})\,{\mathbf {Q}}^{k=2}(t_\mathrm{c}), \end{aligned}$$

with \(t_{k-1}\le t\le t_{k}\). The operators \(\mathrm {A}_{k}\) and \(\mathrm {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_\mathrm{p}\), we get

$$\begin{aligned} {\mathbf {Q}}^{k_\mathrm{p}}(t_\mathrm{p})&=\mathrm {A}_{k_\mathrm{p}}(t_{k_\mathrm{p}-1},t_\mathrm{p}) \,\mathrm {B}_{k_\mathrm{p}-1}(t_{k_\mathrm{p}-1})\,\mathrm {A}_{k_\mathrm{c}-1}(t_{k_\mathrm{p}-2}, t_{k_\mathrm{p}-1})\nonumber \\&\qquad \ldots \mathrm {A}_{3}(t_{2},t_{3})\,\mathrm {B}_{2}(t_{2}) \,\mathrm {A}_{2}(t_\mathrm{c},t_{2})\,{\mathbf {Q}}^{k=2}(t_\mathrm{c}). \end{aligned}$$
(2)

Note that Eq. (2) generalizes Eq. (S1) of Ref. Etienne et al. (2012) to the case in which the rates are time-dependent.

We specify the different terms appearing in the right-hand side of Eq. (2):

  • For the initial vector \({\mathbf {Q}}^{k=2}(t_\mathrm{c})\), we assume that there are no missing species at crown age, that is, \(Q_{m}^{k=2}(t_\mathrm{c})=\delta _{m,0}\).

  • The matrix \(\mathrm {A}_{k}\) corresponds to the dynamics of \(Q_{m}^{k}(t)\) in the time interval \([t_{k-1},t_{k}]\), 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)):

    $$\begin{aligned} \frac{\mathrm{d}Q_{m}^{k}(t)}{\mathrm{d}t}&= \mu _{k+m+1}(m+1)Q_{m+1}^{k}(t)+\lambda _{k+m-1}(m-1+2k)Q_{m-1}^{k}(t)\nonumber \\&\quad -(\lambda _{m+k}+\mu _{m+k})(m+k)Q_{m}^{k}(t), \quad \forall m > 0, \nonumber \\ \frac{\mathrm{d}Q_{0}^{k}(t)}{\mathrm{d}t}&=\mu _{k+1}Q_{1}^{k}(t)-(\lambda _{k}+\mu _{k})kQ_{0}^{k}(t), \, \quad \text {if} \, \, m = 0. \end{aligned}$$
    (3)
    Fig. 2
    figure 2

    An example of how to build a likelihood for a tree with \(k_\mathrm{p}=4\) tips. We start with a vector \({\mathbf {Q}}^{2}(t_\mathrm{c})\) at the crown age. We use \(\mathrm {A}_{k}(t_{k},t_{k-1})\) and \(\mathrm {B}_{k}(t_{k})\) to evolve the vector across the entire tree (on branches and nodes, respectively) up to the present time \(t_\mathrm{p}\) according to \({\mathbf {Q}}^{4}(t_{4})=A_{4}(t_{4},t_{3})B_{3}(t_{3}) A_{3}(t_{3},t_{2})B_{2}(t_{2})A_{2}(t_{2},t_\mathrm{c}){\mathbf {Q}}^{2}(t_\mathrm{c})\). At the present time, the likelihood accounting for \(m_\mathrm{p}\) missing species will be proportional to the \(m_\mathrm{p}\)th component of the vector \(L_{4,m_\mathrm{p}}\propto Q_{m_\mathrm{p}}^{4}(t_\mathrm{p})\)

    Full size image

    The quantity \(Q_{m}^{k_\mathrm{p}}(t_\mathrm{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_{m}^{k}(t)\) on the right-hand side of the ODE system in a matrix \(\mathrm {V}_{k}(t)\). If we do so, the system can be rewritten as

    $$\begin{aligned} \frac{\mathrm{d}{\mathbf {Q}}^{k}(t)}{\mathrm{d}t}=\mathrm {V}_{k}(t)\,{\mathbf {Q}}^{k}(t), \end{aligned}$$

    which has solution

    $$\begin{aligned} {\mathbf {Q}}^{k}(t)=\exp \Bigg (\int _{t_{k-1}}^{t}\!\!\! \mathrm {V}_{k}(s)\,\mathrm{d}s\Bigg )\,{\mathbf {Q}}^{k}(t_{k-1}), \end{aligned}$$

    so that

    $$\begin{aligned} \mathrm {A}_{k}(t_{k-1},t)=\exp \Bigg ({\displaystyle \int _{t_{k-1}}^{t}\!\!\! \mathrm {V}_{k}(s)\,\mathrm{d}s}\Bigg ). \end{aligned}$$
    (4)

  • The matrix \(\mathrm {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\lambda _{k+m}\mathrm{d}t\), so that

    $$\begin{aligned} Q_{m}^{k+1}(t_{k})=\big (\mathrm {B}_{k}(t_{k})\big )_{m,m}Q_{m}^{k} (t_{k})=k\lambda _{k+m}\mathrm{d}t\,Q_{m}^{k}(t_{k}). \end{aligned}$$
    (5)

    The multiplication by \(\lambda _{k+m}\mathrm{d}t\) corresponds to the probability that a speciation occurs in the time interval \([t_{k},t_{k}+\mathrm {d}t]\). 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.

We are then ready to formulate the claim made by Etienne et al. (2012) (in particular, from their Appendix S1, see Eqs. (S2) and (S6) to obtain Eq. 6 and Eqs. (S7–11) to obtain Eq. 7).

Claim 3.1

Consider the diversity-dependent diversification model, given by speciation rates \(\lambda _{n}(t)\) and extinction rates \(\mu _{n}(t)\). The diversification process starts at crown age \(t_c\) with two ancestor species and ends at the present time \(t_\mathrm{p}\), at which a fixed number of species \(m_\mathrm{p}\) are not sampled. A phylogenetic tree is constructed for the sampled species. Then, the likelihood that the phylogenetic tree has \(k_\mathrm{p}\) tips and vector of branching times \({\mathbf {t}} = (t_1,t_2,\ldots ,t_{k_\mathrm{p}-1})\), conditional on the event that both crown lineages survive until the present, is equal to

$$\begin{aligned} L_{k_\mathrm{p},{\mathbf {t}},m_\mathrm{p}}=\frac{Q_{m_\mathrm{p}}^{k_\mathrm{p}}(t_\mathrm{p})}{\left( {\begin{array}{c}k_\mathrm{p}+m_\mathrm{p}\\ m_\mathrm{p}\end{array}}\right) P_\mathrm{c}(t_\mathrm{c},t_\mathrm{p})}. \end{aligned}$$
(6)

The term \(Q_{m_\mathrm{p}}^{k_\mathrm{p}}(t_\mathrm{p})\) in the numerator of this expression is obtained from Eq. (2). The term \(P_\mathrm{c}(t_\mathrm{c},t_\mathrm{p})\) in the denominator, where the subscript c stands for conditioning, is equal to

$$\begin{aligned} P_\mathrm{c}(t_\mathrm{c},t_\mathrm{p}) = \sum _{m=0}^{\infty } \frac{6}{(m+2)(m+3)} \, { \sum _{n=0}^{\infty } \big (\mathrm {A}_{2}(t_\mathrm{c},t_\mathrm{p})\big )_{m,n} \, Q_{n}^{k=2}(t_\mathrm{c}) ,} \end{aligned}$$
(7)

where \(Q_{n}^{k=2}(t_\mathrm{c}) = \delta _{n,0}\).

The structure of the likelihood expression (6) can be understood intuitively. It is proportional to \(Q_{m_\mathrm{p}}^{k_\mathrm{p}}(t_\mathrm{p})\), which in Etienne et al.’s interpretation is the probability that the diversification process generates the phylogenetic tree with \(k_\mathrm{p}\) tips and \(m_\mathrm{p}\) missing species at present time \(t_\mathrm{p}\). The combinatorial factor \(\left( {\begin{array}{c}k_\mathrm{p}+m_\mathrm{p}\\ m_\mathrm{p}\end{array}}\right) \) accounts for the number of ways to select \(m_\mathrm{p}\) missing species out of a total pool of \(k_\mathrm{p}+m_\mathrm{p}\) species. The factor \(P_\mathrm{c}(t_\mathrm{c},t_\mathrm{p})\) is the probability that both ancestor species at crown age \(t_\mathrm{c}\) have descendant species at the present time \(t_\mathrm{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 diversity-independent 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 \(\mu =0\)).

4 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 diversity-independent case, i.e., \(\lambda _{n}(t)=\lambda (t)\) and \(\mu _{n}(t)=\mu (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 \(\lambda (t)\) and extinction rate \(\mu (t)\). The diversification process starts at crown age \(t_c\) with two ancestor species, and ends at the present time \(t_\mathrm{p}\), at which a fixed number of species \(m_\mathrm{p}\) are not sampled. A phylogenetic tree is constructed for the sampled species. Then, the likelihood that the phylogenetic tree has \(k_\mathrm{p}\) tips and vector of branching times \({\mathbf {t}} = (t_1,t_2,\ldots ,t_{k_\mathrm{p}-1})\), conditional on the event that both crown lineages survive until the present, is equal to

$$\begin{aligned} L_{k_{\mathrm{p}},{\mathbf {t}},m_{\mathrm{p}}}^{{\text {(div-indep)}}}&=\frac{(k_\mathrm{p}-1)!}{\left( {\begin{array}{c}k_\mathrm{p}+m_\mathrm{p}\\ m_\mathrm{p}\end{array}}\right) } \,(1-\eta (t_\mathrm{c},t_\mathrm{p}))^{2}\prod _{i=2}^{k_\mathrm{p}-1} \lambda (t_{i})(1-\xi (t_{i},t_\mathrm{p}))(1-\eta (t_{i},t_{\mathrm{p}}))\nonumber \\&\quad \times \sum _{{\mathbf {m}}\mid m_{\mathrm{p}}}\prod _{j=0}^{k_{\mathrm{p}}-1}(m_{j}+1)(\eta (t_{j},t_{\mathrm{p}}))^{m_{j}} , \end{aligned}$$
(8)

where we used the convention \(t_{0}=t_{1}=t_\mathrm{c}\). The components \(m_{j}\) (with \(j=0,1,\dots ,k_\mathrm{p}-1\)) of the vectors \(\mathbf {m}\), in the sum on the second line, are nonnegative integers satisfying \(\sum _{j=0}^{k_\mathrm{p}-1}m_{j}=m_\mathrm{p}\). The functions \(\xi (t,t_\mathrm{p})\) and \(\eta (t,t_\mathrm{p})\) are given by

$$\begin{aligned} \xi (t,t_\mathrm{p})&=1-\frac{1}{\alpha (t,t_\mathrm{p})+\int _{t}^{t_\mathrm{p}} \lambda (s)\,\alpha (t,s)\,\mathrm{d}s}\nonumber \\&=1-\frac{1}{1+\int _{t}^{t_\mathrm{p}}\mu (s)\, \alpha (t,s)\,\mathrm{d}s} \end{aligned}$$
(9)
$$\begin{aligned} \eta (t,t_\mathrm{p})&=1-\frac{\alpha (t,t_\mathrm{p})}{\alpha (t,t_\mathrm{p}) +\int _{t}^{t_\mathrm{p}}\lambda (s)\,\alpha (t,s)\,\mathrm{d}s}\nonumber \\&=1-\frac{\alpha (t,t_\mathrm{p})}{1+\int _{t}^{t_\mathrm{p}}\mu (s) \,\alpha (t,s)\,\mathrm{d}s}, \end{aligned}$$
(10)

with

$$\begin{aligned} \alpha (t,s)=\exp \Big (\int _{t}^{s}(\mu (s')-\lambda (s'))\,\mathrm{d}s'\Big ). \end{aligned}$$

The functions \(\xi (t,t_\mathrm{p})\) and \(\eta (t,t_\mathrm{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_\mathrm{p})\) of realizing a process starting with 1 species at time t and ending with n species at time \(t_\mathrm{p}\), we have that \(\xi (t,t_\mathrm{p}) = P_{0}(t,t_\mathrm{p})\) and \(\eta (t,t_\mathrm{p}) = \frac{P_{n^\star +1}(t,t_\mathrm{p})}{P_{n^\star }(t,t_\mathrm{p})}\) for any \(n^\star > 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),

$$\begin{aligned} L_{k_\mathrm{p},{\mathbf {t}},m_\mathrm{p}}^{\text {(div-indep)}}&=\frac{(k_\mathrm{p}-1)!}{\left( {\begin{array}{c}k_\mathrm{p}+m_\mathrm{p}\\ m_\mathrm{p}\end{array}}\right) } \,(g(t_\mathrm{c},t_\mathrm{p}))^{2}\prod _{i=2}^{k_\mathrm{p}-1}f(t_{i},t_\mathrm{p})\sum _{\mathbf {m}\mid m_\mathrm{p}}\prod _{j=0}^{k_\mathrm{p}-1}(m_{j}+1) (1-g(t_{j},t_\mathrm{p}))^{m_{j}}. \end{aligned}$$
(11)

Etienne et al. (2014), Lambert et al. (2015) specify the functions \(f(t,t_\mathrm{p})\) and \(g(t,t_\mathrm{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,

$$\begin{aligned} f(t,t_\mathrm{p})=\frac{\mathrm{d}g(t,t_\mathrm{p})}{\mathrm{d}t}&=\lambda (t)\,(1-q(t,t_\mathrm{p}))\,g(t,t_\mathrm{p})\\ \frac{\mathrm{d}q(t,t_\mathrm{p})}{\mathrm{d}t}&=-\mu (t)+(\mu (t)+\lambda (t))\,q(t,t_\mathrm{p})-\lambda (t)\,q^{2}(t,t_\mathrm{p}). \end{aligned}$$

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_\mathrm{p}\) are given by \(g(t_\mathrm{p},t_\mathrm{p})=1\) and \(q(t_\mathrm{p},t_\mathrm{p})=0\).

The solution of this system of ODEs can be expressed in terms of \(\eta (t,t_\mathrm{p})\) and \(\xi (t,t_\mathrm{p})\),

$$\begin{aligned} f(t,t_\mathrm{p})&=\lambda (t)\,(1-\xi (t,t_\mathrm{p}))\,(1-\eta (t,t_\mathrm{p}))\\ g(t,t_\mathrm{p})&=1-\eta (t,t_\mathrm{p})\\ q(t,t_\mathrm{p})&=\xi (t,t_\mathrm{p}), \end{aligned}$$

which can be checked using the derivatives of the expressions 10 and 9

$$\begin{aligned} \frac{\partial \eta (t,t_\mathrm{p})}{\partial t}&=-\lambda (t)\,(1-\xi (t,t_\mathrm{p}))\,(1-\eta (t,t_\mathrm{p}))\\ \frac{\partial \xi (t,t_\mathrm{p})}{\partial t}&=-(\mu (t)-\lambda (t)\,\xi (t,t_\mathrm{p}))\,(1-\xi (t,t_\mathrm{p})). \end{aligned}$$

Substituting the functions \(f(t,t_\mathrm{p})\) and \(g(t,t_\mathrm{p})\) into the likelihood expression (11) concludes the proof. \(\square \)

The functions \(\xi (t,t_\mathrm{p})\) and \(\eta (t,t_\mathrm{p})\) are directly related to the functions used by Nee et al. (1994). In particular, the functions they denoted by \(P(t,t_\mathrm{p})\) and \(u_{t}\) correspond in our notation to \(1-\xi (t,t_\mathrm{p})\) and \(\eta (t,t_\mathrm{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_\mathrm{p}=0\), we get

$$\begin{aligned} L_{k_\mathrm{p},{\mathbf {t}},0}^{\text {(div-indep)}} =(k_\mathrm{p}-1)!\,(1-\eta (t_\mathrm{c},t_\mathrm{p}))^{2}\prod _{i=2}^{k_\mathrm{p}-1} \lambda (t_{i})(1-\xi (t_{i},t_\mathrm{p}))(1-\eta (t_{i},t_\mathrm{p})), \end{aligned}$$

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_\mathrm{c},t_\mathrm{p}]\) and one for each interval \([t_{i},t_\mathrm{p}]\) with \(i=2,3,\ldots ,k_\mathrm{p}-1\). Each branch contributes a factor \((1-\xi (t_{i},t_\mathrm{p}))(1-\eta (t_{i},t_\mathrm{p}))\), equal to the probability that the branch starting at \(t_{i}\) has a single descendant species at \(t_\mathrm{p}\). For the two branches originating at \(t_\mathrm{c}\), the factor \((1-\xi (t_{i},t_\mathrm{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 \(\lambda (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_\mathrm{p}\) missing species. For the branch over the interval \([t_{j},t_\mathrm{p}]\), there are \(m_{j}\) such species, each contributing a factor \(\eta (t_{j},t_\mathrm{p})\) to the likelihood. Indeed, \((1-\xi (t_{j},t_\mathrm{p}))(1-\eta (t_{j},t_\mathrm{p}))(\eta (t_{j},t_\mathrm{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 \(\xi (t,t_\mathrm{p})\) and \(\eta (t,t_\mathrm{p})\) in the case of constant rates, \(\lambda (t)=\lambda \) and \(\mu (t)=\mu \),

$$\begin{aligned} \xi (t,t_\mathrm{p})&=\frac{\mu \big (1-\mathrm {e}^{-(\lambda -\mu )(t_\mathrm{p}-t)}\big )}{\lambda -\mu \,\mathrm {e}^{-(\lambda -\mu )(t_\mathrm{p}-t)}}\\ \eta (t,t_\mathrm{p})&=\frac{\lambda \big (1-\mathrm {e}^{-(\lambda -\mu )(t_\mathrm{p}-t)}\big )}{\lambda -\mu \,\mathrm {e}^{-(\lambda -\mu )(t_\mathrm{p}-t)}}. \end{aligned}$$

5 Equivalence for the Diversity-Independent Case

Likelihood formula (8) allows speciation and extinction rates to be arbitrary functions of time, \(\lambda (t)\) and \(\mu (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 \(\partial _{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_{m}^{k}(t)\),

$$\begin{aligned} F_{{k}}(z,t)=\sum _{m=0}^{\infty }z^{m}Q_{m}^{k}(t). \end{aligned}$$
(12)

The set of ODEs satisfied by \(Q_{m}^{k}(t)\), Eq. (3), transforms into a partial differential equation (PDE) for the generating function \(F_{k}(z,t)\),

$$\begin{aligned} \partial _{t} F_{k}(z,t)&= (\mu (t)-z\,\lambda (t))(1-z) \partial _{z} F_{k}(z,t) + k(2z\,\lambda (t)-\lambda (t)-\mu (t)) F_{k}(z,t) \nonumber \\&= c(z,t) \partial _{z} F_{k}(z,t) + k \partial _{z} c(z,t) F_{k}(z,t) , \end{aligned}$$
(13)

with

$$\begin{aligned} c(z,t)=(\mu (t)-z\lambda (t))(1-z). \end{aligned}$$

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}\le t\le 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

$$\begin{aligned} F_{k+1}(z,t_{k})=k\lambda (t_{k})\,F_{k}(z,t_{k}). \end{aligned}$$
(14)

The initial condition at crown age is \(F_{2}(z,t_\mathrm{c})=1\) because \(Q_{m}^{k=2}(t_\mathrm{c})=\delta _{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,

$$\begin{aligned} G(z,s,t)=\sum _{n=0}^{\infty }z^{n}P_{n}(s,t). \end{aligned}$$
(15)

The set of ODEs satisfied by \(P_{n}(s,t)\), Eq. (1), transforms into a PDE,

$$\begin{aligned} \partial _{t} G(z,s,t) = c(z,t) \partial _{z} G(z,s,t). \end{aligned}$$
(16)

Its solution was given by Kendall (1948, Eq. (9)),

$$\begin{aligned} G(z,s,t)=\frac{\xi (s,t)+(1-\xi (s,t) -\eta (s,t))z}{1-z\eta (s,t)}, \end{aligned}$$
(17)

where \(\xi (s,t)\) and \(\eta (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_{m}^{k}(t)\) is given by

$$\begin{aligned} F_{k}(z,t)=H^{2}(z,t_\mathrm{c},t)\prod _{j=2}^{k-1}j\lambda _{j} (t_{j})H(z,t_{j},t) \end{aligned}$$
(18)

with

$$\begin{aligned} H(z,s,t) = \partial _{z} G(z,s,t) = \frac{(1-\xi (s,t))\,(1-\eta (s,t))}{(1-z\,\eta (s,t))^{2}}. \end{aligned}$$
(19)

To prove the lemma, let us suppose that the solution of Eq. (13) is of the form,

$$\begin{aligned} F_{k}(z,t)=C_{k}({\mathbf {t}})\prod _{j=0}^{k-1} \partial _{z}G(z,t_{j},t) \end{aligned}$$
(20)

where we used the convention \(t_{0}=t_{1}=t_\mathrm{c}\) and \(C_{k}({\mathbf {t}})\) is a constant depending on the branching times. This expression can be rewritten as:

$$\begin{aligned} F_{k}(z,t)=C_{k}({\mathbf {t}})\,\frac{1}{k}\sum _{i=0}^{k-1} \partial _{z}G(z,t_{i},t)\prod _{j\ne i,j=0}^{k-1}\partial _{z}G(z,t_{j},t). \end{aligned}$$

The partial derivatives of \(F_{k}\) can now be computed,

$$\begin{aligned} \partial _{z}F_{k}= & {} C_{k}({\mathbf {t}})\sum _{i=0}^{k-1}\partial _{z}^{2} G(z,t_{i},t)\prod _{j\ne i,j=0}^{k-1}\partial _{z}G(z,t_{j},t)\\ \partial _{t}F_{k}= & {} C_{k}({\mathbf {t}})\sum _{i=0}^{k-1}\partial _{t} \partial _{z}G(z,t_{i},t)\prod _{j\ne i,j=0}^{k-1}\partial _{z}G(z,t_{j},t). \end{aligned}$$

We substitute these expressions into the PDE, Eq. (13),

$$\begin{aligned}&\sum _{i=0}^{k-1}\partial _{t}\partial _{z}G(z,t_{i},t) \prod _{j\ne i,j=0}^{k-1}\partial _{z}G(z,t_{j},t)\\&\quad =c(z,t)\sum _{i=0}^{k-1}\partial _{z}^{2}G(z,t_{i},t) \prod _{j\ne i,j=0}^{k-1}\partial _{z}G(z,t_{j},t)\\&\qquad +k\,\partial _{z}c(z,t)\,\frac{1}{k}\,\sum _{i=0}^{k-1} \partial _{z}G(z,t_{i},t)\prod _{j\ne i,j=0}^{k-1}\partial _{z}G(z,t_{j},t). \end{aligned}$$

This equation is satisfied if the following equation is satisfied for every \(i=0,1,\ldots ,k\),

$$\begin{aligned}&\partial _{t}\partial _{z}G(z,t_{i},t)\prod _{j\ne i,j=0}^{k-1}\partial _{z}G(z,t_{j},t)\\&\quad =c(z,t)\,\partial _{z}^{2}G(z,t_{i},t)\prod _{j\ne i,j=0}^{k-1}\partial _{z}G(z,t_{j},t)\\&\qquad +\partial _{z}c(z,t)\,\partial _{z}G(z,t_{i},t)\prod _{j\ne i,j=0}^{k-1}\partial _{z}G(z,t_{j},t). \end{aligned}$$

This is the case if

$$\begin{aligned} \partial _{t}\partial _{z}G(z,t_{i},t)=c(z,t)\,\partial _{z}^{2} G(z,t_{i},t)+\partial _{z}c(z,t)\,\partial _{z}G(z,t_{i},t), \end{aligned}$$

or, equivalently, if

$$\begin{aligned} \partial _{z}\big [\partial _{t}G(z,t_{i},t)\big ] =\partial _{z}\big [c(z,t)\,\partial _{z}G(z,t_{i},t)\big ]. \end{aligned}$$

This is an identity because \(G(z,t_{i},t\)) satisfies Eq. (16).

Next, we verify that the constants \(C_{k}({\mathbf {t}})\) can be determined such that initial conditions (14) are satisfied. This is indeed the case if we take

$$\begin{aligned} C_{k}({\mathbf {t}})=\prod _{j=2}^{k-1}j\lambda (t_{j}). \end{aligned}$$

Introducing the function H(z, s, t) and using \(t_{0}=t_{1}=t_\mathrm{c}\) complete the proof of the lemma.

Next, we use Eq. (18) to derive an explicit expression for the likelihood (6) of Claim 3.1. It will be useful to have explicit expressions for derivatives of the function H(z, s, t). It follows from Eq. (19) that

$$\begin{aligned} \frac{1}{a!}\partial _{z}^{a}\big [H^{b}(z,t_{j},t)\big ]= \left( {\begin{array}{c}a+2b-1\\ a\end{array}}\right) H^{b}(z,t_{j},t)\left( \frac{\eta (t_{j},t)}{1-z\, \eta (t_{j},t)}\right) ^{a}, \end{aligned}$$
(21)

where a and b are positive integers.

To evaluate the numerator of Eq. (6), we have to extract \(Q_{m_\mathrm{p}}^{k_\mathrm{p}}(t_\mathrm{p})\) from the generating function \(F_{k_\mathrm{p}}(z,t_\mathrm{p})\). Using Leibniz’ rule,

$$\begin{aligned} Q_{m_\mathrm{p}}^{k_\mathrm{p}}(t_\mathrm{p})&=\frac{1}{m_\mathrm{p}!}\partial _{z}^{m_\mathrm{p}} \left[ F_{k_\mathrm{p}}(z,t_\mathrm{p})\right] _{z=0}\nonumber \\&=\frac{C_{k_\mathrm{p}}({\mathbf {t}})}{m_\mathrm{p}!}\;\partial _{z}^{m_\mathrm{p}} \left[ \prod _{j=0}^{k_\mathrm{p}-1}H(z,t_{j},t_\mathrm{p})\right] _{z=0}\nonumber \\&=\frac{C_{k_\mathrm{p}}({\mathbf {t}})}{m_\mathrm{p}!}\; \sum _{\mathbf {m}|m_\mathrm{p}}\left( {\begin{array}{c}m_\mathrm{p}\\ m_{0},m_{1},\ldots ,m_{k_\mathrm{p}-1}\end{array}}\right) \prod _{j=0}^{k_\mathrm{p}-1}\partial _{z}^{m_{j}}\left[ H(z,t_{j},t_\mathrm{p})\right] _{z=0}\nonumber \\&=\frac{C_{k_\mathrm{p}}({\mathbf {t}})}{m_\mathrm{p}!}\sum _{\mathbf {m}|m_\mathrm{p}} \frac{m_\mathrm{p}!}{\prod _{i}m_{i}!}\prod _{j=0}^{k_\mathrm{p}-1}(m_{j}+1)!\nonumber \\&\quad \times \left[ H(z,t_{j},t_\mathrm{p})\left( \frac{\eta (t_{j},t_\mathrm{p})}{1-z\, \eta (t_{j},t_\mathrm{p})}\right) ^{m_{j}}\right] _{z=0}\nonumber \\&=C_{k_\mathrm{p}}({\mathbf {t}})\prod _{j=0}^{k_\mathrm{p}-1}H(0,t_{j},t_\mathrm{p}) \sum _{\mathbf {m}|m_\mathrm{p}}\frac{1}{\prod _{i=0}^{k_\mathrm{p}-1}m_{i}!} \prod _{j=0}^{k_\mathrm{p}-1}(m_{j}+1)!\,\eta ^{m_{j}}(t_{j},t_\mathrm{p})\nonumber \\&=\prod _{j=2}^{k_\mathrm{p}-1}j\lambda (t_{j})\,\prod _{j=0}^{k_\mathrm{p}-1}(1- \xi (t_{j},t_\mathrm{p}))(1-\eta (t_{j},t_\mathrm{p}))\nonumber \\&\quad \times \sum _{\mathbf {m}|m_\mathrm{p}} \prod _{j=0}^{k_\mathrm{p}-1}(m_{j}+1)\,\eta ^{m_{j}}(t_{j},t_\mathrm{p})\nonumber \\&=(k_\mathrm{p}-1)!(1-\xi (t_\mathrm{c},t_\mathrm{p}))^{2}(1-\eta (t_\mathrm{c},t_\mathrm{p}))^{2}\nonumber \\&\quad \times \prod _{j=2}^{k_\mathrm{p}-1}\lambda (t_{j})\,(1-\xi (t_{j},t_\mathrm{p}))(1- \eta (t_{j},t_\mathrm{p}))\nonumber \\&\quad \times \sum _{\mathbf {m}|m_\mathrm{p}}\prod _{j=0}^{k_\mathrm{p}-1}(m_{j}+1) \,\eta ^{m_{j}}(t_{j},t_\mathrm{p}). \end{aligned}$$
(22)

To evaluate the denominator of Eq. (6), we have to extract \(Q_{m}^{k=2}(t_\mathrm{p})\) from the generating function,

$$\begin{aligned} Q_{m}^{k=2}(t_\mathrm{p})=\frac{1}{m!}\partial _{z}^{m}\big [ F_{2}(z,t_\mathrm{p})\big ]_{z=0}=\frac{1}{m!}\partial _{z}^{m}\big [ H^{2}(z,t_\mathrm{c},t_\mathrm{p})\big ]_{z=0}. \end{aligned}$$

Substituting into Eq. (7) and using Eq. (21), we get

$$\begin{aligned} P_\mathrm{c}(t_\mathrm{c},t_\mathrm{p})&=\sum _{m=0}^{\infty }\frac{6}{(m+2)(m+3)} \,\frac{1}{m!}\partial _{z}^{m}\big [H^{2}(z,t_\mathrm{c},t_\mathrm{p})\big ]_{z=0}\nonumber \\&=H^{2}(0,t_\mathrm{c},t_\mathrm{p})\sum _{m=0}^{\infty }(m+1)\,\eta ^{m}(t_\mathrm{c},t_\mathrm{p})\nonumber \\&=(1-\xi (t_\mathrm{c},t_\mathrm{p}))^{2}. \end{aligned}$$
(23)

Finally, substituting Eqs. (22) and (23) into the likelihood formula (6) of Claim 3.1,

$$\begin{aligned} L_{k_\mathrm{p},{\mathbf {t}},m_\mathrm{p}}&=\frac{(k_\mathrm{p}-1)!}{ \left( {\begin{array}{c}k_\mathrm{p}+m_\mathrm{p}\\ m_\mathrm{p}\end{array}}\right) }(1-\eta (t_{1},t_\mathrm{p}))^{2} \prod _{j=2}^{k_\mathrm{p}-1}\lambda (t_{j})(1-\xi (t_{j},t_\mathrm{p})) (1-\eta (t_{j},t_\mathrm{p}))\nonumber \\&\quad \times \sum _{\mathbf {m}|m_\mathrm{p}}\prod _{j=0}^{k_\mathrm{p}-1} (m_{j}+1)\,\eta ^{m_{j}}(t_{j},t_\mathrm{p}), \end{aligned}$$
(24)

which is identical to likelihood formula (8). This concludes the proof of the theorem. \(\square \)

6 A Note on Sampling a Fraction of Extant Species

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_\mathrm{p}\), which has also been called \(\rho \)-sampling (Lambert et al. 2015). We use the link with mass extinction to extend the previous formula for n-sampling to the case of \(\rho \)-sampling.

First, we formulate the \(\rho \)-sampling version of Claim 3.1.

Claim 6.1

Consider the diversity-dependent diversification model, given by speciation rates \(\lambda _{n}(t)\) and extinction rates \(\mu _{n}(t)\). The diversification process starts at crown age \(t_\mathrm{c}\) with two ancestor species and ends at the present time \(t_\mathrm{p}\), at which extant species are sampled with probability \(f_\mathrm{p}\). Then, the likelihood of a phylogenetic tree with \(k_\mathrm{p}\) tips and branching times \({\mathbf {t}}\), conditional on the event that both crown lineages survive until the present, is equal to

$$\begin{aligned} L_{k_\mathrm{p},{\mathbf {t}}} = \frac{ P_\mathrm{s}(t_\mathrm{c},{\mathbf {t}},t_\mathrm{p},f_\mathrm{p}) }{ P_\mathrm{c}(t_\mathrm{c},t_\mathrm{p}) }. \end{aligned}$$
(25)

The term \(P_\mathrm{s}(t_\mathrm{c},{\mathbf {t}},t_\mathrm{p},f_\mathrm{p})\) in the numerator, where the subscript s stands for sampling, is equal to

$$\begin{aligned} P_\mathrm{s}(t_\mathrm{c},{\mathbf {t}},t_\mathrm{p},f_\mathrm{p}) = \sum _{m=0}^{\infty } f_\mathrm{p}^{k_\mathrm{p}} \, (1 - f_\mathrm{p})^{m} \, Q_{m}^{k_\mathrm{p}}(t_\mathrm{p}) , \end{aligned}$$
(26)

where \(Q_{m}^{k_\mathrm{p}}(t_\mathrm{p})\) is obtained from Eq. (2). The term \(P_\mathrm{c}(t_\mathrm{c},t_\mathrm{p})\) in the denominator, where the subscript c stands for conditioning, is equal to

$$\begin{aligned} P_\mathrm{c}(t_\mathrm{c},t_\mathrm{p}) = \sum _{m=0}^{\infty } \frac{6}{(m+2)(m+3)} \, Q_{m}^{k=2}(t_\mathrm{p}) , \end{aligned}$$
(27)

where \(Q_{m}^{k=2}(t_\mathrm{p})\) is again obtained from Eq. (2).

Next, we establish as a reference the likelihood formula for \(\rho \)-sampling in the diversity-independent case.

Proposition 6.1

Consider the diversity-independent diversification model, given by speciation rates \(\lambda (t)\) and extinction rates \(\mu (t)\). The diversification process starts at crown age \(t_\mathrm{c}\) with two ancestor species, and ends at the present time \(t_\mathrm{p}\), at which extant species are sampled with probability \(f_\mathrm{p}\). Then, the likelihood of a phylogenetic tree with \(k_\mathrm{p}\) tips and branching times \({\mathbf {t}}\), conditional on the event that both crown lineages survive until the present, is equal to

$$\begin{aligned} L_{k_\mathrm{p},{\mathbf {t}}}^{\text {(div-indep)}} =(k_\mathrm{p}-1)!\,(1-{\widetilde{\eta }}(t_\mathrm{c},t_\mathrm{p}))^{2} \prod _{i=2}^{k_\mathrm{p}-1}\lambda (t_{i})(1-{\widetilde{\xi }}(t_{i},t_\mathrm{p})) (1-{\widetilde{\eta }}(t_{i},t_\mathrm{p})).\qquad \quad \end{aligned}$$
(28)

The functions \({\widetilde{\xi }}(t,t_\mathrm{p})\) and \({\widetilde{\eta }}(t,t_\mathrm{p})\) are given by

$$\begin{aligned} {\widetilde{\xi }}(t,t_\mathrm{p})&=1-\frac{f_\mathrm{p}}{\alpha (t,t_\mathrm{p}) +f_\mathrm{p}\int _{t}^{t_\mathrm{p}}\lambda (s)\,\alpha (t,s)\,\mathrm{d}s}\nonumber \\&=1-\frac{f_\mathrm{p}}{f_\mathrm{p}+(1-f_\mathrm{p})\,\alpha (t,t_\mathrm{p})+f_\mathrm{p} \int _{t}^{t_\mathrm{p}}\mu (s)\,\alpha (t,s)\,\mathrm{d}s} \end{aligned}$$
(29)
$$\begin{aligned} {\widetilde{\eta }}(t,t_\mathrm{p})&=1-\frac{\alpha (t,t_\mathrm{p})}{\alpha (t,t_\mathrm{p})+f_\mathrm{p}\int _{t}^{t_\mathrm{p}}\lambda (s)\,\alpha (t,s)\,\mathrm{d}s}\nonumber \\&=1-\frac{\alpha (t,t_\mathrm{p})}{f_\mathrm{p}+(1-f_\mathrm{p})\, \alpha (t,t_\mathrm{p})+f_\mathrm{p}\int _{t}^{t_\mathrm{p}}\mu (s)\,\alpha (t,s)\,\mathrm{d}s}, \end{aligned}$$
(30)

with

$$\begin{aligned} \alpha (t,s)=\exp \Big (\int _{t}^{s}(\mu (s')-\lambda (s'))\,\mathrm{d}s'\Big ). \end{aligned}$$

Proof

We use the equivalence between \(\rho \)-sampling and a mass extinction, see Ref. Nee et al. (1994), Eq. (31). We introduce a modified extinction rate \(\mu (t)\) containing a delta function just before the present,

$$\begin{aligned} {\widetilde{\mu }}(t)=\mu (t)-\ln f_\mathrm{p}\ \delta (t-t_\mathrm{p}). \end{aligned}$$
(31)

The likelihood formula is then obtained by setting \(m_\mathrm{p}=0\) in Eq. (8), while evaluating the functions \(\xi (t,t_\mathrm{p})\) and \(\eta (t,t_\mathrm{p})\) with the modified extinction rate \({\widetilde{\mu }}(t,t_\mathrm{p})\). This establishes Eq. (28); it remains to be proven that the modified functions \({\widetilde{\xi }}(t,t_\mathrm{p})\) and \({\widetilde{\mu }}(t,t_\mathrm{p})\) are given by Eqs. (29) and (30). This follows by noting that the modified version \({\widetilde{\alpha }}(t,t_\mathrm{p})\) of the function \(\alpha (t,t_\mathrm{p})\) appearing in Eqs. (9) and (10) satisfies

$$\begin{aligned} {\widetilde{\alpha }}(t,t_\mathrm{p})&=\exp \Big (\int _{t}^{t_\mathrm{p}} ({\widetilde{\mu }}(s)-\lambda (s))\,\mathrm{d}s\Big )\\&=\frac{1}{f_\mathrm{p}}\,\exp \Big (\int _{t}^{t_\mathrm{p}}(\mu (s) -\lambda (s))\,\mathrm{d}s\Big )=\frac{1}{f_\mathrm{p}}\,\alpha (t,t_\mathrm{p}), \end{aligned}$$

while \({\widetilde{\alpha }}(t,s)=\alpha (t,s)\) if \(s<t_\mathrm{p}\). \(\square \)

We are then ready to establish the following result.

Theorem 6.1

Claim 6.1 holds for the diversity-independent case.

Proof

We use again the equivalence between \(\rho \)-sampling and a mass extinction; see Eq. (31). Due to Theorem 5.1, likelihood formula (6) is valid for the diversity-independent case. Hence, we can derive the corresponding likelihood formula for \(\rho \)-sampling by introducing the modified extinction rate \({\widetilde{\mu }}(t)\), and setting \(m_\mathrm{p}=0\) in the likelihood formula for n-sampling.

The introduction of the modified extinction rate \({\widetilde{\mu }}(t,t_\mathrm{p})\) corresponds to applying an additional operator to the vector \({\mathbf {Q}}^{k_\mathrm{p}}(t_\mathrm{p})\) at the present time. In particular, the modified vector \({\widetilde{{\mathbf {Q}}}}^{k_\mathrm{p}}(t_\mathrm{p})\) is given by

$$\begin{aligned} {\widetilde{{\mathbf {Q}}}}^{k_\mathrm{p}}(t_\mathrm{p})=\mathrm {C}(f_\mathrm{p}) \,{\mathbf {Q}}^{k_\mathrm{p}}(t_\mathrm{p}), \end{aligned}$$
(32)

where the operator \(\mathrm {C}(f_\mathrm{p})\) corresponds to the following ODE, acting in a small time interval \([t_\mathrm{p}-\varepsilon ,t_\mathrm{p}]\) before the present,

$$\begin{aligned} \frac{d{\widetilde{Q}}_{m}^{k_\mathrm{p}}(t)}{\mathrm{d}t}&= \big (\mu -\frac{1}{\varepsilon }\ln f_\mathrm{p}\big )(m+1) {\widetilde{Q}}_{m+1}^{k_\mathrm{p}}(t)+\lambda (m-1+2k_\mathrm{p}) {\widetilde{Q}}_{m-1}^{k_\mathrm{p}}(t)\\&\quad -\big (\lambda +\big (\mu -\frac{1}{\varepsilon }\ln f_\mathrm{p}\big )\big )(m+k_\mathrm{p}){\widetilde{Q}}_{m}^{k_\mathrm{p}}(t), \end{aligned}$$

where we added a delta peak to the extinction rate, Eq. (31), in the ODE satisfied by \({\mathbf {Q}}^{k}(t)\), Eq. (3). In the limit \(\varepsilon \rightarrow 0\), the terms in \(\frac{1}{\varepsilon }\) dominate, so that

$$\begin{aligned} \frac{d{\widetilde{Q}}_{m}^{k_\mathrm{p}}(t)}{\mathrm{d}t}=-\frac{1}{\varepsilon }\ln f_\mathrm{p}\,(m+1){\widetilde{Q}}_{m+1}^{k_\mathrm{p}}(t)+\frac{1}{\varepsilon }\ln f_\mathrm{p}\,(m+k_\mathrm{p}){\widetilde{Q}}_{m}^{k_\mathrm{p}}(t). \end{aligned}$$

This can be rewritten in matrix form as

$$\begin{aligned} \frac{d{\widetilde{{\mathbf {Q}}}}^{k_\mathrm{p}}(t)}{\mathrm{d}t}=\frac{1}{\varepsilon } \mathrm {W}(f_\mathrm{p})\,{\widetilde{{\mathbf {Q}}}}^{k_\mathrm{p}}, \end{aligned}$$

where the operator \(\mathrm {W}(f_\mathrm{p})\) is an infinite-dimensional matrix with components

$$\begin{aligned} \mathrm {W}_{m,n}(f_\mathrm{p})={\left\{ \begin{array}{ll} \ln f_\mathrm{p}\,(m+k_\mathrm{p}) &{} \text {if }m=n\\ -\ln f_\mathrm{p}\,(m+1) &{} \text {if }m=n-1\\ 0 &{} \text {otherwise.} \end{array}\right. } \end{aligned}$$

Hence, the operator \(\mathrm {C}(f_\mathrm{p})\), which is also an infinite-dimensional matrix, is equal to

$$\begin{aligned} \mathrm {C}(f_\mathrm{p})=\exp \Bigg (\int _{t_\mathrm{p}-\varepsilon }^{t_\mathrm{p}} \frac{1}{\varepsilon }\mathrm {W}(f_\mathrm{p})\,\mathrm{d}s\Bigg )=\exp \big (\mathrm {W}(f_\mathrm{p})\big ). \end{aligned}$$

We need the row \(m=0\) to evaluate the likelihood, which is equal to

$$\begin{aligned} \mathrm {C}_{m=0,n}(f_\mathrm{p})=f_\mathrm{p}^{k_\mathrm{p}}\,(1-f_\mathrm{p})^{n}. \end{aligned}$$

We are then ready to evaluate likelihood formula (6) with the modified extinction rate. Setting \(m_\mathrm{p}=0\), we get

$$\begin{aligned} L_{k_\mathrm{p},{\mathbf {t}}}=\frac{{\widetilde{Q}}_{0}^{k_\mathrm{p}} (t_\mathrm{p})}{P_\mathrm{c}(t_\mathrm{c},t_\mathrm{p})}. \end{aligned}$$

Recall that the conditioning probability \(P_\mathrm{c}(t_\mathrm{c},t_\mathrm{p})\) is not affected by the process of sampling extant species. We get

$$\begin{aligned} L_{k_\mathrm{p},{\mathbf {t}}}&=\frac{\big (\mathrm {C}(f_\mathrm{p}) \,{\mathbf {Q}}^{k_\mathrm{p}}(t_\mathrm{p})\big )_{m=0}}{P_\mathrm{c}(t_\mathrm{c},t_\mathrm{p})}\\&=\frac{\sum _{n=0}^{\infty }\mathrm {C}_{m=0,n}(f_\mathrm{p}) \,Q_{n}^{k_\mathrm{p}}(t_\mathrm{p})}{P_\mathrm{c}(t_\mathrm{c},t_\mathrm{p})}\\&=\frac{\sum _{n=0}^{\infty }f_\mathrm{p}^{k_\mathrm{p}}\,(1-f_\mathrm{p})^{n} \,Q_{n}^{k_\mathrm{p}}(t_\mathrm{p})}{P_\mathrm{c}(t_\mathrm{c},t_\mathrm{p})}\\&=\frac{P_\mathrm{s}(t_\mathrm{c},{\mathbf {t}},t_\mathrm{p},f_\mathrm{p})}{P_\mathrm{c}(t_\mathrm{c},t_\mathrm{p})}, \end{aligned}$$

which is identical to Eq. (28). This ends the proof. \(\square \)

Note that Eq. 26 is equal to \(f_\mathrm{p}^{k_\mathrm{p}} F_{k_\mathrm{p}} (1 - f_\mathrm{p}, t_\mathrm{p})\) which provides an alternative route to prove Claim 6.1 (Manceau et al. 2019).

Finally, we give the expressions for the functions \({\widetilde{\xi }}(t,t_\mathrm{p})\) and \({\widetilde{\eta }}(t,t_\mathrm{p})\) in the case of constant rates, \(\lambda (t)=\lambda \) and \(\mu (t)=\mu \),

$$\begin{aligned} {\widetilde{\xi }}(t,t_\mathrm{p})&=\frac{f_\mathrm{p}\,\mu +((1-f_\mathrm{p}) \lambda -\mu )\,\mathrm {e}^{-(\lambda -\mu )(t_\mathrm{p}-t)}}{f_\mathrm{p}\, \lambda +((1-f_\mathrm{p})\lambda -\mu )\,\mathrm {e}^{-(\lambda -\mu )(t_\mathrm{p}-t)}}\\ {\widetilde{\eta }}(t,t_\mathrm{p})&=\frac{f_\mathrm{p}\,\lambda \big (1-\mathrm {e}^{-(\lambda -\mu )(t_\mathrm{p}-t)}\big )}{f_\mathrm{p}\,\lambda +((1-f_\mathrm{p})\lambda -\mu )\,\mathrm {e}^{-(\lambda -\mu )(t_\mathrm{p}-t)}}, \end{aligned}$$

which are identical to Eqs. (4) and (5) in the paper by Stadler (2013).

7 The Diversity-Dependent Case Without Extinction

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.

Proposition 7.1

Consider the diversity-dependent model without extinction, given by speciation rates \(\lambda _{n}(t)\). The diversification process starts at crown age \(t_\mathrm{c}\) with two ancestor species, and ends at the present time \(t_\mathrm{p}\), at which all extant species are sampled. Then, the likelihood of a phylogenetic tree with \(k_\mathrm{p}\) tips and branching times \({\mathbf {t}}\) is equal to

$$\begin{aligned} L_{k_\mathrm{p},{\mathbf {t}},0}^{\mathrm{(no-extinct)}}= (k_\mathrm{p}-1)!\,\prod _{i=2}^{k_\mathrm{p}-1}\lambda _{i}(t_{i}) \,\prod _{j=2}^{k_\mathrm{p}}\exp \left( \!-j\int _{t_{j-1}}^{t_{j}}\!\! \lambda _{j}(s)\,\mathrm{d}s\right) , \end{aligned}$$
(33)

where we used the convention \(t_{1}=t_\mathrm{c}\) and \(t_{k_\mathrm{p}}=t_\mathrm{p}\).

Proof

Equation (33) follows from Eqs. (2.4) and (2.5) in Ref. Rabosky and Lovette (2008), by noting that \(\xi _{i}\) in their notation corresponds to

$$\begin{aligned} \exp \Bigg (\!-\sum _{j=i}^{k_\mathrm{p}}\int _{t_{j-1}}^{t_{j}}\!\!\lambda _{j}(s)\,\mathrm{d}s\Bigg ) \end{aligned}$$

in our notation. \(\square \)

Note that in the case without extinction likelihood conditioning has no effect.

Theorem 7.1

Claim 3.1 holds for the diversity-dependent case without extinction.

Proof

To evaluate likelihood expression (6), we have to solve the ODE for \(Q_{m}^{k}(t)\), Eq. (3). Because species cannot become extinct and because all extant species are sampled, every species created during the process is represented in the phylogeny, i.e., there are no missing species. Hence, only the \(m=0\) component of \({\mathbf {Q}}^{k}(t)\) is different from zero. The ODE simplifies to

$$\begin{aligned} \frac{\mathrm{d}Q_{0}^{k}(t)}{\mathrm{d}t}=-k\lambda _{k}(t)\,Q_{0}^{k}(t), \end{aligned}$$

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_{0}^{k}(t_{k-1})\) at \(t_{k-1}\) , the solution is

$$\begin{aligned} Q_{0}^{k}(t)=Q_{0}^{k}(t_{k-1})\,\exp \Bigg (\!-\!k \int _{t_{k-1}}^{t}\!\!\!\lambda _{k}(s)\,\mathrm{d}s\Bigg ). \end{aligned}$$

At branching time \(t_{k}\), variable \(Q_{0}^{k}(t_{k})\) is transformed into variable \(Q_{0}^{k+1}(t_{k})\),

$$\begin{aligned} Q_{0}^{k+1}(t_{k})=k\lambda _{k}(t_{k})\,Q_{0}^{k}(t_{k}). \end{aligned}$$

Using the initial condition at crown age \(t_\mathrm{c}\), \(Q_{0}^{k=2}(t_\mathrm{c})=1\), we get

$$\begin{aligned} Q_{0}^{k_\mathrm{p}}(t_\mathrm{p})=\prod _{i=2}^{k_\mathrm{p}-1}i\lambda _{i}(t_{i}) \;\prod _{j=2}^{k_\mathrm{p}}\exp \Bigg (\!-j\int _{t_{j-1}}^{t_{j}}\!\!\lambda _{j}(s)\,\mathrm{d}s\Bigg ). \end{aligned}$$

Substituting into Eq. (6) yields the desired result. \(\square \)

8 Concluding Remarks

We have shown here that for the diversity-independent, but time-dependent birth–death 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 and Etienne 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 diversity-dependent 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.