The Fossilized Birth–Death Range Process

The FBD range process (illustrated in Fig. 1) was introduced in Stadler et al. (Reference Stadler, Gavryushkina, Warnock, Drummond and Heath2018). It begins with a single lineage at the origin time x ₀ and each lineage has an instantaneous speciation (or branching) rate λ and extinction rate μ, assuming a budding mode of species origination. Fossils are sampled along each lineage with rate ψ and extant species are sampled with probability ρ. We sample n species, of which m are extinct and l are extant. The total number of sampled fossil occurrences is denoted k. A given species i branches at time b_i in the phylogeny of sampled species (i.e., the tree in which all lineages leading to nonsampled species are removed) and goes extinct at time d_i. If species i survives to the present, we set d_i = 0. The first (or oldest) appearance of species i is sampled at time o_i and the last (or youngest) appearance at time y_i. The interval between o_i and y_i represents the stratigraphic range of species i. If a species is sampled only once (i.e., the species is a singleton), then o_i = y_i, and if species i is sampled only once at the present, then o_i = y_i = d_i. Figure 1A shows an example set of sampled species and how they relate to the true underlying tree. This is referred to as the “complete tree.”

Figure 1. The fossilized birth–death (FBD) range process. A, The complete FBD tree, including the observed data. Diamonds represent observed samples of fossils or extant species. Points labeled o_i and y_i are the first and last appearances of range i, respectively. B, The extended sampled stratigraphic ranges, illustrating the information used to calculate the likelihood. b_i and d_i are the species attachment and extinction times, respectively. The attachment time b_i is the time at which species i attaches to other lineages in an incompletely sampled tree, which is not necessarily equivalent to speciation time of species i. For example, the attachment time b ₂ for species 2 is not equivalent to the true species origin time shown in A, instead it is the origin time of its nonsampled ancestor, species 3. b ₅ = x ₀ is the origin time. In total we sample k = 5 occurrences and n = 4 ranges, with l = 2 extant and m = 2 extinct ranges. The dashed line highlights the number of possible attachment points (black circles) used to calculate γ_i, which accounts for each possible topology, given the value of b_i. γ_i is equivalent to the number of lineages that coexist with species i at time b_i. The exception is the oldest species, where b_i represents the origin and γ_i is always one. The number of possible attachment points for each species is γ₁ = 2, γ₂ = 2, γ₄ = 1, and γ₅ = 1.

To calculate the probability of observing a set of stratigraphic ranges, given the FBD parameters λ, μ, and ψ, we need to know the interval between the branching time b_i and the first appearance o_i, as well as the interval between o_i and extinction time d_i. Note that the interval between b_i and o_i is equivalent to a ghost lineage (Norell et al. Reference Norell, Novacek, Wheeler and Novacek1992). The interval between o_i and d_i is referred to as the “extended sampled stratigraphic range.” An example of the extended sampled stratigraphic ranges is shown in Figure 1B. Note that in the extended sampled stratigraphic ranges, the time at which species i attaches to the tree at time b_i may not correspond to the true speciation (origination) time of i in the complete tree, if species sampling is incomplete. Following Stadler et al. (Reference Stadler, Gavryushkina, Warnock, Drummond and Heath2018), we refer to the branching time b_i as the “attachment time.”

In reality, we typically only have information about o_i and y_i and do not know the underlying topology or the true branching times b_i and extinction times d_i. To account for uncertainty in the underlying topology, we need to define γ_i as the total number of coexisting lineages at time b_i, which represents all possible points at which species i can attach to the tree (see Fig. 1B). This value allows us to analytically integrate over all possible tree topologies (see Stadler et al. Reference Stadler, Gavryushkina, Warnock, Drummond and Heath2018). To account for uncertainty in the times, we can augment our data such that b_i > o_i and d_i < y_i and marginalize over all possible speciation and extinction times (b_i, d_i) using Markov chain Monte Carlo (MCMC). Summarizing ${\cal D} = \lpar {k\comma \;\;l\comma \;\;\lcub {b_i\comma \;\;d_i\comma \;\;o_i} \rcub \comma \;\;i = 1\comma \;\;\ldots \comma \;\;n} \rpar $, the probability of the augmented extended sampled stratigraphic ranges, as derived in Stadler et al. (Reference Stadler, Gavryushkina, Warnock, Drummond and Heath2018), is

