6 Exploring the congruence class

In theory there are infinitely many speciation and extinction rate functions that are within a congruence class, i.e., all rate functions that have exactly the same likelihood. It is clearly impossible to enumerate all rate functions. Instead of trying to explore literally all rate functions, we will sample from these rate functions instead. This idea is very similar to all sampling based methods. For example, a posterior distribution with a 95% credible interval between 0.01 and 0.02 contains infinitely many values, but we are we comfortable today with simply sampling values proportional to their posterior probability and looking/plotting the distribution.

6.1 Assumptions on rates

There are several plausible assumptions that can be made about diversification rates and how these rates vary. First, you could make assumptions about the magnitude of rate variation, and specifically what are plausible minimum and maximum rates. Clearly, speciation and extinction rates should not have been smaller than 0, so that could be a safe assumption.

Depending on your study system, you could also assume that rates were never larger than, say, 10 events per lineage per million years. For example, if you would have a speciation rate of 10 per lineage per million years, then you would expect every lineage to speciate on average every 0.1 million years. If there would be no extinction, then this would leave the group with just above 22,000 species after only one million years. So perhaps a rate of 10 could really be seen as an upper bound in systems such as vertebrates and plants. However, you could argue that there is extinction, and we should better restrict the net-diversification rates to be at most 10. We get into some more of these restrictions in another section.

The second assumption is about how rates vary over time. One could argue that if diversification rates were low, say 42 million years ago, then you would assume that the diversification rates were also low 43 million years ago. This would mean that diversification rates are autocorrelated and do not vary completely arbitrarily over time.

In the next subsections, we explore a few different options of rate variation over time in the next sections, which include:

• HSMRF (autocorrelated) rates
• Exponentially decreasing rates with stochastic variation
• Linearly increasing rates with stochastic variation
• Episodic constant models (with \(n\) episodes)
• Uncorrelated (independent) rates

With the built-in functions sample.basic.models and sample.rates, users have freedom to explore these, and many other forms of rate variation. If these functions are not sufficient, as shown in previous sections, it is easy to customize the sampling to your liking.

6.2 Sample alternative extinction rates

6.2.1 Temporally autocorrelated

Here, we generate a variety of alternative extinction rates using HSMRF models (Magee et al., 2020). When sampling from models, we get random walks, or temporally autocorrelated, extinction rates, and the HSMRF tends to produce random walks that look like piecewise constant functions. The rates produced here have no overall trend, the median change from past to present is 0. We use the function sample.basic.models to construct these rates, which by default assumes that rates are between 0 and 10, though here we change these boundaries slightly to bound the rates to be between 0 and 1. Our motivation for these values was, that (i) having a maximal expected extinction time per lineage of 1.0 million years seems reasonable from the fossil record, and (ii) we estimated a maximal speciation rate of 0.4, which is clearly smaller than 1.0.

times <- seq(0, max(my_model$time), length.out = 100)
extinction_rate_samples <- function() {
   sample.basic.models( times = times,
                        model="MRF",
                        max.rate=1)
}

samples <- sample.congruence.class(
  my_model,
  num.samples=10,
  rate.type="extinction",
  sample.extinction.rates=extinction_rate_samples)

We can plot the samples using

plot(samples)

Figure 6.1: Rate functions assuming HSMRF-distributed extinction rates

p <- summarize.trends(samples,
                      threshold = 0.02)
plot(p)

Figure 6.2: Congruence class assuming linearly changing extinction rates.

6.2.2 Exponentially decreasing

Here, we generate a variety of alternative extinction rates using HSMRF models (Magee et al., 2020) with an exponential trend. We otherwise make the same assumptions as above. To help the rates we simulate fit within 0 and 1, we set the fold change average to be 2 instead of the default 3.

times <- seq(0, max(my_model$time), length.out = 100)
extinction_rate_samples <- function() {
   sample.basic.models( times = times,
                        model="exponential",
                        max.rate=1,
                        fc.mean=2)
}

samples <- sample.congruence.class(
  my_model,
  num.samples=10,
  rate.type="extinction",
  sample.extinction.rates=extinction_rate_samples)

We can plot the samples using

plot(samples)

Figure 6.3: Rate functions assuming exponential trends in the extinction rate

p <- summarize.trends(samples, threshold = 0.02)
plot(p)

Figure 6.4: Congruence class assuming exponentially changing extinction rates

6.2.3 Linearly increasing

Here, we generate a variety of alternative extinction rates using GMRF models (Magee et al., 2020) with a linear trend. Unlike the exponential section, we here assume that rates have increased, that the increase has been monotonic (there were no decreases, where in the exponential decreases above there were occasional increases), and that the stochastic variability comes from a GMRF. We otherwise make the same assumptions as above.

times <- seq(0, max(my_model$time), length.out = 100)
extinction_rate_samples <- function() {
   sample.basic.models( times = times,
                        direction="increase",
                        model="linear",
                        monotonic=TRUE,
                        MRF.type="GMRF",
                        max.rate=1,
                        fc.mean=2)
}

samples <- sample.congruence.class(
  my_model,
  num.samples=10,
  rate.type="extinction",
  sample.extinction.rates=extinction_rate_samples)

