Enforced symmetry: the necessity of symmetric waxing and waning

Temporal dynamics of geographic range size and occupancy

Systematic studies of the geographic expansion of taxa over their lifespan were spurred by Willis (1922), who postulated a constant increase of the taxon’s range size over its duration. This trajectory, known as the age and area hypothesis, has been met with many counterexamples. Among them, the concept of taxon cycle (Dillon, 1966; Ricklefs & Cox, 1972; Ricklefs & Cox, 1978) postulated phases of range expansion, establishment of maximum range, fragmentation due to local extinctions, and a final demise. Neither of these trajectories has been found to hold universally (Gaston, 1998; Pigot, Owens & Orme, 2012). Empirical studies revealed such a variety of temporal trends in geographic range sizes, both within and between taxonomic levels, that Gaston (1998) proposed that each species might show an idiosyncratic trajectory that does not allow to infer on any general pattern. Nonetheless, the duration of a taxon and its geographic range are reciprocally influencing each other, depending on which phase of the taxon’s duration is considered (Foote et al., 2008).

But apparent idiosyncrasy of trajectories of individual taxa need not mean that there are no shared mechanisms determining the relationship between range size and age of a taxon. On the contrary, this is consistent with the trajectories being produced by a stochastic process. It has been shown that a wide range of empirical patterns in the age-area relationships observed in diverse groups of organisms can be produced by a single stochastic process (Pigot, Owens & Orme, 2012). This approach took into account two essential parameters: (1) the relationship between the range size and the probabilities of extinction and origination, and (2) the asymmetry of range division during a speciation event (Gaston, 1998). Both parameters vary between taxa substantially and their variation yields all types of range size-age relationships (Pigot, Owens & Orme, 2012). Moreover, all types of trajectories considered in classifications outlined by Gaston (1998) and Liow & Stenseth (2007) can be produced by one stochastic process and their relative frequencies will vary according to the parameters of speciation dynamics and range area division. Notably, averaged simulated range size trajectories from origin to extinction show a symmetrical waxing and waning pattern even though this pattern is not even produced once as an individual trajectory in the simulation (Pigot, Owens & Orme, 2012). Therefore, the underlying stochastic process cannot be reconstructed from empirical trajectories which it produced.

Although the above discussion focuses on geographic range size, the same considerations are—with respect to average patterns—valid for other measures of geographic distribution, such as latitudinal range and various MESs employed in ecology and palaeobiology (for example  Jernvall & Fortelius, 2004; Raia et al., 2006; Foote, 2007; Foote et al., 2007).

Temporal dynamics of diversity

As a descriptive pattern, symmetry emerges in compilations of averaged diversity trajectories across various taxonomic and temporal scales (Foote, 2007; Foote et al., 2007; Quental & Marshall, 2013). It is also manifested in the averaged CG (center of gravity) measure introduced by Gould, Gilinsky & German (1987), which corresponds to the center of mass of the diversity of a taxon plotted against its relative lifespan. For most of the Phanerozoic and when averaged across multiple taxa, the CG metric is indistinguishable from 0.5 (Gould, Gilinsky & German, 1987; Plotnick & Wagner, 2018), in concordance with symmetric waxing and waning.

When phylogenetic relationships are not known, patterns of diversity consistent with those observed in the fossil record are modeled as a random walk (Cornette & Lieberman, 2004; Žliobaite, Fortelius & Stenseth, 2017). Evolutionary relationships can be represented using branching (or birth-death) processes, which also allow reproducing empirical patterns (Raup et al., 1973; Gould, Gilinsky & German, 1987; Nee, 2006; Foote, 2007; Quental & Marshall, 2009). Both models are widely adopted to test macroevolutionary hypotheses, which commonly requires constraining the parameters of the process, the origination and extinction rates. To this aim, statistical methods are necessary to generalize over realized rates seen in individual trajectories and to estimate the underlying rates of the common process generating these trajectories (Foote, 1988; Newman & Sibani, 1999). It should be noted that there are other methods available to estimate diversification rates which do not address commonalities in these rates within groups and thus avoid the problem of averaging the trajectories (Alroy, 2014; Foote, 2001; Silvestro et al., 2014; Sakamoto, Benton & Venditti, 2016; Liow & Finarelli, 2014). The approach using diversity trajectories is outlined here.

