Introducing a general class of species diversification models for phylogenetic trees

Phylogenetic trees are types of networks that describe the temporal relationship between individuals, species, or other units that are subject to evolutionary diversification. Many phylogenetic trees are constructed from molecular data that is often only available for extant species, and hence they lack all or some of the branches that did not make it into the present. This feature makes inference on the diversification process challenging. For relatively simple diversification models, analytical or numerical methods to compute the likelihood exist, but these do not work for more realistic models in which the likelihood depends on properties of the missing lineages. In this article, we study a general class of species diversification models, and we provide an expectation‐maximization framework in combination with a uniform sampling scheme to perform maximum likelihood estimation of the parameters of the diversification process.

Phylogenetic trees are types of networks that describe the temporal relationship between individuals, species, or other units that are subject to evolutionary diversification. Many phylogenetic trees are constructed from molecular data that is often only available for extant species, and hence they lack all or some of the branches that did not make it into the present. This feature makes inference on the diversification process challenging. For relatively simple diversification models, analytical or numerical methods to compute the likelihood exist, but these do not work for more realistic models in which the likelihood depends on properties of the missing lineages. In this article, we study a general class of species diversification models, and we provide an expectation-maximization framework in combination with a uniform sampling scheme to perform maximum likelihood estimation of the parameters of the diversification process.