We can plot the samples using

plot(samples)

Figure 6.5: Congruence class assuming linearly changing extinction rates.

p <- summarize.trends(samples, threshold = 0.02)
plot(p)

Figure 6.6: Congruence class assuming linearly changing extinction rates.

6.2.4 Episodic

Here, we generate a variety of alternative extinction rates without stochastic noise. Without noise (noisy=FALSE), all changes will be in the same direction, here decreases. The model we generate is an episodic model with 5 rates. Here, we consider decreasing rates with larger changes from past to present than the previous section. To help rates stay within the [0,1] bounds we have chosen, we set the rate median at present to be slightly smaller than the default 0.1.

times <- seq(0, max(my_model$time), length.out = 100)
extinction_rate_samples <- function() {
   sample.basic.models( times = times,
                        direction="decrease",
                        model="episodic5",
                        noisy=FALSE,
                        max.rate=1,
                        fc.mean=4,
                        rate0.median=0.05)
}

samples <- sample.congruence.class(
  my_model,
  num.samples=10,
  rate.type="extinction",
  sample.extinction.rates=extinction_rate_samples)

We can plot the samples using

plot(samples)

p <- summarize.trends(samples, threshold = 0.02)
plot(p)

Figure 6.7: Congruence class assuming episodicly changing extinction rates.

6.2.5 Uniform independent

The previous functions used sample.basic.models to construct a variety of pre-defined hypotheses. ACDC also implements sample.rates to allow more flexibility in user-defined rate functions (and of course the user is free to define their own functions as well). Here we show this with IID rates. We assume that the extinction rates could have had any value between 0 and 1 (you could easily pick different numbers here). Our motivation for these values was, that (i) having a maximal expected extinction time per lineage of 1.0 million years seems reasonable from the fossil record, and (ii) we estimated a maximal speciation rate of 0.4, which is clearly smaller than 1.0.

We furthermore assume that each extinction rate is independent of the previous time interval. That means that we can model any rate function, even completely crazy zig-zagging functions.

coarse_times <- seq(0, max(my_model$time), length.out = 10)
rsample_extinction <- function(n) runif(n,0,1)
extinction_rate_samples <- function() {
   sample.rates( times = coarse_times,
                 rsample=rsample_extinction,
                 rsample0=NULL,
                 autocorrelated=FALSE)
}

samples <- sample.congruence.class(
  my_model,
  num.samples=10,
  rate.type="extinction",
  sample.extinction.rates=extinction_rate_samples)

We can plot the samples to inspect them visually

plot(samples)

Figure 6.8: Rate functions for the congruence class with proposed alternative extinction rates whose values are drawn from an independent uniform distribution for each episode.

Next we can plot the directional trends

p <- summarize.trends(samples, threshold = 0.02)
plot(p)

Figure 6.9: Congruence class assuming independent and uniformly distributed extinction rates.

6.3 Sample alternative speciation rates

It is also possible to sample alternative speciation rates. Since the congruence class requires a shared \(\lambda_0\) for all models in the congruence class, we must constrain our sampler to always start on the same speciation rate at the present.

6.3.1 Temporally autocorrelated

Here we choose to sample speciation rates according to a horseshoe Markov random field distribution.

times <- seq(0, max(my_model$time), length.out = 100)
speciation_rate_samples <- function() {
   sample.basic.models( times = times,
                        rate0 = my_model$lambda(0.0),
                        model = "MRF",
                        MRF.type = "HSMRF",
                        max.rate = 1.0)
}

Next, we can sample ten models from the congruence class.

samples <- sample.congruence.class(
  my_model,
  num.samples=10,
  rate.type="speciation",
  sample.speciation.rates = speciation_rate_samples)

We can plot the samples using

plot(samples)

Figure 6.10: Rate functions for the congruence class. The reference model is given in black. The alternative speciation rates are proposed by sampling from a temporally autocorrelated (HSMRF) model.

Notice that all of our congruent models have a negative extinction rate near the present. This becomes more evident if we plot the extinction rates on its own panel and adjust the axis limits.

p <- plot(samples)[[2]] + 
  coord_cartesian(ylim = c(-1, 1)) +
  xlab("time before present") +
  ylab("rate")
plot(p)

Figure 6.11: Inferred congruent extinction rates, zoomed in for clarity.

The models depicted in 6.11 do not make sense, since rates are only defined on the positive numbers. Similarly to before, we can evaluate the directional changes in the extinction rate

p <- summarize.trends(samples, threshold = 0.02, rate_name = "mu")
plot(p)

Figure 6.12: A summary of directional trends through time across the congruent models.

Notice how the reference model (top row) is flat throughout the entire time span. Yet, when we propose alternative speciation rates, the inferred (extinction) rates depict a pattern of two decreases and one increase in the near present. We can also inspect the directional trends for the net-diversification rate (\(\delta(t) = \lambda(t) - \mu(t)\))

p <- summarize.trends(samples, threshold = 0.02, rate_name = "delta")
plot(p)

Figure 6.13: Directional trends in the net-diversification rate. The congruence class was constructed using alternative speciation rates, and all congruent models have negative extinction rates.