(1)$$\eqalign{& \hskip -.5pc f\lsqb {{\cal D}\vert {{\rm \lambda }\comma \;{\rm \mu }\comma \;{\rm \psi }\comma \;{\rm \rho }\comma \;x_0} } \rsqb \cr & = \displaystyle{{{\rm \psi }^k{\rm \mu }^m{\rm \rho }^l{\lpar {1-{\rm \rho }} \rpar }^{n-m-l}} \over {{\rm \lambda }\lpar {1-p\lpar {x_0} \rpar } \rpar }}\mathop \prod \nolimits_{i = 1}^n \lambda \gamma _i\displaystyle{{{\tilde{q}}_{asym}\lpar {o_i} \rpar q\lpar {b_i} \rpar } \over {{\tilde{q}}_{asym}\lpar {d_i} \rpar q\lpar {o_i} \rpar }}\comma }$$

with

(2)$$\eqalign{p\lpar t \rpar & = 1 \cr & + \displaystyle{{-\lpar {{\rm \lambda }-{\rm \mu }-{\rm \psi }} \rpar + c_1_{{{e^{{-}c_{1^t}}\lpar {1-c_2} \rpar -\lpar {1 + c_2} \rpar } \over {e^{{-}c_{1^t}}\lpar {1-c_2} \rpar + \lpar {1 + c_2} \rpar }}}} \over {2{\rm \lambda }}}\comma \;}$$

(3)$$\tilde{q}_{{\rm asym}}\lpar t \rpar : = \sqrt {e^{{-}t\lpar {{\rm \lambda } + {\rm \mu } + {\rm \psi }} \rpar }q\lpar t \rpar } \comma $$

(4)$$q\lpar t \rpar = \displaystyle{{4e^{{-}c_{1^t}}} \over {{\lsqb {e^{{-}c_{1^t}\lpar {1-c_2}\rpar + \lpar {1 + c_2} \rpar} } \rsqb }^2}}\comma \;$$

(5)$$\eqalign{& c_1 = \left\vert {\sqrt {{\lpar {{\rm \lambda }-{\rm \mu }-{\rm \psi }} \rpar }^2 + 4{\rm \lambda \psi }} } \right\vert \comma \;\cr & c_2 ={-}\displaystyle{{{\rm \lambda }-{\rm \mu }-2{\rm \lambda \rho }-{\rm \psi }} \over {c_1}}.} $$

Because the FBD model incorporates the fossil- and extant species-sampling processes explicitly, all observed occurrences, including all extant samples, contribute to the likelihood calculation and may be included in the analysis.

The Birth–Death Sampling Process Assuming Complete Species Sampling

To examine the impact of incomplete species sampling we compared our approach with a birth–death model based on Keiding (Reference Keiding1975),

(6)$$f\lsqb {{\cal D}\vert {{\rm \lambda }\comma \;{\rm \mu }\comma \;{\rm \psi }} } \rsqb = {\rm \lambda }^B{\rm \mu }^D{\rm \psi }^ke^{-\lpar {{\rm \lambda } + {\rm \mu } + {\rm \psi }} \rpar S}\comma \;$$