It is based on the intuition that analyzing a large set of diversity trajectories should reveal their shared properties. Therefore the symmetry emerging from averaging diversity trajectories across taxa has been interpreted as indication that both extinction and origination rates are equal and time-homogeneous (Quental & Marshall, 2013; Nee, 2006). From this, it also follows that rates calculated based on the averaged, symmetrical pattern have mirrored dynamics in the sense that the origination rate declines at the same pace as the extinction rate rises and the final value of one equals the initial value of another (see model A of Gilinsky & Bambach (1987) and Quental & Marshall (2013)). Moreover, the symmetric pattern has been used to argue that, at any moment in the Earth’s history, half of the taxa are diversifying and half are in decline (Quental & Marshall, 2010). There is, however, a discrepancy between the averaged pattern and the frequency of types of trajectories, for example individual symmetrical trajectories form only a marginal proportion of the datasets producing averaged symmetrical patterns (Quental & Marshall, 2013). The distribution of changes in the origination and extinction rates is skewed towards more common increases in the origination rates during the diversification phase of clades, although it is symmetrical during their decline phase. Similarly, origination declines are more common than extinction declines in diversity trajectories of marine taxa (Gilinsky & Bambach, 1987). This raises caution that the average symmetrical waxing and waning pattern does not allow us to conclude equal or time-homogeneous rates of the underlying branching process.

Predicting extinction risk

Inferring age-dependent extinction risk is a question parallel to the age-area relationship discussed in the previous point; it is also among the questions where paleobiological data may aid conservation decisions. It has been shown across numerous clades and timescales that occupancy, geographic range and other MESs are among key predictors of extinction risk (Finnegan, Payne & Wang, 2008; Kiessling & Aberhan, 2007; Payne & Finnegan, 2007; Finnegan et al., 2015). Kiessling & Kocsis (2016) showed that the change in occupancy through time is a better predictor of extinction risk than standing occupancy. This brings into light the possibilities offered for conservation by the fossil record and the trajectories of species presence through time (Harnik et al., 2012). The potential of this approach can be compromised by averaging trajectories of the MESs across entire groups to predict extinction risk in higher level clades (Tietje & Kiessling, 2013).

The way data are processed

The claim of this paper is directly linked to the way data are processed. We will refer to the following data-processing procedure to generate averaged trajectories:

1. Take the trajectories describing the presence of different taxa from their origination to their extinction

2. Shift each trajectory so the origination of the taxon is at time zero (x-shift)

3. Rescale each trajectory so the extinction of the taxon is at time one (x-scale)

4. Rescale the amplitude of all trajectories (y-scale)

5. Average all trajectories

This procedure is based on (Foote, 2007; Foote et al., 2007) and shown in Fig. 1.

One part of this data-processing procedure that will increase symmetry and similarity is averaging. This is demonstrated by Donsker’s invariance principle (Klenke, 2008, p. 474), which states that averaging different stochastic processes can lead to the same, symmetric stochastic process (section “Symmetry and Similarity by Averaging”). Another influence is the blurring effect of background noise that is being introduced into the analysis by the x-shift and the x-scale. This combines different sources of noise from different time intervals and of different taxa, which can drown possible signals (section “The Effect of Noise”). The last effect discussed here that can increase symmetry is conditioning. This describes the effect of already knowing the extinction and the origination, which puts an upper and a lower bound on the developments of the taxa in their lifespan, thereby making strongly asymmetric developments impossible (section “Symmetry by Conditioning”).

Measuring symmetry

Symmetry and asymmetry are easily tangible by intuition. To make our results quantifiable and reproducible, we define a measure of asymmetry named QuAsy (Quantified Asymmetry).

