# In richelbilderbeek/razzo: The Error if Nature is MBD

df$n_mb_taxa <- df$f_mb * df$n_taxa knitr::kable(head(df))  ## 4. MB regimes: calculate$q$We pick a desired MB regime. There are three MB regimes, in which the co-occuring speciation events are: • few events of strong intensity • some events of intermediate intensity • many events of modest intensity Starting from two lineages and$q = 1$, the number of (extant and extinct) lineages$N$will be $$N = 2^{n_{\nu}}$$ As we know the highest number of taxa$N$we need: max_n <- max(df$n_mb_taxa)
print(max_n)


We can calculate the minimal number of events we need, rounded up to the nearest integer:

$$N = 2^{n_{\nu}}$$

$$n_{\nu} = \left \lceil{\frac{log(N)}{log(2)}}\right \rceil$$

min_n_events <- ceiling(log(max_n) / log(2))
print(min_n_events)


So, to be able to arrive at the highest needed number of taxa, one needs at least r min_n_events events.

We set the number of co-occurring speciation events, $n_{\nu}$, as:

n_nu_eventses <- c(min_n_events, min_n_events * 2, min_n_events * 4) # Gollumese plural


Now we know the amount of taxa created by the MB process $N_{\mathbb{M}}$, the number of nu events, and that we start with two taxa, we can calculate the intensity of an MB event $q$:

$$N_{\mathbb{M}} = 2 (1.0 + q)^{n_{\nu}}$$

Because we start with two taxa, $N_{\mathbb{M}}$ can never become less than 2.

In the current context, we known the number of MB taxa that must be generated $N_{\mathbb{M}}$ and the number of co-occurring speciation events $n_{\nu}$. All we need to do is calculate $q$.

This we can solve this analytically in Maxima:

solve(N  = 2 * ((1 + q) ^ n), q);


results in:

$$q = - \frac{2^{1/n_{\nu}} - N_{\mathbb{M}}^{1/n_{\nu}}}{2^{1/n_{\nu}}}$$

Or plotting it in Maxima:

f(n, N) := rhs(solve(N = 2 * ((1 + q) ^ n), q)[1]);
plot3d(f(n, N), [n, 9, 16], [N, 15, 500], [zlabel, "q"]);


Converting this to an R function:

calc_q <- function(n_nu_events, n_mb_taxa) {
testit::assert(n_nu_events >= 0)
testit::assert(n_mb_taxa >= 0)
# From the equation, zero MB taxa results in strange values for q.
# Just use q == zero instead
if (n_mb_taxa < 2) {
return (0.0)
}

q <- - ((2.0 ^ (1.0 / n_nu_events)) - (n_mb_taxa ^ (1.0 / n_nu_events))) /
(2.0 ^ (1.0 / n_nu_events))
testit::assert(q >= 0.0)
testit::assert(q <= 1.0)
q
}


As we know the number of (extinct and extant) taxa $N$, we can calculate this:

df <- expand.grid(
n_extant_taxa = n_extant_taxas,
f_extinct = f_exinctses,
f_mb = f_mbs,
n_nu_events = n_nu_eventses
)
df$n_taxa <- df$n_extant_taxa / (1.0 - df$f_extinct) df$n_mb_taxa <- df$f_mb * df$n_taxa


Calculating the $q$s:

calc_q(n_nu_events = 9, n_mb_taxa = 15)
df$q <- NA for (i in seq_along(df$q)) {


## 6. speciation rate $\lambda$

The number of (extinct and extant) species $N_{\mathbb{B}}$ created is dependent on the initial number of lineages (which is two), the crown age $T$ and per-lineage speciation rate $\lambda$:

$$N_{\mathbb{B}} = 2 \cdot T^{\lambda}$$

Solving this for per-lineage speciation rate $\lambda$:

$$log(\frac{N_{\mathbb{B}}}{2}) = \lambda \cdot log(T)$$

$$\frac{log(\frac{N_{\mathbb{B}}}{2})}{log(T)} = \lambda$$

$$\lambda = \log_T(\frac{N_{\mathbb{B}}}{2})$$

calc_lambda <- function(n_bd_taxa, crown_age) {
testit::assert(n_bd_taxa >= 0)
testit::assert(crown_age >= 0)
# From the equation, zero BD taxa results in strange values for lambda.
# Just use lambda == zero instead
if (n_bd_taxa < 2) {
return (0.0)
}
lambda <- log(n_bd_taxa / 2, base = crown_age)
testit::assert(lambda >= 0.0)
lambda
}


Calculate the $\lambda$s:

df$lambda <- NA for (i in seq_along(df$lambda)) {
df$lambda[i] <- calc_lambda(n_bd_taxa = df$n_bd_taxa[i], crown_age = crown_age)
}


## 7. number of extinct taxa

The number of taxa gone extinct $N_{\dagger}$ has a simple ralationship to the total number of (extinct and extant) taxa $N$ and the extant taxa $N_{\star}$:

$$N = N_{\star} + N_{\dagger}$$

$$N_{\dagger} = N - N_{\star}$$

So the number of extinct taxa $N_{\dagger}$ is:

for (i in seq_along(df$nu)) { df$nu[i] <- df$n_nu_events[i] / crown_age } knitr::kable( head( dplyr::select( df, c("n_extant_taxa", "n_mb_taxa", "n_nu_events", "q", "nu") ) ) )  ## Conclusion There you have it, all simulation parameters! From • a desired number of extant taxa ($n$) • a desired fraction of taxa that goes extinct ($f_{\dagger}$) • a desired fraction of taxa formed by the MB process ($f_{\mathbb{M}}$) • a desired MB regime ($n_{\nu}$) we derive the MBD parameters • speciation rate$\lambda$• extinction rate$\mu$• rate at which a co-occuring speciation event in triggered$\nu$• fraction of species that speciate during such a co-occurring speciation event$q$knitr::kable( dplyr::select( df, c( "n_extant_taxa", "f_extinct", "f_mb", "n_nu_events", "lambda", "mu", "nu", "q" ) ) )  ## Data verification testit::assert(all(df$n_taxa == df$n_extant_taxa + df$n_extinct_taxa))


richelbilderbeek/razzo documentation built on Feb. 18, 2020, 7:24 p.m.