where B is the total number of speciation (birth) events, D is the total number of extinction (death) events, and S is the sum of species durations $S = \;\mathop{\sum}_\limits_{i = 1}^n b_i-d_i$ (Foote and Miller Reference Foote and Miller2007). We refer to this as the BD model. This model allows for the inclusion of all observed fossil samples but assumes that each species was sampled at least once. However, the inclusion of extant singletons will lead to biased estimates, as modeling incomplete species sampling through time is required to account for the increase in diversity at the present if all extant samples are included. Thus, extant singletons should be excluded from the analysis. Like the FBD range model defined earlier, the BD model also assumes a constant rate of fossil sampling. The BD model described here is similar to the birth–death model implemented in the program PyRate (Silvestro et al. Reference Silvestro, Salamin and Schnitzler2014a), with the only difference being the PyRate BD model allows for rate variation across lineages.

Per-Taxon Rates and the Boundary-Crosser Method

Traditional, summary statistics-based approaches used in paleontology to estimate speciation and extinction rates are broadly based on quantifying the proportion of first and last appearances sampled during discrete geologic intervals. Following (Foote Reference Foote2000), fossils sampled during a given interval t_i may be categorized as:

1. 1. singletons: taxa confined to the interval, that is, taxa whose first and last appearance are both within the interval t_i (FL) (regardless of the number of occurrences found within that interval);

2. 2. bottom boundary crossers: taxa that cross the bottom boundary, that is, the boundary older than interval t_i, and make their last appearance during the interval (bL);

3. 3. top boundary crossers: taxa that make their first appearance during the interval and cross the top boundary, that is, the boundary more recent than interval t_i (Ft); and

4. 4. range-through taxa: taxa that range through the entire interval, crossing both the top and bottom boundaries (bt).

Examples of each category are shown in Figure 2. The first three categories represent species that must be sampled within interval t_i, while the fourth category (bt) can include taxa that are sampled before and after the interval, but not necessarily within the interval.

Figure 2. Boundary-crosser versus three-timer taxon categories illustrated for a given interval t_i. Boundary-crosser categories: Ft taxa first appear and bL taxa last appear within the interval going forward in time, while the first and last appearance of FL taxa is confined to the same interval and bt taxa appear before and after the interval. The boundary-crosser method uses the bL, Ft, and bt taxa only in the estimation of λ and μ. Three-timer categories: two-timers are sampled immediately before and within bin t_i (2t_i) or within and immediately after t_i (2t_{i +1}). Three-timers are sampled immediately before, within, and after t_i (3t), and part-timers are sampled immediately before and after but not within t_i (pt) (note all three-timers are also two-timers). Singletons (i.e., taxa sampled from a single interval or FL taxa) are excluded before analysis. Image adapted from Alroy (Reference Alroy2010, Reference Alroy2014).

The most straightforward approach to estimating speciation and extinction rates for a given interval using this information is to use the proportion of first and last appearances observed during the that interval,

(7)$$\eqalign{& {\rm \lambda } = \left({\displaystyle{{N_{FL} + N_{Ft}} \over {N_{{\rm tot}}}}} \right)\times \displaystyle{1 \over {\Delta t_i}}\comma \;\cr & {\rm \mu } = \left({\displaystyle{{N_{FL} + N_{bL}} \over {N_{{\rm tot}}}}} \right)\times \displaystyle{1 \over {\Delta t_i}}\comma \;} $$

where N denotes the number of representatives for a given category and N _tot =  N _FL +  N _bL + N _Ft + N _bt and Δt_i is interval length. This equation relies on Δt_i being a very short time interval (i.e., the probability of a single species speciating twice within this interval is negligible). Then, for the speciation process, the probability that a species does not speciate within time interval Δt_i is λΔt_i, which can be estimated by ${{N_{FL} + N_{Ft}} \over {N_{{\rm tot}}}}$. This approach is taken from Foote (Reference Foote2000: eqs. 12 and 13) and is an extension of theory presented by Raup (Reference Raup1985). Estimates are defined as “per-taxon origination or extinction rates”; we refer to this as the “per-taxon rates method.” This approach makes no attempts to mitigate the effects of incomplete sampling and allows for the inclusion of all fossil occurrences, including singletons. However, because rates are calculated for intervals before the present, extant singletons do not contribute to the estimation of rates (this also applies to the boundary-crosser and three-timer methods described below).