A function f is axis symmetric (with respect to the y axis) if f(x) = f(−x). The symmetric part of said function is the largest symmetric function smaller than f and is given by f_sym = min{f(x), f(−x)}. The asymmetrical part of f is then given by f_asym = f − f_sym, allowing to defined the QuAsy as (1)${\displaystyle QuAsy\left(f\right)=\int f_{asy}\left(x\right)\mathrm{d}x=\int f\left(x\right)-min\left\{f\left(x\right),f\left(-x\right)\right\}\mathrm{d}x,}$

see Fig. 2 for an example. For binned data, the QuaSy is given by (2)${\displaystyle QuAsy\left(f\right)=\sum _{i=1}^n\left[f_i-min\left\{f_i;f_{n-i}\right\}\right],}$ where f_i is the value of the function in the ith bin. This measure of asymmetry is slightly weaker than a seminorm, and it is zero if and only if the function is symmetric. Its properties are described in Appendix S1.

Symmetry and similarity by averaging

One way to generate symmetry out of randomness and reduce distinguishability is averaging. Here, we take distinguishability as the ability to statistically distinguish hypotheses by the sets of trajectories generated and similarity as the opposite of distinguishability. The increase of symmetry and similarity by averaging is best known from the central limit theorem, which states that the averaged sum of random variables converges towards the symmetric standard normal distribution (Laplace, 1812). This is even the case for the sum of random variables with different distributions (Lindeberg, 1922). Similar results exist for stochastic processes. One such result is Donsker’s invariance principle (Klenke, 2008, p. 474)). It states that averaged sums of stochastic processes will converge to Brownian motion, a stochastic process that is commonly named the Wiener process in probability theory. Donsker’s invariance principle can be seen as an infinite dimensional version of the central limit theorem, with Brownian motion replacing the standard normal distribution. Just like the standard normal distribution, Brownian motion is highly symmetric:

• It is self-similar in the sense that rescaled versions of Brownian motion are again Brownian motion (Klenke, 2008, p. 455)

• It is time reversible in the sense that if (W_t)_(t∈[0,1]) is a Brownian motion process, then so is (W₁ − W_1−t)_(t∈[0,1]).

Note that the reversibility of Brownian motion implies the axis symmetry of symmetric waxing and waning, as it is for example displayed in Fig. 2 of Foote (2007). This is because reversibility means that it is indistinguishable whether the measures of interest were recorded from origination to extinction or vice versa, which implies that the averaged trajectories look the same, independently of whether they are plotted from origination to extinction or from extinction to origination. In other words, averaging trajectories that are reversible results in symmetry along a vertical axis through the midpoint of the taxon’s life span. The idea of reversibility and its link to axis symmetry is further discussed in Appendix S2.

Donsker’s invariance principle is not only valid for a symmetric random walk, but also for the much broader class of stochastic processes with expectation value zero and finite variance (see Klenke (2008) p. 474 for technical details). Donsker’s invariance principle demonstrates two different effects of averaging:

1. Averaging increases similarity, since if averaged, different stochastic processes converge to the same limiting process

2. Averaging increases symmetry, since it can turn asymmetrical stochastic processes into symmetrical stochastic processes

So if the trajectories of any type of presence MESs of taxa are generated by some type of stochastic process, they will be affected by this averaging-induced increase in symmetry (see Figs. 3 and 4) and similarity.

It must be noted that Donsker’s invariance principle is based on the weak averaging factor $\frac{1}{\sqrt{n}}$ as it is used in the central limit theorem, which guarantees that the result will be a stochastic process. Many limiting theorems, such as the weak and strong law of large numbers (Klenke, 2008, ch. 5), are based on the stronger averaging factor $\frac{1}{n}$, making the limit a number, and not a distribution. Using this stronger averaging factor will therefore make the averaged trajectory smoother and emphasize the symmetry even more.

The effect of noise

Another effect that reduces distinguishability is the combination of noise from different sources, which can make it hard to detect signals.

