Parameter Identifiability of a Multitype Pure-Birth Model of Speciation

Species diversification models are used in Biology to make inferences about historical speciation and extinction rates over the time since a group of species, or taxa, evolved from a common ancestor. By providing information on rates of speciation and extinction, inference with these models seeks to give insight into the evolutionary dynamics leading to the present diversity of life. These models have a long history, starting with Yule’s constant-rate pure-birth model [18], and a fairly large literature has developed.

Diversification models describe a process beginning with a single lineage at some time in the past, which as time progresses may speciate or go extinct. When a speciation occurs the edge bifurcates into two edges, with the number of lineages increasing by 1. When an extinction occurs, the lineage ends, and the number of lineages decreases by 1. After either event, the process continues forward, independently on all lineages, producing a growing tree structure until the present time is reached. This tree, which has both topological and metric structure, constitutes an observation. (In applications, it may be necessary to consider the reconstructed tree, which is obtained by removing all tree edges with no descendents at the present [12, 5].)

Two basic sorts of these models have found common use in empirical studies. In the first, the speciation and extinction rates are functions of time, and apply to all taxon lineages present at any moment. This can be thought of as modeling exogenous factors, such as environmental conditions, that affect all taxa in the tree identically. Since all lineages behave in the same probabilistic way at any moment, it is not hard to show that the exact branching pattern of the tree-structure is irrelevant, with all the information in a tree observation being captured by the number of lineages as a function of time. Thus the work on time-dependent birth-death models by Kendell is foundational [7].

In the second sort of diversification model, which we call the multitype birth-death tree model, lineages are assigned one of a finite number of types at each moment, with the model’s speciation and extinction rates dependent only on the type. Over time, however, species may change types at fixed switching rates. This models endogenous factors, such as a particular biological trait a taxon may possess, including, for instance, a morphological feature, behavior, or whether a particular gene is present and active in an organism. A given type might correlate with faster or slower speciation than another, and/or affect the extinction rate. For these models the branching structure of a tree observation does matter, as taxa present at a given time may each have different types, and thus different tendencies to speciate or go extinct.

The Binary State-specific Speciation and Extinction (BiSSE) model formalized the multitype framework for biological applications [10]. Multitype (MuSSE) and quantitative-type (QuaSSE) variants of the model were subsequently proposed [4]. Although these works assumed the type is observed for the extant taxa at the leaves of a tree, we consider the multitype birth-death tree model with no type information observable for any lineage at any time, as type observations are unnecessary for our results. Indeed, the usefulness of these models to infer correlation between observed types and diversification rates from data with type information for extant taxa has been called into question [14].

Many other diversification models have been proposed, combining or extending these basic frameworks, with [16] offering one review. New variants continue to be developed, e.g., [11, 17, 15, 3].

When these models are used for inference, the data is taken as a single tree assumed to show the true evolutionary relationships of the taxa. (In practice, this tree itself must be inferred, usually from sequence data using phylogenetic and/or phylogenomic methods which we do not discuss here.) Multiple trees which one can reasonably hypothesize were generated with the same parameter values are simply not available. If the tree is sufficiently large, researchers hope it provides enough information to infer the speciation and extinction parameters reasonably well. More precisely, it has been implicitly assumed that the inference is statistically consistent, in the sense that as the number of taxa increases toward infinity (i.e., the tree grows larger), the probability of inferring model parameters arbitrarily close to the generating ones approaches 1. Establishing such a result, however, requires showing identifiability of the model parameters: A distribution derived from an observation of a single tree has a limit, as the number of taxa approaches infinity, that uniquely determines all parameter values.

Of course a full proof of the statistical consistency of a particular estimator requires additional arguments. For instance, the standard results on the consistency of maximum likelihood assume the availability of multiple independent samples, and therefore cannot be applied. Leroux’s result on the consistency of maximum likelihood inference from a single sequence of observation from a Hidden Markov Model [8] is analogous to what is need for applications of these diversification models. Nonetheless, establishing parameter identifiability is the first step toward this goal.