More commonly applied methods in paleontology use strategies to eliminate or minimize biases created by incomplete sampling. The boundary-crosser method (Foote Reference Foote2000: eqs. 22 and 23) uses the above categories but does not consider singletons (category 1) in the estimation of speciation and extinction rates,

(8)$$\eqalign{& {\rm \lambda } ={-}ln\left({\displaystyle{{N_{bt}} \over {N_{Ft} + N_{bt}}}} \right)\times \displaystyle{1 \over {\Delta t_i}}\comma \;\cr & {\rm \mu } ={-}ln\left({\displaystyle{{N_{bt}} \over {N_{bL} + N_{bt}}}} \right)\times \displaystyle{1 \over {\Delta t_i}}.} $$

The derivation of these equations is straightforward. For example, in the case of the speciation process, the probability that a species does not speciate within time interval $\Delta $t_i is $e^{-\lambda \Delta t_i}$, which can be estimated by ${{N_{bt}} \over {N_{Ft} + N_{bt}}}$. This approach can only include taxa that have been sampled during more than one interval and thus incorporates fewer data overall.

The Three-Timer Method

The three-timer method goes further in an effort to account for incomplete sampling by defining an alternative set of categories and applying a sampling correction (Alroy Reference Alroy2008; Alroy et al. Reference Alroy, Aberhan and Bottjer2008). The three-timer categories are defined as:

1. 1. two-timers: taxa sampled immediately before and within a bin (2t_i) or taxa sampled within and immediately after (2t_{i+ 1});

2. 2. three-timers: taxa sampled immediately before, within, and after a bin (3t_i); and

3. 3. part-timers: taxa sampled immediately before and after but not within the bin (pt_i), implying at least one unsampled lineage.

Notation and definitions follow Alroy (Reference Alroy2008, Reference Alroy2010), and examples are shown in Figure 2. The part-timer sampling probability is defined as,

(9)$$P_S = \displaystyle{{N_{3t}} \over {N_{3t} + N_{\,pt}}}\comma \;$$

where N _3t and N_pt may be summed across the entire dataset. The three-timer speciation and extinction rates are defined as

(10)$$\eqalign{& \lambda = \left({ln\left({\displaystyle{{N_{2t_{i + 1}}} \over {N_{3t_i}}}} \right)+ {\rm ln}\lpar {P_S} \rpar } \right)\times \displaystyle{1 \over {\Delta t_i}}\comma \;\cr & \mu = \left({ln\left({\displaystyle{{N_{2t_i}} \over {N_{3t_i}}}} \right)+ {\rm ln}\lpar {P_S} \rpar } \right)\times \displaystyle{1 \over {\Delta t_i}}.} $$

To obtain tree-wide estimates of λ and μ using all three paleontological methods, we summed the number of samples that belong to each taxon category across all intervals, rather than averaging rates calculated separately for each interval. This relies on the assumption that all time intervals (i.e., Δt) are of equal length, which we control for in our simulation design. This choice works in favor of these methods, as it increases the potential information used to calculate rates: some intervals may contain sampled species but do not have sufficient information to recover rate estimates (e.g., an interval may contain boundary crossers but no range-through taxa, or two-timers but no three-timers). We cannot estimate rates for such intervals, but we can use the sampled species in the summation used to calculate tree-wide rates. Note that the sampling parameter P_S used by the three-timer method accounts for incomplete sampling in the preceding (i.e., older) interval in the case of speciation and in the following interval (i.e., younger) in the case of extinction. Because P_S does not apply to extant species sampling, the interval before the present is excluded from tree-wide estimates of extinction rates using this method.

Note that in contrast to the phylogenetic model–based approaches, these three paleontological methods do not require information about the total number of sampled fossil occurrences (i.e., parameter k). Instead, these methods only require information about whether a given taxon was present or absent during each interval. We used simulated data to assess and compare the performance of alternative approaches to estimating tree-wide speciation and extinction rates.