For example, occupancy trajectories of individual taxa are influenced by a variety of factors, like environmental conditions or biotic interactions, affecting a given taxon throughout its lifespan. Both the nature and timing of these factors may differ among taxa. By rescaling the trajectories of a group of taxa to start at zero and to end at one (x-shift and x-scale in the data processing procedure described in section “The Way Data are Processed”), signals that were originally in different time intervals are then placed in the same bin between zero and one. Analyzing these rescaled trajectories combines all the different signals from different sources and different effects on different taxa in different times. This leads to a strong background noise that can drown any potential signal.

As an example, take any real valued random variable Y_n describing the result of a statistical analysis after n samples have been evaluated. A common example for such a random variable is the sample mean, given by (3)${\displaystyle Y_n=\frac{1}{n}\sum _{i=1}^n{X}_i,}$ used as an estimator for the expectation value and as the basis for the construction of confidence intervals. In a more general setting, the task of Y_n can be to identify the underlying probability distribution by the realizations it generates.

Now, assume that the statistical analysis is well behaving in the sense that its quality increases with sample size and the corresponding random variable converges towards some deterministic value a. For applications, this raises the question about the number of samples required for the results of the statistical analysis to be reliable in the sense of displaying only fluctuations below some threshold around the limit a. We will show that the number of samples necessary for a reliable result increases as noise is introduced. This is formalized by adding perturbation random variables Z_n to the results of the statistical analysis. The new, perturbated analysis is described by the random variables (4)${\displaystyle {\tilde{Y}}_n=Y_n+Z_n.}$ In Appendix S2, we show that even in the best case scenario where the Z_n converge to zero, meaning that noise gets weaker as sample size increases, the perturbated random variables ${\tilde{Y}}_n$ fluctuate more than the original random variables Y_n. This implies that the introduction of noise increases the number of samples necessary to obtain results of any chosen reliability. Conversely, for a fixed number of samples, results derived from a noisy scenario are less reliable than the ones from a scenario free of noise.

Note that no further assumptions on the type of statistical analysis, the type of noise or how they change with sample size were made.

Symmetry by conditioning

Another effect that increases symmetry is the incorporation of the knowledge of the already known origination and extinction times of taxa. In probability theory, this is known as conditioning, and is based on the conditional probability. One example for conditioning can be found in Bayesian statistics, where the posterior is defined as the prior, conditioned on the observation. This incorporation of knowledge leads to a re-evaluation of the probabilities and the feasibility of trajectories, thereby putting boundary conditions on the trajectories of taxa:

• The MESs of a taxon cannot decline too much in the early part of its life, since otherwise it would go extinct

• The MESs of a taxon cannot increase too much in the second half of its life, since it has to go extinct at a fixed point

This leads to a necessary increase of the presence in the first half of a taxon’s life and a necessary decrease of the presence in the second half of a taxon’s life, which are both based on conditioning.

As an example for the effects of conditioning, we give a basic example using a random walk (Y_n)_n∈ℕ (see Foote (2007),foote2014). The random walk is defined by (5)${\displaystyle Y_n=\sum _{i=1}^n{X}_i,}$ where (6)${\displaystyle P\left(X_i=1\right)=p\text{and}P\left(X_i=-1\right)=1-p.}$ Here, p is a value between zero and one that we will call the upward transition probability and X_i is the change in position at the ith step. The random walk defined here is asymmetric in the sense that p can differ from 0.5 and is time-homogeneous in the sense that the transition probabilities do not change through time.

It can be shown that the trajectories of this random walk, conditioned on origination and extinction, are always symmetric after they have been averaged, independent of the choice of p. For the case p = 0.6, this is displayed in Figs. 5 and 6.

This result can be shown as follows: first, extinction can only occur every second time step. Second, it is sufficient to show that the averaged trajectories conditioned on extinction at 2n and origination at 0 are symmetric. This is because symmetry is not affected by x-scaling (see the section “The Way Data are Processed”).