Recent work has shown that the first type of diversification model, with time-dependent rates, does not in fact have identifiable parameters, despite its widespread use by empiricists [9]. This result, which holds even if one allows for identification to be based on arbitrarily many independent tree observations with the same underlying rate parameters, was compellingly illustrated by construction of examples of wildly different rate functions producing identical tree distributions. An instance of this lack of identifiability had in fact appeared earlier, in an argument in which speciation rates were modified and extinction rates set to zero without changing the model distribution [12].

Little work, however, has addressed identifiability questions for multitype birth-death tree models. The strongest results on parameter identifiability for a pure birth model focus on a tree’s topological features but assume the types of both leaf nodes and their parents are observed [13]. In biological applications, however, the type of a leaf of the tree may be observable, but the type of the parent nodes is virtually never known. Thus no identifiability result relevant to typical data analyses has been produced.

One might hope that the analysis of multitype birth-death tree models would be simpler than for a time-dependent rate model, as its parameter space is finite dimensional. On the other hand, while trees produced by the time-dependent rate models can be summarized by the counts of lineages through time with no loss of information, this is not true for the multitype models. Effectively extracting information from a tree with both topological and metric structure requires a new approach.

In this paper, we investigate parameter identifiability of the multiype pure-birth tree (MPBT) model with any finite number of types. We thus restrict extinction rates for all classes to be zero. This model has also been called the multitype Yule model [13]. We assume only that the metric tree is observable, with no information on the types either at points internal to the tree or at the leaves. More formally, we establish generic identifiability of parameters up to label swapping. “Generic” means the result holds if we exclude parameters lying in a measure-zero subset of the parameter space. We give an explicit characterization of such a measure-zero exceptional set, as the zero set of a certain polynomial. “Up to label swapping” means that there are certain symmetries of the parameter space, arising from interchanging types, so that their corresponding speciation and switching rates are also interchanged, that have no effect on the model’s behavior. Generic identifiability up to label swapping is often the strongest form of identifiability that holds in models with hidden variables [1], and since we treat the types as unobservable, its appearance here is not surprising.

Our arguments draw on several earlier works. The first is [2] on Multitype Continuous Time Markov Branching Processes. In fact, these models and the MBPT model have the same underlying structure. But much of the classical branching process literature allows only for observing type counts over time, and not for observing the tree structure indicating the branching of specific lineages. The MPBT model, in contrast, treats the tree structure as observable, with type information hidden. Thus while providing an important tool in this work, the results of [2] are not immediately applicable to the MPBT model.

The second result crucial to our work is a general theorem on identifiability up to label swapping of parameters of a mixture model of product distributions [1]. In applying this to the MPBT model, we consider the joint distribution of edge lengths around a node on a uniformly-at-random chosen edge of a random tree, as the random tree grows arbitrarily large. Due to conditional independence of edge lengths conditioned on the type of the shared node, this joint distribution takes the form of a mixture distribution (over types) of product distributions. Although additional work is necessary to show parameter identifiability, this theorem is a crucial ingredient in our argument.

Although we do not address the multitype birth-death tree model with non-zero extinction rates here, we believe that our approach provides a pathway toward a more general result.

This paper is structured as follows. In Section 2 we provide a more formal definition of the MPBT model, and begin its analysis by deriving formulas related to the generation of a single edge in the tree in Section 3. Section 4 uses the results in [2] to obtain asymptotic results on the distribution of types across lineages in the tree at times increasingly distant from the root of the tree. Then, in Section 5, we bring these ingredients together, and apply the theorem of [1] to obtain our main results. Concluding remarks appear in Section 6.

In this section we formalize the Multitype Pure-Birth Tree model, in a form useful for our analysis.

Let $m$ be a positive integer denoting the number of types, and denote the set of types by $[m]=\{1,2,\dots,m\}$.

The parameter space of the MBDT model with $m$ types is all 3-tuples $(\boldsymbol{\pi},\boldsymbol{\lambda},S)$ described as follows:

A root distribution $\boldsymbol{\pi}=(\pi_{1},\pi_{2},\dots,\pi_{m})$, with $\pi_{i}\geq 0$, $\sum_{i}\pi_{i}=1$ gives probabilities $\pi_{i}$ of type $i$ being chosen for the tree root. A vector $\boldsymbol{\lambda}=(\lambda_{1},\lambda_{2},\dots,\lambda_{m})$ with non-negative entries gives speciation rates $\lambda_{i}$ for type $i$. An $m\times m$ matrix $S=(s_{ij})$ with non-negative off-diagonal entries and rows summing to 0 gives scalar type switching rates $s_{ij}$ from type $i$ to type $j$, $i\neq j$. Note that $S$ is determined by the $m^{2}-m$ independent scalar switching rates.

For parameters $(\boldsymbol{\lambda},S)$, let $D=\operatorname{diag}(\boldsymbol{\lambda})$ and $U=S-D$, so that the edge process $E_{\tau}$ has Markov rate matrix

For technical reasons we impose the following assumption, which is also biologically plausible.

Another ingredient of our approach to establishing the identifiability of MPBT model parameters is an analysis of an associated classical branching process, in which only the type counts are observed. More specifically, it records the number of edges of the tree which have each type as a function of time, but retains no information on the topology of the tree. We call this the type counting process, and in this section use established results to determine the asymptotic behavior of the relative frequencies of each type.

The asymptotics of the relative frequencies follow from results of [2] on multitype continuous-time Markov branching processes, specifically Theorems 1 and 2 of that work, which are paraphrased below as Theorem 4.5. Such a model can be described as a process where individuals of type $i$ live an exponentially-distributed length of time (whose rate only depends on type) and on death may be replaced by individuals of any type according to a distribution over $(\mathbb{Z}^{\geq 0})^{m}$.

To place the type counting process of the MPBT model into this framework, both speciation and change in type are viewed as deaths. Speciation results in replacement by 2 individuals of the same type, and change in type results in replacement by an individual of a different type. Since a speciation “death” of a type $i$ individual occurs with rate $\lambda_{i}$, and a type change “death” of a type $i$ individual followed by replacement with type $j$ occurs with rate $s_{ij}$, the combined rate of death for type $i$ is $\lambda_{i}+\sum_{j\neq i}s_{ij}$. When a death occurs, it is a speciation with probability

and a change to type $j$ with probability

Basic properties of the type counting process are summarized in the following.

We introduce yet another matrix defined in terms of the MPBT model parameters, as its leading eigenvalue and corresponding eigenvector plays a large role in the counting process’s behavior.

The matrix $A$ is the infinitesimal generator of the conditional expectation of the $N_{i}$s. More precisely,

with

where $\mathbf{e}_{i}$ is the $i$-th standard basis vector.

We will shortly show $\omega$ and $\mathbf{u}$ are uniquely determined, under an additional assumption.

Key properties of the counting process follow from the following more general theorem on classical branching processes.

Using the distributions of edge lengths and relative frequencies of each type of edge in a tree at a given time found in Sections 3 and 4, we are ready to establish identifiability of the MPBT parameters. To do so, we consider an asymptotic joint distribution of the lengths of 3 edges around a common node in the tree (see Figure 2). We seek to show that from this distribution the model parameters $(\boldsymbol{\lambda},S)$ can be determined, up to label swapping.

Due to the conditional independence of the lengths of three edges sharing a common node, given that node’s type, this distribution is a mixture of product distribution, with the mixing distribution and the components of the products closely related to distributions previously computed. This structure allows for the application of the following theorem, to obtain unmixed distributions of edge lengths conditioned on the type of the parental node. Thus even though we have no observation of type at any point in the tree, we can extract a distribution that is conditioned on type.

The following is a variant of Theorem 8 of [1], with the hypotheses modified as discussed on p. 3116 of that paper.

To apply this theorem, we make a further technical assumption, denoting the vector of $1$s by $\mathbf{1}$.

While the role of this assumption in our arguments will be clear in our proofs of Lemma 5.5 and Theorem 5.6 below, to understand its implications concretely, consider first the case $m=2$. Then