INTRODUCTION
Evolutionary relationships of species are commonly described by phylogenetic trees or, in more general scenarios, by phylogenetic networks (Ragan, 2009). A phylogenetic tree is a hypothesis on how species or other biological units have diversified over time. It is usually described by a binary tree whose nodes are ordered in time. Phylogenetic relationships can be inferred from a variety of sources such as morphology and behaviors of species, biochemical pathways, DNA, and protein sequences (Lemey, Salemi, & Vandamme, 2009), both from extant, that is, living species or from extinct species through ancient DNA or the fossil record. However, data on extinct species are often incomplete and only accurate molecular phylogenies of extant species are available. In this article, we consider such phylogenetic trees as primary observations. Even though they lack extinct lineages, they are believed to contain information on how species diversified and hence they have been used to answer fundamental questions, such as "does diversity affect diversification?" (Cornell, 2013;Etienne et al., 2012), "what is the effect of environmental and ecological interactions on evolutionary dynamics?" (Barraclough, 2015;Ezard, Aze, Pearson, & Purvis, 2011;Lewitus & Morlon, 2017), "how does biodiversity vary spatially?" (Goldberg, Lancaster, & Ree, 2011;Mittelbach et al., 2007), and "what traits play a key role in species diversification?" (FitzJohn, Maddison, & Otto, 2009;Lynch, 2009;Paradis, 2005), to name just a few.
To help to answer these questions, specific mathematical models have been developed that can infer various parameters from phylogenetic diversification pattern (Morlon, 2014). Most current approaches have started to use likelihood-based methods to perform inference on phylogenetic trees FitzJohn et al., 2009;Ricklefs, 2007;Stadler, 2011). Although statistically principled, in each of these models, a new method to compute the likelihood needs to be developed. These models often rely on describing the macroevolutionary process by coupled ordinary differential equations-the so-called master or Kolmogorov equations-and these quickly become intractable as model complexity increases, particularly due to the lack of data on extinct species (Höhna, Stadler, Ronquist, & Britton, 2011;Ricklefs, 2007).
Alternative ways to deal with Kolmogorov equations have been used since the 1950s in fields outside evolutionary biology. These methods have used point process theory (Daley & Vere-Jones, 2007;Serfozo, 1990), which does not solve Kolmogorov equations directly but employs Gillespie-type simulations that were introduced in the context of chemical reaction modeling (Gillespie, 1976(Gillespie, , 1977. A single Gillespie simulation represents an exact sample from the probability mass function that is the solution of the system, thus allowing for stochastic optimization methods to maximize the likelihood (Tijms, 1994).
In this article, we present a first step for a general inference procedure of a general species diversification model. In Section 2, we describe a general diversification process based on a generalized linear model (GLM) description of a nonhomogeneous point process. This model can be used to describe many alternative evolutionary hypotheses. In Section 3, we introduce an expectation-maximization (EM) algorithm to optimize the likelihood under incomplete information, namely, the extinct lineages. We present a data augmentation algorithm, involving stochastic simulation combined with an importance sampler, to perform the E-step. We provide a proof-of-concept by comparing our inference with that obtained using direct likelihood calculations. In Section 4, we apply our method to the diversification of a small clade of Vangidae, consisting of a group of medium-sized birds living in Madagascar. Our aim is to discover Rampal Etienne and Ernst C. Wit jointly contribute to last authorship. whether the evolutionary record supports more the diversity dependence hypothesis  or the phylodiversity hypothesis (Castillo, Verdú, & Valiente-Banuet, 2010), for which no direct likelihood computation exists. Finally in Section 5, we provide directions for future extensions of the method that are needed to allow evolutionary biologists to routinely apply our approach to larger phylogenetic trees to study general diversification dynamics in a unified framework.

A GENERAL DIVERSIFICATION MODEL
We define a phylogenetic tree x = ( , t, a) on a time interval [0, T] as a functional object described by three components: a binary vector of event types (speciation or extinction), a vector of continuous event times t, and a network configuration object a, describing which species speciated or went extinct at each event time. We model the shape and structure of the tree by means of a collection of point processes, in this case, a set of dynamical nonhomogeneous Poisson processes (NHPP) where speciation and extinction of species are random events that happen within a time interval [0, T]. Figure 1 shows an example of a phylogenetic tree with three speciation events and one extinction event.
In this article, we assume that the process starts at time t 0 = 0 with a single species b 1 . At this stage, the tree is subject to two Poisson processes: a potential speciation of species b 1 and a potential extinction of species b 1 . Both processes are assumed to have a waiting time with time-continuous rates b 1 (t) and b 1 (t), respectively. In the time-homogeneous case, the waiting time for the first event to occur is therefore an exponential with rate b 1 + b 1 . More in general (Daley & Vere-Jones, 2007), the probability density for the process x to have a single species up tox time t 1 and a speciation event exactly at time t 1 is given by (1) F I G U R E 1 Phylogenetic tree with four events: three speciation events and one extinction event. Each branch represents a species If indeed a speciation occurs, the process continues with four NHPPs: two potential speciations and two potential extinctions. This is repeated until the present time T, unless the tree dies out before then. We consider a general scenario where at time t each of the N t present species b has its own speciation rate b (t) and extinction rate b (t) defined as a linear function via link function h, where c bjt is one of j = 1, … , m possible covariates of species b at time t affecting the speciation and/or extinction processes. Our entire process is therefore governed by the parameter set Typically, we will consider the logarithmic link function h = log, but Equation (2) can be trivially modified by choosing for h any monotonous increasing function that maps (0, ∞) onto R. The class of statistical models satisfying these specifications are an extension of the well-known GLMs (Dobson & Barnett, 2008).
This GLM extension to phylogenetic trees spans a very broad spectrum of possibilities for evolutionary biologists to test hypotheses and integrate their species diversification data. Diversification rates can be influenced by individual attributes, typically called traits, environmental factors, such as average temperature, by the composition of the diversifying clade itself or of its local ecological community. In the literature, a range of models have been explored, where diversification rates are assumed to be constant (Nee, May, & Harvey, 1994), change through time (Rabosky & Lovette, 2008), depend on diversity , on individual traits (Freckleton, Phillimore, & Pagel, 2008;Paradis, 2005) or other factors (Morlon, 2014). In order to test realistic models, we are interested in flexible rates that are able to change dynamically through all those factors simultaneously. For example, the speciation rate of species b at time t could also depend on other species' traits.
Mathematically, the method allows the inclusion of any set of covariates that might be interesting to incorporate for evolutionary biologists; however, full information on individual covariates, such as traits, are rarely available-especially not on the missing species. One way to deal with this is by including an extra augmentation step and simulating full information of traits on augmented trees (Hoehna et al., 2019). Another option is to use observable proxies related to, for example, trait diversity, such as different forms of phylogenetic diversity. These present interesting direction for future work.

MLE INFERENCE WITH MCEM USING IMPORTANCE SAMPLING
The loglikelihood of a full tree including extinct branches x ∈  involving a total of M events by extrapolating from (1) can easily be shown to be given by where be added to the likelihood, if the final event time t M does not correspond to the present T. For the case when diversification rates are stepwise constant, this reduces to the solutions in Wrenn (2012) and Reynolds (1973). When the full phylogenetic tree and the covariates at all times are given, we can directly maximize the loglikelihood function (3) to obtain the maximum likelihood estimates of the parameters (Paradis, 2005) and perform model selection to determine what factors are important for diversification. In practice, however, we almost never observe the full phylogenetic tree, but only a tree with the extant species.

Difficulties of MLE estimation and an MCEM algorithm
Let us denote  as the space of ultrametric trees (Gavryushkin & Drummond, 2016), that is, time-calibrated trees without extinct lineages and (y) as the space of all full trees that, when pruning all extinct species, lead to the ultrametric tree y ∈ . Then the log likelihood of an observed, extant species only tree y is given by the integral of the likelihood (3) over all possible full trees, However, because of the complexity of the space (y) a closed-form solution for Equation (4) is not available in most cases (Gavryushkin, Whidden, & Matsen, 2016), making inference, or in particular, direct MLE computations difficult or impossible.
A typical method for likelihood maximization under incomplete data is the application of the EM algorithm (Dempster, Laird, & Rubin, 1977), considering the information about the extinct species as a missing data problem. In the EM algorithm, a sequence { (s) } of parameter values are generated by iterating the following two steps: This algorithm is run iteratively until convergence is reached. Under certain regularity conditions (Dempster et al., 1977), the point of convergence can be shown to be the MLE for the incomplete data problem, that is, maximizing y ( ).
As in the case of Equation (4), the calculation of Q( | * ) does not have a closed-form due to the complexity of the space (y), so approximations are needed. To perform this task, we use Monte Carlo integration (Wei & Tanner, 1990), where given a set of sampled trees x 1 , … , x p from an importance sampler distribution g(x|y, ) we approximate Q( | * ) by where the importance weights are defined as w i = f X,Y (x i ,y| * ) g X|Y (x i |y, * ) , using the law of conditional probabilities to obtain the proportional expression.
In the M-step, we optimize (5) via numerical methods and the Hessian is calculated and represents the Fisher information matrix H −1 . Assuming that the errors given by the EM algorithm are independent of the Monte Carlo errors, the standard errors for the MCEM algorithm are defined as where −H −1 i,i corresponds to the diagonal components of the information matrix giving the EM error, VMCE is the variance of the MC error, and N EM is the number of MCEM iterations considered for estimation. Note that if the EM is run long enough, the second term in (6) goes to zero, making the information matrix the decisive value for standard errors on MCEM algorithm (McLachlan & Krishnan, 2007). With the standard errors we can construct confidence intervals for the parameters and test hypotheses about the significance of covariates of interest.

A simple importance sampler
To sample trees we propose a tree augmentation algorithm that samples independently the three components of the tree: event types, event times, and species allocations. The algorithm is shown in Figure 2.

Step 1. Generate event times and number of extinctions
The number of extinct species d and 2d missing event times, that is, speciations and extinctions of these d missing species are sampled uniformly in the following manner: 1. Sample the number of missing species d uniformly from the discrete space {0, … , M e }, where M e is a predefined ceiling, such that the probability of more than M e extinctions is extremely unlikely. 2. Sample 2d branching times uniformly from the continuous space (0, T] and then sort them. The probability of sampling a set of 2d unobserved event times t e = (t e 1 , … , t e 2d ) for a tree of dimension d is Note that this scheme samples the dimension of the tree uniformly, but the size of the space of trees grows in a factorial way with the dimension of the tree. This means that the sample size required to obtain a robust Monte Carlo approximation of the integral (4) must be large. This is a limitation of this importance sampler, and hence it is only reliable when many extinctions are unlikely.

3.2.2
Step 2. Generate event types We simulate a binary event chain e = ( e 1 , … , e 2d ) assigning either S (speciation) or E (extinction) to each event time. This chain is subject to the rule that the number of Es up to any point in F I G U R E 2 The three components of our phylogenetic tree augmentation algorithm

Input: Observed Phylogeny
Step 1: Draw number of events (2d), and then 2d event times Step 2: Draw event types (Dyck word) Step 3 : Allocate missing species the chain should be less than or equal to the number of Ss in the chain up to that point. The set of allowed chains is known in the mathematical literature as the set of Dyck words and several methods for sampling Dyck words have been developed (Kasa, 2010). Furthermore, given a number of events 2d, the number of possible Dyck words is known as the Catalan number (Zvonkin, 2014), By uniformly sampling a Dyck word e of length 2d, the probability of a specific event sequence is given by g events ( e ) = 1∕C d .

3.2.3
Step 3. Species allocation Given the missing event times and missing event types we can perform the tree allocations by sampling a parent species of each missing speciation and by defining which species, that is, the parent species or the inserted "new species," becomes extinct at the extinction event. To sample uniformly we just need to count the number of possible trees in agreement with the event times t e = (t e 1 , … , t e 2d ) and event types. This number, n( e 2d , t e 2d ), can be calculated by starting with n( e 0 , t e 0 ) = 1 and applying the following rules when going from root to tips in the phylogenetic tree: • For each unobserved speciation event at t e i , that is, e i = S, update n( e i , t e i ) in the following way, where N o t − is the number of observed branches just before t and N e t − is the number of unobserved branches just before t. Note that events on observed branches count twice compared with those on unobserved branches. Intuitively, this accounts for the two eventualities following an unobserved speciation on an observed branch: either the first or the second daughter species is observed (the other one is unobserved), while for a speciation on an unobserved branch, both daughter species are unobserved. A more formal argument justifying the factor of two is provided by Laudanno, Haegeman, and Etienne (2019).
• For each unobserved extinction event at t e i , that is, e i = E, update n( e i , t e i ) in the following way, As we sample uniformly, the probability for each possible allocation a e of the d missing species at the missing event times t e with Dyck word e in the tree of extant species x obs is then given by g allocation (a e ) = 1 n( e 2d ,t e 2d ) .

Sampling probability of a uniformly augmented tree
The uniform sampling probability of the augmented tree x unobs = (d, t e , e , a e ) is then given by From this equation, we can see how the dimension of the tree space plays an important role. For this reason, the uniform importance sampler becomes less efficient when many extinctions are likely. On the other hand, the uniform sampling scheme allows for easy implementation and quick computation, thereby making it suitable as a default sampler. Note: The first column is the mean of the effective sample size over the 1,000 iterations considered. The last row is the MLE directly calculated by computing the likelihood . Estimated values are for the linear DD model with 1 = 0 , 2 = ( 0 − 0 )∕K and 3 = 0 . 2012). In this model, speciation rates depend on diversity of the phylogenetic tree at that point.

Checking performance by comparing with direct ML
We consider the diversification model with rates where N t is the number of extant species (diversity) at time t and = { 0 , 0 − 0 K , 0 } are model parameters. This model is a special case of our general modeling framework, defined in (2). We perform the MCEM routine on a clade of Malagassy birds, the so-called Vangidae clade shown in Figure 3, which has been analyzed in Jønsson et al. (2012). We replicated the routine several times with different sample sizes to observe the impact of sample size on estimation and the robustness of the method.
In Table 1 we show six replicates corresponding to three pairs with different sample size orders. We drop the first 1,000 iterations as burn-in, and use the next 1,000 MCEM iterations for parameters estimation, reporting the mean value and the standard error from Equation (6). We observe that for small sample sizes (Replicates 1 and 2), the estimation is poor. For the cheapest setup, the mean effective sample size (ESS) is approximately 37 and this does not seem enough to sample in spaces with a substantial number of missing species. In this scenario, the MCEM estimates are not robust. As sample size increases, we see that inference becomes more and more accurate and matches the MLE procedure by Etienne et al. (2012). MCEM iteration Initial speciation rate (log) MCEM parameter estimation for three different sample sizes F I G U R E 4 MCEM applied to the tree of Figure 3 under the LDD diversification model. Evolution of the estimate of the first parameter, the initial speciation rate 1 = 0 through EM iterations. We plot six replicates: two for three different sample sizes. For better visualization we cut higher sample sizes at iteration 2,000 and 2,500 These replicates are also summarized in Figure 4, where we show a visualization of the dynamical MCEM parameter estimation for log 0 corresponding to the logarithm of the initial speciation rate at stem age. The dashed black line indicates the true MLE. We see in all six cases that estimations go quickly to the true MLE with a stable behavior after a couple of 100 iterations. To visually compare biases and variation through different sample sizes, we show the replicates for small sample sizes until the 2,000th and 2,500th MCEM iteration. We clearly see that for higher sample sizes, bias and variation decrease.
Note that the ESS is between 30% and 40% in these cases. An efficient importance sampler with 100% ESS is a priority for future publications in order to apply the method to larger phylogenetic trees.

DIVERSITY-DEPENDENCE: DIVERSITY OR PHYLODIVERSITY?
Phylodiversity is defined as the total branch length of extant species of a tree, and it has been proposed as an alternative to diversity in conservation ecology (Faith, 1992). Figure 5 shows phylodiversity and diversity through time for a simple example tree.
As an illustration of the flexibility of our method, we now consider a model similar to diversity-dependence introduced in the previous section, but with dependence on phylodiversity P t instead of N t . Diversity-dependence has been detected in a Vangidae clade (Jønsson et al., 2012) and we would like to extend the analysis to check if phylodiversity-dependence (LPD) is a more suitable factor in diversification of Vangidae than diversity-dependence (LDD). In addition to these two models, which both assume linear dependence of speciation rate on diversity or phylodiversity, we consider the exponential diversity dependence (EDD) and exponential phylogenetic diversity (EPD) models. The exponential models use the log-link function common in the statistical literature, rather than the identity link suggested by the evolutionary biology literature. Table 2 shows the parameter definitions for the four models tested on the phylogenetic tree of the Vangidae.
Note: All models assume constant extinction rate and have three parameters to be estimated. Abbreviations: EPD, exponential phylogenetic diversity; LDD, linear diversity-dependence; LPD, linear phylodiversity-dependence.
We performed the MCEM routine for each of the four diversification models, obtaining the ML estimates of the parameters and calculating Monte Carlo estimation for the likelihood function and the corresponding AIC values (Wit, Heuvel, & Romeijn, 2012). Interestingly, we found that phylodiversity models do not performs better than ordinary diversity models, but there is an improvement of the exponential diversity-dependence model over the linear DD model. Table 3 shows the inference results for each of the four diversification models.
To get an idea of the computational cost of the method we include, next to Table 3, a plot of computing times (for one MCEM iteration) as a function of Monte Carlo sample size for PD and DD models starting at their respective MLE values reported in the table. The values are average of 100 replicates performed in an ordinary computer. From the plot we can see that for our example tree, each iteration takes a couple of minutes for large Monte Carlo sample size, which means that the whole routine should take few hours at most. We also see that the computing times increases linearly with the MC sample size.
We conclude that the best model in this analysis is an EDD model with parameters 1 = 2.58(0.96), 2 = −1.02(0.25), 3 = 0.04(0.03), suggesting an exponential decreasing speciation rate with an exponential decay constant close to 1, given by 2 . We found an initial speciation rate of approximately 4.85 species per million years which decreases until 0.03 at the present time. This indeed suggests that the diversification process of this Vangidae clade in Madagascar has slowed down dramatically over the past 10 million years. Moreover, the extinction rate of 0.04 species per million years suggests that the clade has now reached a stable diversification behavior, whereby any further speciations will tend to be offset by extinctions. TA B L E 3 Parameter estimation of the four diversity-dependent models of Table 2  Note: Next to the table we see the plot of average computing times per MCEM iteration (in seconds) for DD and PD models at their respective MLE. Abbreviations: EPD, exponential phylogenetic diversity; LDD, linear diversity-dependence; LPD, linear phylodiversity-dependence.

DISCUSSION
We have presented a flexible method for testing a broad variety of diversification models in phylogenetic analysis and provided some simple examples. This is a first step toward a robust general methodology to identify potential factors in diversification processes from phylogenetic trees. The unobserved extinct species turn the inference problem naturally into a problem that can be approached by means of an EM algorithm. Given the complexity of the E-step, a Monte Carlo importance sampler has been proposed, involving a uniform importance sampler. Given the computational simplicity both in terms of sampling and calculation of uniform samplers, this may be a convenient option for small sized trees, where more sophisticated importance samplers, involving the underlying nonhomogenous Poisson processes, would not necessarily improve efficiency. As in the case of Vangidae clade where a few missing species are likely, we found that the uniform importance sampler leads to accurate estimation. However, the performance of our uniform importance sampler deteriorates as the dimension of the phylogenetic tree increases. In order to apply this method on high-dimensional trees, a more efficient importance sampler should be carefully chosen. This we will leave for future work. Current approaches perform inference by means of likelihood maximization, which requires that formulas for the likelihood must be derived on a case-by-case basis. Here, we consider a general class of models that include an augmentation step inside an EM algorithm, thereby avoiding direct likelihood calculation and thus allowing inference for a wide variety of diversification models.
In principle, in cases when full information of covariates is still missing after the augmentation step, extensions of the augmentation procedure are possible. However, this is beyond the scope of the current article.
Moreover, to increase efficiency alternatives to MCEM algorithms may be considered, such as the stochastic approximation version of the EM algorithm (SAEM; Delyon, Lavielle, & Moulines, 1999) or a Bayesian approach (Richardson & Green, 1997). In both cases, the algorithm could make use of the previous MC samples, thereby improving efficiency at some computational cost.
Even though in this article we only refer to the context of a diversification process of ecological species, a phylogenetic tree is used in many other fields to describe other kinds of processes, such as language evolution (Greenhill, Atkinson, Meade, & Gray, 2010) and cultural diversification