For some fixed 2n, every path from origination at 0 to extinction at 2n consists of 2n steps. Since all feasible trajectories start with presence zero and end with presence zero, and the magnitude of change is the same in each step, half of the 2n steps must increase and half must decrease. Therefore every such trajectory has the probability pⁿ(1 − p)ⁿ. If one of these trajectories is mirror symmetric, it will only contribute symmetrically to the averaged trajectory. If it is not mirror symmetric, it has a mirrored counterpart that has the same probability of occurring. So the asymmetry of this trajectory will in the long run be canceled out by the asymmetry of its mirrored counterpart. This shows that in the long run, all contributions to the averaged trajectory are either symmetric or will be balanced by their counterparts, making the averaged trajectory symmetric itself.

This demonstrates that for the model of an asymmetrical time-homogeneous random walk, conditioning makes it impossible to infer any information about p from the symmetry or asymmetry of the averaged trajectories. In Appendix S4.3, a more general line of argument shows that not only symmetry, but many properties derived from temporal trajectories will systematically change between the unconditioned and the conditioned case, which can lead to erroneous inferences (see Appendix S4.4).

Avoiding noise

In Appendix S3, we suggest a method that avoids averaging, one of the three effects discussed in this paper. This naturally raises the questions whether there are methods that avoid the other effects. The introduction of noise is directly linked to the combination of information from different time intervals and thereby unavoidable when questions about large scale processes are tackled. The chances are that this effect will lessen as both more data and higher resolutions of data are available in the future. That reduction of noise has an impact on results was demonstrated by Liow & Stenseth (2007), where a symmetric averaged trajectory becomes asymmetric as data with lower quality is removed.

The generality of averaging

The effects of averaging can be strong and are not dependent on any model assumption. Especially two points contribute to this.

First, averaging condenses information about probability distributions. In descriptive statistics, this is a desired property, but information about important properties of the distribution such as skewness and the median are inevitably lost.

Second, averaging is deeply linked to most of the limit theorems in probability theory. What unifies all these limit theorems is the convergence of an average over some objects to an invariant object, which is not dependent on the properties of the objects that were averaged in the first place. Some famous examples of such invariant objects are the standard normal distribution (central limit theorem (Klenke, 2008, p. 319)(Lindeberg, 1922; Laplace, 1812)), Brownian motion (Donsker’s invariance principle, see the section “Symmetry and Similarity by Averaging” and (Klenke, 2008, p. 474)) or zero (weak and strong laws of large numbers (Klenke, 2008, p.108-119)). This tendency to converge to an invariant object will automatically make averages appear more similar and reduce their distinguishability.

The properties of averaging can also be passed down. This is for example displayed by the center of gravity (CG) defined by Gould, Gilinsky & German (1987). On average, it displays values that indicate symmetry. This is simply a reflection of the fact that the underlying averaged trajectories are symmetrical as a result of the points discussed in the section “Three Ways to Increase Symmetry and Similarity”.

Conditioning of models

Conditioning formalizes the incorporation of additional knowledge into a model and is therefore common in paleontology, where potentially extinct taxa are examined. It alters the behaviour of a model by fundamentally changing the probability distribution that determine the probability of the temporal trajectories to be observed (Appendices S4.2 and S4.3). One example for this was given by Budd & Mann (2018), who showed that the analysis of extant clades favors clades with an early boost of diversification. This displays that in contrast to averaging (which is commonly part of a statistical procedure) conditioning is part of the question asked, making its effects independent of the statistical approach used.

By mathematical necessity, functions extracting information from trajectories will have a different expectation value for the conditioned and the unconditioned model (Appendices S4.3). Commonly a set of models is given and the task is to identify the model at hand based on the observed trajectories. The result above implies that trying to identify the unconditioned models based on conditioned trajectories will lead to misidentifications (Appendices S4.4, also containing an example).

For these reasons there will be no solution to avoid the effects of conditioning such as better data or new statistical methods. Conditioning completely changes the model used, and statistical methods must be adapted to this new model to avoid misinterpretations of paleontological data. This is especially important when both conditioned and unconditioned models are used, for example when geological data is used to predict developments in the near future.