so

The non-singularity of $M$ thus is equivalent to $\lambda_{1}\neq\lambda_{2}.$ That these speciation rates would need to be different for parameters to be identifiable is intuitively clear, since otherwise type changes governed by $S$ would have no impact on the structure of the uncolored tree.

For general $m$, Assumption 3 is equivalent to the non-vanishing of $\det M$, a degree $\sum_{i=1}^{m-1}i=\binom{m}{2}$ polynomial in the $m^{2}$ independent entries of $\boldsymbol{\lambda},S$. Its non-vanishing thus excludes an algebraic variety of codimension 1, a set of Lebesque measure 0 in the unrestricted parameter space. An explicit calculation in the $m=3$ case shows the polynomial to be an irreducible polynomial in the $\lambda_{i}$ and $s_{ij}$, $i\neq j$.

The non-vanishing of $\det M$ always requires that the vector $\boldsymbol{\lambda}=-U\mathbf{1}$ not be a multiple of $\mathbf{1}$ (so that the first two columns of $M$ are linearly independent), and hence that not all $\lambda_{i}$ are the same. However, the additional restrictions it imposes on the parameters are more opaque to intuition without considering special cases.

For instance, when $m=3$, if all the $s_{ij}$ are equal, so the type switching behavior is identical for all types, the polynomial simplifies considerably, and factors as

Non vanishing of the polynomial, then requires that the three $\lambda_{i}$ be distinct, as one would expect is needed for identifiability, for otherwise several types would behave identically. However, for other choices of the $s_{ij}$, two of the $\lambda_{i}$ can be equal without the polynomial vanishing.

Next, we define the joint edge length distribution for several edges of a tree.

The three edge lengths used in the definition of $G_{T,t}$ are depicted in Figure 2, for $t=T/2$. The conditioning in the definition of $G_{T,t}$ ensures it only considers edges in which the edge process has led to speciation, that is, the edge processes for the parental and child edges are not truncated.

This immediately gives that $G_{\infty}$ is a finite mixture of product distributions.

In order to apply Theorem 5.1 to $G_{\infty}$, we need to verify that some of the univariate distributions in its decomposition above are linearly independent. To do so, the following lemma is needed.

We now introduce an additional assumption, which holds for generic parameters.

We now arrive at our main result.

Theorem 5.6, and Remark Remark, show that parameters $(\boldsymbol{\lambda},S)$ of the MPBT model can be identified from an asymptotic distribution as the tree depth grows, whether or not the number of sampled trees grows. Although this is not sufficient to conclude that estimation of parameters by maximum likelihood (ML) from a single tree, as suggested in [10], is statistically consistent, it does at least indicate that is a possibility. A similar question on ML inference of parameters for a hidden Markov model from a single sequence of observations was addressed by [8], with the consistency of ML estimation established as the sequence length goes to infinity.

For applications, it would be highly desirable to extend our identifiability result to a model incorporating constant extinction rates for each type. In most biological settings, the obtainable “data,” however is not the tree with edges stopping at extinction events, but rather the pruned tree in which all edges with no extant descendants are removed.

For a single type, parameter identifiability of a model with pruning was essentially considered in [12], where it was shown that the lineages-through-time function’s rate of change allowed the speciation and extinction rates to be determined, by separately considering the time regimes much earlier than the tree tips, and near the tree tips. An analysis combining the insight from [12] with the mixture distribution framework used in this work might be successful in showing parameters can be recovered from a single large tree observation for the multitype model.

Another interesting identifiability question for multitype tree models concerns what information on parameters is contained in the tree topology alone, or from weaker metric information than precise branch lengths. While our analysis depends heavily on metric features of the tree, that of [13] required no metric information. However, it did use type observations at the tips of the tree, and at their parental nodes. While types at tree tips may be observed in some biological studies, types of the parental nodes are generally not observable, as data is generally collected only from the taxa extant at the present. Even if ancient DNA or other trait data from earlier times is available, it is unlikely to be from the time of the last speciation.

ESA and JAR were supported in part by NSF Grant DMS-2051760.