---
title: "CRABS: Congruent Rate Analyses in Birth-death Scenarios"
date: '`r Sys.Date()`'
author: 
  - Sebastian Höhna
  - Bjørn T. Kopperud
  - Andrew F. Magee
  - Jérémy Andréoletti
output:
  rmarkdown::html_vignette: 
    toc_depth: 1
  bookdown::html_document2: 
    fig_caption: yes
    toc_depth: 2
    toc_float: true
  bookdown::pdf_document2:
    fig_caption: yes
    keep_tex: yes
    toc: yes
    toc_depth: 2
    includes:
      in_header: preamble.tex
  bookdown::gitbook:
    lib_dir: lib
    fig_caption: yes
    toc_depth: 3
fontsize: 12pt
documentclass: article
classoption:
vignette: |
  %\VignetteKeyword{Documentation}
  %\VignetteIndexEntry{Congruence Classes}
  %\VignettePackage{CRABS}
  %\VignetteEngine{knitr::rmarkdown}
bibliography: literature.bib
csl: methods-in-ecology-and-evolution.csl
pandoc_args:
  - --filter=pandoc-eqnos
---


<!-- % revise margins -->
<!-- \setlength{\headheight}{0.0in} -->
<!-- \setlength{\topmargin}{0.0in} -->
<!-- \setlength{\headsep}{0.0in} -->
<!-- \setlength{\textheight}{8.65in} -->
<!-- \setlength{\footskip}{0.35in} -->
<!-- \setlength{\oddsidemargin}{0.0in} -->
<!-- \setlength{\evensidemargin}{0.0in} -->
<!-- \setlength{\textwidth}{8.5in} -->
<!-- \setlength{\parskip}{6pt} -->
<!-- \setlength{\parindent}{0pt} -->


```{r knitr-options, include=FALSE}
library(knitr)
opts_chunk$set(
  fig.pos = "H",
  fig.align='center',
  concordance=TRUE,
  fig.width=4, 
  fig.height=4, 
  fig.path='figures/',
  warning=FALSE,
  message=FALSE,
  tidy=FALSE
)

options(width=60)
set.seed(53079239)
# install R/CRABS package if necessary:
#if(!require("CRABS")) install.packages("CRABS", repos="https://cran.us.r-project.org")
```

# Overview

This is a vignette, or long-form documentation, accompanying the R-package CRABS: Congruent Rate Analyses in Birth-death Scenarios.

The main purpose of `CRABS` is to explore congruence classes of macroevolutionary diversification rates.
Specifically, this tutorial has two main aims/applications:

* For a given speciation and extinction rate, you can compute the alternative speciation and extinction rates via the pulled net-diversification rate 
* For a given speciation and extinction rate, which might have been estimated from a phylogenetic tree or provided simply for exploration purposes, you can explore the entire congruence class given some rate assumptions (e.g., type of change or min and max rates).

The first use case corresponds to testing specific hypotheses while the second use case corresponds to exploring the congruence class for general patterns/trends.
We will explain these in detail below with an empirical example.



# Installation

To get started with exploring congruence classes, you need to have our package `CRABS` installed.
First, make sure that you have a recent version of R installed.
You can install the stable version of `CRABS` from CRAN:
```{r inst-cran, eval=FALSE}
install.packages("CRABS")
```

We host the developmental version of `CRABS` on github. If you don't already have it, the `devtools` R-package makes it easy to install directly from github. Using devtools makes it easy to install `CRABS` directly from GitHub.
```{r load-devtools, eval=FALSE}
install.packages("devtools")
library(devtools)
install_github("afmagee/CRABS")
```
<!-- % Sebastian: This is only for now until we have the package properly on GitHub -->
```{r local-install-CRABS, echo=FALSE}
install.packages("..",type="source",repos=NULL)
```


You can then load the `CRABS` package using the `library` function:
```{r load-CRABS}
library(CRABS)
```
Now you are ready to use `CRABS` to explore congruence classes in macroevolutionary diversification rates.

<!-- Dependencies for the vignette -->
```{r vignette-dependencies}
library(dplyr)
```



# Features

##  Overview of functions in `CRABS`

Here we give you an overview of the main functions available in `CRABS`.


| Function      | Description  |
| :---        |    :----    |
| `create.model`      | Creates a diversification model object, which consists of the so-called reference speciation and extinction rates.        |
| `congruent.models` | Generates alternative models that are congruent with a provided reference model(s). The user needs to explicitly supply a list of alternative extinction and/or speciation rate functions. |
| `sample.basic.model`   |  Samples a rate function under specified assumptions. More information in section \@ref(samplebasic)   |
| `sample.basic.models.joint` | Jointly samples speciation and extinction trajectories. More information in section \@ref(samplejoint) |
| `sample.rates`  | Helps you to randomly sample either a new speciation or extinction rate function given some constraints and/or expected behavior (e.g., autocorrelated vs independent rates) |
| `sample.congruence.class` | Randomly draws a set of samples from the congruence class. The user provides a function that draws extinction rates, and/or a function that draws speciation rates.  | 
| `sample.congruence.class.posterior` | Uses `sample.basic.models` to samples a set of congruent models for each sample in the posterior distribution | 
| `summarize.trends` |  Calculates the slope in the rate function, and plots a summary of whether the rate was increasing, decreasing or flat for each episode |
| `summarize.posterior` | Similar to `summarize.trends`, but instead iterating over all congruent samples for each posterior sample. This plot is sorted |
| `full.plot.regularity.thresholds` | Plots the rate functions after filtering them to only retain the most regular (smooth) trajectories. |


You can consult the R documentation for each of these by placing a question mark in front of the name. For example, try calling `?create.model` or `?sample.basic.models` or `?sample.basic.models.joint` in R. 
We use these functions throughout the sections in the vignette.
For `sample.basic.models` and `sample.basic.models.joint`, we have separate sections explaining the most important flags and use cases because they might represent more complex but important functions of `CRABS`.

## `sample.basic.models` {#samplebasic}
This function is designed to enable easy generation of commonly employed diversification-rate scenarios, for either speciation or extinction rates. For sampling these rates jointly, see the following section \@ref(samplejoint).

The basic models are specified by the argument `model`: exponentially increasing/decreasing rates, linearly decreasing rates, episodic (or piecewise constant rates), and Markov random field (MRF) models (constant rates are also possible).
The direction of change (increase or decrease) is controlled by the argument `direction`.
The function will fit a piecewise linear rate function, using the argument `times` to represent a vector for the discrete time steps. 
The number of episodes is equal to the length of the `times` vector.
So, for example, the code `sample.basic.models()`, which uses all defaults, samples exponentially decreasing rates.

By default, the exponential, linear, and episodic models are sampled with an additional random component given by a MRF.
This can be stopped by setting `noisy=FALSE`.
The argument `MRF.type` allows you to control the type of Markov random field model, the default is "HSMRF" for Horseshoe Markov random fields [@Magee2020], which produces jumpier changes.
The alternative, "GMRF" for Gaussian Markov random field, produces smoother changes.
This MRF noise allows rates to go up or down, such that `direction` gives the overall trend across many sampled rates rather than a guarantee for any one rate.
To change this, and force rates to go one direction while allowing stochastic noise from an MRF, set `monotonic=TRUE`.

To have no trends, and pure MRF models, specify `model="MRF"`.
Note that this means that setting `model="MRF"` and `noisy=FALSE` amounts to a constant-rate model (by disabling both trends and stochastic variability).

The argument `fc.mean` dictates (more or less) the average amount of change from first to last rate as a fold change.
Specifically, `fc.mean` must be a value greater than 1 (1.00...01 is allowed), and `direction` determines the direction of change.
If `direction="increasing"`, then `fc.mean` defines first/last, while if `direction="decreasing"`, `fc.mean` defines last/first.
The default is 3, which is biologically plausible but not the universally best option.
Fold changes are drawn from a Gamma distribution with rate 1.25, offset 1, and a shape parameter that gives the desired mean.
This means, that every drawn single rate function has it's own fold change.

All rate curves are sampled by rejection sampling to guarantee that the rates are all between `min.rate` and `max.rate`.
By default, rates are assumed to be in [0,10] which on a million-year time scale is probably the extreme end of plausible variation.
This rejection sampling means that the actual average change may not quite be what it is set to be.
For example if one sets the model to be a decrease through time, sets the rate at present to be 0.5, and the maximum rate to be 0.75, clearly the average fold change could not be a larger than 1.5x increase.
Be aware, setting rejection limits that make it difficult to sample trajectories can make this function (which is normally quite fast) take a very long time!

The initial rate is drawn from a lognormal distribution with median `rate0.median` and log(sigma) = `rate0.logsd`, which is to say a Lognormal(log(`rate0.median`),`rate0.logsd`) distribution.
The default is rates with median 0.1 and a 95% range covering [0.01,1].
Note that the argument `rate0` overrides drawing random starting values if it is set.
This allows sampling speciation rates, as the initial speciation rate must be the same in all replicates.


## `sample.basic.models.joint` {#samplejoint}
This function is designed to jointly generate speciation and extinction rates that are congruent with the original diversification scenario, while inherently preventing negative extinction rates.

`sample.basic.models.joint` primarily relies on generating a trajectory for a parameter $p$, which serves as a reference for balancing the speciation and extinction rates. When $p=0$, the speciation rate is constant (barring negative extinction rates) and the time-varying pattern is wholly reflected in the extinction rate. Conversely, when $p=1$, the extinction rate is constant, and the time-varying pattern is reflected in the speciation rate. Each unique $p$-trajectory is tied to a distinct pair of speciation and extinction rate trajectories, given a starting extinction rate at present.

Like in the `sample.basic.models` function, the `MRF.type` argument ("HSMRF" or "GMRF") controls the type of Markov random field used for adding stochastic noise to the trajectories, although here it corresponds to the $p$-trajectory rather than the speciation or extinction trajectories themselves. The scale of the stochastic noise can be adjusted using `mrf.sd.scale`, which defaults to 1.0 and can be set to 0.0 for constant $p$-trajectories.

The congruence class in which to sample is specified through the pulled diversification rate function (`p.delta`), together with the speciation rate at the present. The extinction rate at present can either be specified or drawn randomly from a lognormal distribution, similarly to `sample.basic.models`.

Bounds for the rates can also be imposed using `min.lambda`, `max.lambda`, `min.mu`, and `max.mu` (which is usually unnecessary with this approach). Similarly, `min.p` and `max.p` can be used to set the lower and upper bounds of the p trajectory, allowing for finer control over the balance between speciation and extinction rates in the generated scenarios.


# Creating the congruence class

Here, we will start with some theoretical background (although quite simplified) on congruence class. Our main aim is to motivate, guide and explain our choices for how to explore the congruence class.


## Rate transformations

One important aspect of the congruence class is that if we know the pulled rates, i.e., the pulled speciation or net-diversification rates, then we only need to know additionally either the new speciation or extinction rate to solve for the missing rate function.
That means, we can compute the pulled diversification rate for our reference model, i.e., using our empirically estimated speciation and extinction rates from the previous sections, and then specify whichever alternative speciation or extinction rate function.
Then we solve the pulled diversification rate equation for this new rate function.

The pulled diversification rate is defined as
\begin{equation}
	r_p(t) = \lambda(t) - \mu(t) + \frac{1}{\lambda(t)} \frac{\text{d}\lambda(t)}{\text{d}t}
	(\#eq:pulled-diversification)
\end{equation}

As you may notice, the latter part of this equation entails the derivative of the speciation rate function.
This derivative is only known for simple rate models, for example the constant rate process, linear rates or exponential rates. 
For more complex models, it is often impractical or impossible to find the derivative analytically.
In order to solve (eq. \@ref(pulled-diversification)) for $mu(t)$, the rearrangement of $\mu(t) = \lambda(t) - r_p(t) + \frac{1}{\lambda(t)}\frac{d\lambda}{dt}$ is trivial, and we use the finite difference method to approximate the derivative $\frac{d\lambda}{dt}$ evaluated at time $t$.
In order to get $\lambda(t)$, we use an ODE (Ordinary Differential Equation) solver to compute the solution numerically. More specifically, we use a Runge-Kutta method of order five, implemented in Fortran. We use the function `radau()` in the R-package `deSolve` as an interface to access the ODE solver in R.

<!-- Computationally, it is difficult if not impossible to work with arbitrary continuous rate functions. -->
<!-- Here we follow a common approach to use piecewise-linear rates as an approximation for an arbitrary continuous rate function. -->
<!-- You could create as many intervals as you want, and in our examples we use 1000. -->

<!-- Since we are using piecewise linear rates, we get (at the interval times) the pulled diversification rate function as -->
<!-- \begin{equation} -->
<!-- 	r_p(t_i) = \lambda_i - \mu_i + \frac{1}{\lambda_i} \frac{\lambda_i - \lambda_{i-1}}{\Delta t} -->
<!-- 	(\#eq:pulled-diversification-linear) -->
<!-- \end{equation} -->
<!-- This function is much simpler to compute and is used in our package. -->

<!-- Next, if we know the pulled diversification rate and the new extinction rate, then we can compute the speciation rate as -->
<!-- \begin{equation} -->
<!-- 	\lambda_i = \frac{\sqrt{4\times \lambda_{i-1}\times \Delta t + (r_p(t_i) \times \Delta t+\mu_i \times \Delta t - 1)^2} + r_p(t_i)\times \Delta t + \mu_i \times \Delta t - 1}{2\times \Delta t} -->
<!-- 	(\#eq:lambdasolve) -->
<!-- \end{equation} -->
<!-- Equation \@ref(eq:lambdasolve) was obtained by solving the pulled diversification rate \@ref(eq:pulled-diversification-linear) for $\lambda_i$.  -->
<!-- Conversely, if we know the speciation rate and the pulled diversification rates, then we can solve Equation \@ref(eq:pulled-diversification-linear) for $\mu$ and get -->
<!-- \begin{equation} -->
<!-- 	\mu_i = \lambda_i - r_p(t_i) + \frac{1}{\lambda_i} \frac{\lambda_i - \lambda_{i-1}}{\Delta t} -->
<!-- 	(\#eq:musolve) -->
<!-- \end{equation} -->

These are our fundamental pieces about the congruence classes.
We have the reference model, from which we compute the pulled diversification rate.
Then you can either specify alternative speciation or extinction rate functions, for which we then compute the matching extinction or speciation rate function respectively, or you can explore a space of speciation or extinction rate functions.
We will walk through examples of these in the following sections.


## Slopes and $\Delta t$

Given a congruence class, you may want to either visualize or summarise the models that belong to the congruence class.
Visualisation of the alternative models is done in a hopefully intuitive way so we explain it later.
Summarizing alternative models is somewhat less established.
In our view, the main objective of these diversification functions is to estimate if speciation and/or extinction rates have changed over time.
A change in rate over time is mathematically represented by the slope of the rate function; a large slope means more/steeper change in less time, while a small value for the slope corresponds to rather constant rates.
Therefore, we summarize the models by computing the slope of the rate function.
Then, we assess each time interval independently if the slope was larger than a specified threshold.
The sign of the slope shows the direction of change, either increasing or decreasing.

We will come back to the summaries later when we created samples from the congruence class.



## The model object

In order to create a congruence class, we require some information about the extinction and speciation rate functions in at least one of the models in the congruence class. 
In the following examples we will use the results from an episodic birth-death model fitted in RevBayes (see Section \@ref(revbayes)), but we note that `CRABS` can work with any rate estimates. 
To simplify this section, we provide the results in a convenient format:
```{r load-revbayes}
data("primates_ebd")
head(primates_ebd)
```

We also provide an example for fitting birth-death models in TESS (Section \@ref(tess)) and TreePar (Section \@ref(treepar)). 
If you wish to use these results instead, simply call `data("primates_ebd_treepar")` or `data("primates_ebd_tess")`. 
After we have estimated our birth-death rates, we transform the data into piecewise linear rate functions in R
```{r revbayes-piecewise-rates}
times <- primates_ebd[["time"]]

lambda <- approxfun(times, primates_ebd[["lambda"]])
mu <- approxfun(times, primates_ebd[["mu"]])
```

Next, we create a congruence class object which contains all the transformed rate functions (i.e., the net-diversification rate and relative extinction rate) as well as the pulled rates (i.e., the pulled net-diversification and the pulled speciation rate).
```{r fig.cap="The diversification rates in the reference model.", fig.width = 6, fig.height = 6}
max_t <- max(primates_ebd[["time"]])
times_fine <- seq(0, max_t, length.out = 1000)

my_model <- create.model(lambda, 
                         mu, 
                         times_fine)
plot(my_model)
```
The `my_model` includes all the rate functions.
You could plot these yourself, or you can use our simple plotting function.

Now that we have our reference model and congruence class functions, we can explore alternative rate functions within the congruence class.



# Testing alternative hypotheses

In this section we focus on specific alternative hypotheses for diversification rates. 
The first part of our congruence class exploration focuses on specific alternative models.
That is, we first start by providing a set of alternative extinction rate functions and compute the matching speciation rate functions.
Then, we turn the approach around and specify some alternative speciation rate functions for which we compute the matching extinction rate functions.

## Testing alternative extinction rates {#alt-ext}

In this subsection we focus entirely on alternative hypotheses of extinction rates for which we compute the resulting speciation rates.


### Constant extinction rates

Let us assume that the extinction rate was constant, but at a range of different values.
Specifically, let's assume that the true extinction rate was $\mu(t)=\{0.0,0.025,0.05,0.075,0.1\}$.
For our models we need to create extinction rate functions with these values.
```{r }
mu_vals <- seq(0,0.1,0.025)
mu1 <- list()
for (i in 1:length(mu_vals)) {
  mu1[[i]] = local({
    mu = mu_vals[i]
    function(t) { rep(mu,length(t)) }
  })
}
```
Now we compute the new speciation rate given the extinction rates.
```{r eval=TRUE}
alt_models <- congruent.models(my_model, mus = mu1)
```
Let's have a look at these *alternative models*.
```{r congruent-models-mu-constant, fig.cap = "Diversification rate functions within the congruence class assuming different constant extinction rate functions. The dashed lines show the different alternatives that we explored here, and the solid lines represent our reference model. Here we notice that the overall pattern of speciation rates as well as net-diversification rates is fairly robust to alternative constant extinction rates within the congruence class. Only the relative extinction rate changes is value, which is not surprising since we explored different extinction rates.", fig.height = 4, fig.width = 4}
plot( alt_models )
```
Next, we can inspect the summary of the trends in the speciation rate through time
```{r summary-mu-constant, fig.cap = "Summary of directional trends in the congruence class.", fig.height = 3, fig.width = 6}
p <- summarize.trends(alt_models, threshold = 0.02)
plot(p)
```
Notice that there is variation in the slope of the speciation rate ($\Delta\lambda$) throughout the entire time span of the phylogenetic tree. Yet, when we only consider rate changes that are smaller or greater than our threshold of $\epsilon = 0.02 \Delta\lambda$ per million years, it becomes evident that the constructed congruent models are in complete agreement of the two increases in speciation rate, as well as the more recent decline in speciation rate.


### Linearly increasing extinction rates

Next, we will explore different linear rate functions for the extinction rate.
Let us assume that the extinction rate was linear with a range of different slopes.
Specifically, let's assume that the slope was $\alpha=\{0.0,0.0005, \ldots, 0.002\}$.
We played around with these $\alpha$ values to get reasonable rate functions.
You will need to modify these for your dataset.
Note that in `CRABS`, rate functions are defined in backwards time, so while these look like decreasing extinction rates, they are in fact increasing extinction rates when you consider time forward from past to present.

We create extinction rate functions with these $\alpha$ values.
```{r }
mu_min = 0.01
mu_vals <- seq(0,0.002,0.0005)
mu1 <- list()
for (i in 1:length(mu_vals)) {
  mu1[[i]] = local({
    mu_slope = mu_vals[i]
    function(t) { mu_min + mu_slope * (max_t-t) }
  })
}
```
Now we compute the new speciation rate given the extinction rates.
```{r eval=TRUE}
alt_models <- congruent.models(my_model, mus = mu1)
```
Let's have a look at these *alternative models*.
```{r congruent-models-mu-linear-inc, fig.cap = "Diversification rate function for different linear extinction rates within the same congruence class. The dashed lines show the different alternatives that we explored here, and the solid lines represent our reference model. We observe that the speciation and net-diversification rates are robust to different linear extinction rates within the congruence. As expected, the relative extinction rate changes because we chose different extinction rates.", fig.height = 4, fig.width = 4}
plot(alt_models)
```
Next, we can inspect the summary of the trends in the speciation rate through time
```{r summary-mu-linear, fig.cap = "Summary of directional trends in the congruence class.", fig.height = 3, fig.width = 6}
p <- summarize.trends(alt_models, threshold = 0.02)
plot(p)
```
Similar to the example with constant alternative extinctionrate models, there is some variation in the trends, but the same conclusion remains. Namely, for our threshold of $\epsilon = 0.02 \Delta\lambda$ per million years, it is evident that the constructed congruent models are in complete agreement of the two increases in speciation rate, as well as the more recent decline in speciation rate.

### Exponential increasing extinction rates

Next, similar to the linear rate functions we explore some exponential rate functions.
We start here with the exponentially increasing rates.
Recall that the extinction rate is a function of time such that $\mu(0)$ is the rate at the present and $\mu(1)$ is the rate 1 time unit previously, thus what appear to be exponential increases are in fact decreases.
Specifically, let's assume that the growth rate was $\alpha=\{0.0,0.005, \ldots, 0.02\}$.
Then, we define the extinction rate function as
```{r }
mu_max = 0.015
mu_vals <- seq(0,0.02,0.005)
mu1 <- list()
for (i in 1:length(mu_vals)) {
  mu1[[i]] = local({
    mu_factor = mu_vals[i]
    function(t) { mu_max * exp(mu_factor*(max_t-t)) }
  })
}
```
Now we compute the new speciation rate given the extinction rates.
```{r eval=TRUE}
alt_models <- congruent.models( my_model, mus = mu1 )
```
Let's have a look at these *alternative models*.
```{r congruent-models-mu-exp-inc, fig.cap = "Diversification rate functions for exponentially increasing extinction rates within the congruence class. The dashed lines show the different alternatives that we explored here, and the solid lines represent our reference model. As in our previous example, the speciation rate and net-diversification rate are robust to the difference choices of extinction rates", fig.height = 4, fig.width = 4}
plot( alt_models )
```

Next, we can inspect the summary of the trends in the speciation rate through time
```{r summary-mu-exponential, fig.cap = "Summary of directional trends in the congruence class.", fig.height = 3, fig.width = 6}
p <- summarize.trends(alt_models, threshold = 0.02)
plot(p)
```
Similar to the two previous examples, there is some variation in the trends, but the same conclusion remains. Namely, for our threshold of $\epsilon = 0.02 \Delta\lambda$ per million years, it is evident that the constructed congruent models are in complete agreement of the two increases in speciation rate, as well as the more recent decline in speciation rate.


### Two-epoch extinction rates

Let us assume that the extinction rate changed abruptly at a pre-defined change-point.
We'll try to fix the change-point to the estimated rapid increase in diversification.
Thus, we are modeling the extinction rate as a 2-epoch model.
Furthermore, let us try a range of different factors of increase/decrease.
Specifically, let's assume that the rate increase/decrease was by a factor of $\alpha=\{1,2,4\}$.
Again, we create the list of extinction rate functions.
```{r }
mu_min = 0.05
t_shift = 5.5
alpha_vals <- seq(-2,2,1)
mu1 <- list()
for (i in 1:length(alpha_vals)) {
  mu1[[i]] = local({
    alpha = alpha_vals[i]
    function(t) { ifelse( t < t_shift, mu_min*2^alpha, mu_min) }
  })
}
```
Now we compute the new speciation rate given the extinction rates.
```{r eval=TRUE}
alt_models <- congruent.models( my_model, mus = mu1 )
```
Let's have a look at these *alternative models*.
```{r congruent-models-mu-shift, fig.cap = "Diversification rate functions for 2-epoch extinction rate functions within the congruence class. The solid line shows the reference model and the dashed lines show the alternative models. We observe no noticeable effect on the overall trend of the derived speciation rates.", fig.height = 4, fig.width = 4}
plot( alt_models )
```
Next, we can inspect the summary of the trends in the speciation rate through time
```{r summary-mu-twoepoch, fig.cap = "Summary of directional trends in the congruence class.", fig.height = 3, fig.width = 6}
p <- summarize.trends(alt_models, threshold = 0.02)
plot(p)
```
Again, for our threshold of $\epsilon = 0.02 \Delta\lambda$ per million years, it is evident that the constructed congruent models are in complete agreement of the two increases in speciation rate, as well as the more recent decline in speciation rate.

### Reversed trend extinction rates

As a last example, let us explore the case where the extinction rate had the opposite trend to the estimated speciation rate.
That is, if the speciation rate was going up, now we model that the extinction rate goes down.
This should produce the same pattern in the net-diversification rate.
We create this reversed/opposite trend extinction rate function by subtracting the observed speciation rate from some arbitrary baseline (which is necessary to make sure that all extinction rates are positive).








```{r }
mu_scaling <- c(1.1, 1.2, 1.5, 2.0)
mu_0s <- max(lambda(my_model$times)) * mu_scaling

mu_prime <- list()
for (i in seq_along(mu_scaling)){
  mu_prime[[i]] <- local({
    mu_0 <- mu_0s[i]
    function(t) mu_0 - lambda(t)
  })
}
```
Now we compute the new speciation rate given the extinction rates.
```{r eval=TRUE}
alt_models <- congruent.models(my_model, mus = mu_prime)
```
Let's have a look at these *alternative models*.
```{r congruent-models-mu-rev-trend-add, fig.cap = "Diversification rate functions for our reversed trend additive extinction rate function. The solid line shows the reference model and the dashed lines show the alternative model. Interestingly, the reversed model where the extinction rate goes down instead of the speciation rate going up is not exactly in the congruence class, as we see that the net-diversification rates are not identical. Instead, we see that the speciation rate now drops first just to still increase towards the present.", fig.height = 4, fig.width = 4}
plot( alt_models )
```

In Figure \@ref(fig:congruent-models-mu-rev-trend-add) we used an additive trend by subtracting the rate.
However, we can also use a multiplicative model by subtracting the log speciation rate and then back-transforming the extinction rate from the log-scale.
```{r }
mu_min = 0.0075
mu1 <- function(t) exp(log(mu_min) - log(my_model$lambda(t)))
```
Now we compute the new speciation rate given the extinction rates.
```{r eval=TRUE}
alt_models <- congruent.models( my_model, mus = mu1 )
```
Let's have a look at these *alternative models*.
```{r congruent-models-mu-rev-trend-mult, fig.cap = "Diversification rate functions for our reversed multiplicative trend extinction rate function. The solid line shows the reference model and the dashed lines show the alternative model. Interestingly, the derived speciation rate is qualitatively identical except that the speciation rate first decreases before sharply increasing towards the present.", fig.height = 4, fig.width = 4}
plot( alt_models )
```
Since there are just two models here, it is not necessary to plot the trend summary. It is obvious that the speciation rate still exhibits a sharp increase followed by a downturn toward the present.

This concludes our brief exploration of different common shapes of extinction rates.
In summary, the speciation and net-diversification rates are mostly robust to our relatively minor variations in extinction rates.
Specifically, the observed pattern of a sharp increase of the speciation rate towards the present is robust to all variation of extinction rates we tested.
Only if we would change the extinction rates in extreme and biologically unrealistic ways, e.g., by allowing the extinction rate to become implausible large, then the extinction rates can overwhelm our observed diversification rate pattern.


## Testing alternative speciation rates

Let us now look at the congruence class from the angle of the speciation rates.
Can we also specify different speciation rate function and obtain valid extinction rates?

### Constant speciation rates

Let us assume that the speciation rate was constant.
However, the congruence class as defined by @Louca2020 dictates that $\lambda(0) = \lambda'(0)$, assuming the same present day sampling fraction for both models.
Thus, there is only one possible constant rate speciation model within the congruence class with $\lambda'(t) = \lambda(0)$
Specifically, let's assume that the true speciation rate was $\lambda(t) = \lambda_0$.
We follow exactly the same approach as with the extinction rates.
Hence, we start by defining the speciation rate function.
```{r }
lambda1 <- function(t) my_model$lambda(0.0)
```
Now we compute the new extinction rate given the speciation rates.
```{r eval=TRUE}
alt_models <- congruent.models( my_model, lambda = lambda1)
```
Let's have a look at these *alternative models*.
```{r congruent-models-lambda-constant, fig.cap = "Diversification rate functions for a constant speciation rates within the congruence class. The solid line shows our reference model and the different dashed line shows the alternative model. Most interesting, we observe that for the constant speciation rate function we compute an extinction rate that was negative towards the present.", fig.height = 4, fig.width = 4}
plot( alt_models )
```

In this example we observe that for the constant speciation rate function the derived extinction rate function was negative at the near present.
A negative extinction rate is clearly unrealistic.
Thus, we can conclude that there is no constant rate speciation rate function within the congruence class.

### Exponential increasing/decreasing speciation rates

Let us assume that the speciation rate was exponentially increasing or decreasing.
Furthermore, let us try a range of different factors of increase/decrease.
Specifically, let's assume that the growth rate was $\alpha=\{-0.04,-0.03,\ldots, 0.04\}$.
Again, we create the list of speciation rate functions.
```{r }
lambda_present = lambda(0)
alpha_vals <- c(seq(0.0,0.04,0.02),seq(-0.5,-0.1,0.2))
lambda1 <- list()
for (i in 1:length(alpha_vals)) {
  lambda1[[i]] = local({
    alpha = alpha_vals[i]
    function(t) { lambda_present * exp(alpha*t) }
  })
}
```
Now we compute the new extinction rate given the speciation rates.
```{r eval=TRUE}
alt_models <- congruent.models( my_model, lambda = lambda1)
```
Let's have a look at these *alternative models*.
```{r congruent-models-lambda-exp, fig.cap = "Diversification rate functions for exponentially increasing/decreasing speciation rate functions within the congruence class. The solid line shows the reference model and the dashed lines show the alternative models. Again, we observe that for all speciation rate functions that the derived extinction rate functions are negative towards the present and/or at some other point in time.", fig.height = 4, fig.width = 4}
plot( alt_models )
```
We can also plot the directional trends in the extinction rate through time:
```{r lambda-exp-trendmu, fig.cap = "Congruence class assuming exponentially changing speciation rates.", fig.width = 5, fig.height = 3}
p <- summarize.trends(alt_models, threshold = 0.02, rate_name = "mu")
plot(p)
```

Our results here for the exponentially increasing speciation functions show again that the derived extinction rates are negative towards the present and/or at another point in time.
Because of the initial condition that all speciation rate functions within the congruence class must have the same value at time $t=0$ (the present), it is very challenging to find other speciation rate functions that have valid extinction rates (i.e., $\mu(t)\leq 0$).
This indicates that the congruence class is much more restrictive for the space of valid speciation rates.

# 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.

## 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.

## Independent or joint rate sampling

CRABS offers 3 approaches to sample rates:

1. generate speciation rates, then construct congruent extinction rates;
2. generate extinction rates, then construct congruent speciation rates;
3. generate a $p$-trajectory, which corresponds to a unique pair of congruent speciation and extinction rates.
Independent rate sampling (1 and 2) can be done using `sample.basic.models` or `sample.rates`, while joint sampling (3) can be done using `sample.basic.models.joint`.

The choice between independent and joint rate sampling depends on the assumptions you're willing to make. Independent sampling allows you to explore specific classes of scenarios for either speciation and extinction rates, while the joint sampling approach offers more flexibility by imposing less constraints on the rates. @Andreoletti2023 showed that the joint sampling approach is better able to recover trajectories where the congruence class contains signal for temporal variations both in speciation and extinction rates (which by construction cannot be captured with independent sampling of either rates).

The 3 approaches are presented in more details bellow.

## Sample alternative extinction rates

### Temporally autocorrelated

Here, we generate a variety of alternative extinction rates using HSMRF models [@Magee2020].
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.

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

```{r sample-extinction-HSMRF, cache=TRUE, results='hide'}
samples <- sample.congruence.class(
  my_model,
  num.samples=10,
  rate.type="extinction",
  sample.extinction.rates=extinction_rate_samples)
```
We can plot the samples using
```{r cg-extinction-HSMRF-plot, fig.align='center', fig.cap = "Rate functions assuming HSMRF-distributed extinction rates"}
plot(samples)
```
```{r eval=TRUE, fig.cap = "Congruence class assuming linearly changing extinction rates.", fig.height = 3, fig.width = 6}
p <- summarize.trends(samples,
                      threshold = 0.02)
plot(p)
```

As we can see, the alternative models are not in perfect agreement about the entire history of the rates.
Nonetheless, the increases and decrease near the present are preserved, indicating robustness to this form of extinction-rate variation.

### Exponentially decreasing

Here, we generate a variety of alternative extinction rates using HSMRF models [@Magee2020] 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.

```{r sample-function-extinction-exponential}
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)
}
```

```{r sample-extinction-exponential, cache=TRUE, results='hide'}
samples <- sample.congruence.class(
  my_model,
  num.samples=10,
  rate.type="extinction",
  sample.extinction.rates=extinction_rate_samples)
```
We can plot the samples using
```{r cg-extinction-exponential-plot, fig.cap = "Rate functions assuming exponential trends in the extinction rate"}
plot(samples)
```
```{r cg-extinction-exponential, fig.cap = "Congruence class assuming exponentially changing extinction rates", fig.height = 3, fig.width = 6}
p <- summarize.trends(samples, threshold = 0.02)
plot(p)
```
As we can see, the alternative models are not in perfect agreement about the entire history of the rates.
Nonetheless, the increases and decrease near the present are preserved, indicating robustness to this form of extinction-rate variation.

### Linearly increasing

Here, we generate a variety of alternative extinction rates using GMRF models [@Magee2020] 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.

```{r sample-function-extinction-linear}
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)
}
```

```{r sample-extinction-linear, cache=TRUE, results='hide'}
samples <- sample.congruence.class(
  my_model,
  num.samples=10,
  rate.type="extinction",
  sample.extinction.rates=extinction_rate_samples)
```
We can plot the samples using
```{r congruence-class-extinction-linear, fig.cap = "Congruence class assuming linearly changing extinction rates."}
plot(samples)
```
```{r cg-extinction-linear, fig.cap = "Congruence class assuming linearly changing extinction rates.", fig.width = 6, fig.height = 3}
p <- summarize.trends(samples, threshold = 0.02)
plot(p)
```

As we can see, the increases and decrease near the present are preserved, indicating robustness to this form of extinction-rate variation.
Compared to the exponential and HSMRF models above, there is even better agreement between the alternative models.

### 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.

```{r sample-function-extinction-episodic}
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)
}
```

```{r sample-extinction-episodic, cache=TRUE, results='hide'}
samples <- sample.congruence.class(
  my_model,
  num.samples=10,
  rate.type="extinction",
  sample.extinction.rates=extinction_rate_samples)
```
We can plot the samples using
```{r cg-extinction-episodic-plot}
plot(samples)
```
```{r extinction-episodic, fig.cap = "Congruence class assuming episodically changing extinction rates.", fig.width = 6, fig.height = 3}
p <- summarize.trends(samples, threshold = 0.02)
plot(p)
```

As we can see, the increases and decrease near the present are preserved, indicating robustness to this form of extinction-rate variation.
The agreement between the alternative models is comparable to the linear scenarios above, and higher than for the exponential and HSMRF models.

### Uniform independent

The previous functions used `sample.basic.models` to construct a variety of pre-defined hypotheses.
`CRABS` 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.

```{r sample-function-extinction-uniform}
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)
}
```

```{r sample-extinction-uniform}
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
```{r uniform-mu-sample, fig.cap = "Rate functions for the congruence class with proposed alternative extinction rates whose values are drawn from an independent uniform distribution for each episode."}
plot(samples)
```


Next we can plot the directional trends
```{r extinction-uniform-spaghetti, fig.cap = "Congruence class assuming independent and uniformly distributed extinction rates.", fig.width = 6, fig.height = 3}
p <- summarize.trends(samples, threshold = 0.02)
plot(p)
```

As we can see, the alternative models present many distinct histories of the rates through time.
Nonetheless, even in this extreme scenario, the increases and decrease near the present are preserved.


### Brownian motion with a trend

It is also possible to specify fully customized rate sampling functions, and use these to sample the congruence class in `CRABS`. Suppose we thought that the extinction rate was tiny when the clade emerged, but steadily rose in a linear fashion towards the present. Additionally, suppose we allow for some random variation around this linear increase. This can be implemented as a Brownian motion with a trend, as follows.

```{r sample-bm-trend-decreasing}
med_times <- seq(0, max(my_model$time), length.out = 30)

extinction_rate_samples <- function(){
  noise <- cumsum(rnorm(30, mean = 0, sd = 0.03))
  y <- -0.008*med_times + noise + 0.6
  y <- ifelse(y < 0.001, 0.001, y)
  f <- approxfun(med_times, y)
}
```
For simplicity, we restrict rates that are very small or negative to be a minimum of 0.001 units. As before, we can use `CRABS` to sample from the congruence class using this form of proposed extinction rates, and plot the rates:

```{r bm-trend-spaghetti, fig.cap = "Rate plots for the congruent models that were sampled assuming an extinction rate modeled as a Brownian motion with a linearly increasing trend"}
samples <- sample.congruence.class(
  my_model,
  num.samples=15,
  rate.type="extinction",
  sample.extinction.rates=extinction_rate_samples)
plot(samples)
```

The summary of trends across these models contain quite a bit of variation. However, the three main signals in the primate reference model are all recovered: two increases and one decrease in the recent part of the phylogeny  \@ref(fig:bm-trend-summary).

```{r bm-trend-summary, fig.cap = "A summary of directional trends through time across the congruent models. The proposed extinction rate is modeled as a Brownian motion with a trend that increases linearly towards the present.", fig.width = 6, fig.height = 3}
p <- summarize.trends(samples, threshold = 0.02)
plot(p)
```


## 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.

### Temporally autocorrelated

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

```{r sample-speciation-hsmrf}
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.
```{r sample-extinction, cache=TRUE, results='hide'}
samples <- sample.congruence.class(
  my_model,
  num.samples=10,
  rate.type="speciation",
  sample.speciation.rates = speciation_rate_samples)
```

We can plot the samples using
```{r sample-extinction-rateplot, fig.cap = "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."}
plot(samples)
```

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.
```{r sample-speciation-mus, fig.height=5, fig.width=5, fig.cap = "Inferred congruent extinction rates, zoomed in for clarity."}
p <- plot(samples)[[2]] + 
  coord_cartesian(ylim = c(-1, 1)) +
  xlab("time before present") +
  ylab("rate")
plot(p)
```

The models depicted in \@ref(fig:sample-speciation-mus) 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
```{r sample-extinction-trendplot, fig.cap = "A summary of directional trends through time across the congruent models.", fig.width = 6, fig.height = 3}
p <- summarize.trends(samples, threshold = 0.02, rate_name = "mu")
plot(p)
```

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)$)

```{r netdiv-trends, fig.cap = "Directional trends in the net-diversification rate. The congruence class was constructed using alternative speciation rates, and all congruent models have negative extinction rates.", fig.width = 5, fig.height = 3}
p <- summarize.trends(samples, threshold = 0.02, rate_name = "delta")
plot(p)
```

## Jointly sample alternative speciation and extinction rates

Here, we simultaneously generate a variety of congruent speciation and extinction rates [@Andreoletti2023].
We use the function `sample.basic.models.joint` to construct these rates. As explained in section \@ref(samplejoint), this function first samples $p$-trajectories -- controlling the degree to which the signal in the congruence class is driven more by speciation rate variations ($p=1$) or extinction rate variations ($p=0$) -- and then constructs the corresponding speciation and extinction rates.

The procedure for drawing, centering and rescaling a $p$-trajectory is explained in @Andreoletti2023.

The initial value $p_0$ is sampled from a Beta distribution. By default, we use $\beta_1=0.3$ and $\beta_2=0.3$ to favor starting points near $p=0$ or $p=1$. Below we set $\beta_1=0.5$ and $\beta_2=0.2$ to place more emphasis on values close to $p=1$ (constant extinction), ensuring that we properly explore these scenarios. The speciation rate at present is fixed to the reference scenario, while the extinction rate at present is drawn from a log-normal distribution.

### Constant joint trajectory

First, we set the scaling factor of the Markov Random Field (`mrf.sd.scale`) to $0$, such that we obtain constant $p$-trajectories.

```{r, fig.cap = "Congruent diversification scenarios, jointly sampled with constant p-trajectories.", fig.width=6, fig.height=4}
joint_rate_samples <- function(){sample.basic.models.joint(
  times = seq(0, max(primates_ebd$time), length.out=300),
  p.delta = my_model$p.delta,
  mrf.sd.scale = 0,
  beta.param = c(0.5,0.2),
  lambda0 = my_model$lambda(0),
  mu0.median = my_model$mu(0))}

samples_primates_CT <- sample.congruence.class(
  my_model,
  num.samples=20,
  rate.type="joint",
  sample.joint.rates=joint_rate_samples)

plot(samples_primates_CT)
```

### HSMRF joint trajectory

Then, we revert to the default scaling factor ($\sigma_{\text{MRF}}=1$) and use the HSMRF sampler for $p$-trajectories.

```{r, fig.cap = "Congruent diversification scenarios, jointly sampled with HSMRF p-trajectories.", fig.width=6, fig.height=4}
joint_rate_samples <- function(){sample.basic.models.joint(
  times = seq(0, max(primates_ebd$time), length.out=500),
  #times = my_model$times,
  p.delta = my_model$p.delta,
  MRF.type = "HSMRF",
  beta.param = c(0.5,0.2),
  lambda0 = my_model$lambda(0),
  mu0.median = my_model$mu(0))}

samples_primates_HST <- sample.congruence.class(
  my_model,
  num.samples=20,
  rate.type="joint",
  sample.joint.rates=joint_rate_samples)

plot(samples_primates_HST)
```

### HSMRF joint trajectory with global regularity filter

Finally, we sample 1000 HSMRF trajectories, with the aim of selecting the smoothest ones (as a proxy for simplicity or parsimony).

```{r}
joint_rate_samples <- function(){sample.basic.models.joint(
  times = seq(0, max(primates_ebd$time), length.out=500),
  p.delta = my_model$p.delta,
  MRF.type = "HSMRF",
  beta.param = c(0.5,0.2),
  lambda0 = my_model$lambda(0),
  mu0.median = my_model$mu(0))}

samples_FHST <- sample.congruence.class(
  my_model,
  num.samples=1000,
  rate.type="joint",
  sample.joint.rates=joint_rate_samples)
```

The resulting trajectories are filtered according to a given regularity penalty (by default a $\mathcal{L}_1$ norm) and several thresholds.

```{r, fig.cap = "Congruent diversification scenarios, jointly sampled with HSMRF p-trajectories and selected based on their regularity, for decreasingly stringent thersholds.", fig.width=7, fig.asp=0.5}
full.plot.regularity.thresholds(samples_FHST, filtering_fractions = c(0.05, 0.2, 0.9))
```

Other penalties can be specified, among the following ("L1", "L2" and "L1_derivative"):

- The $\mathcal{L}_1$ norm on all the step-wise variations in both speciation $(\lambda_i)$ and extinction rates $(\mu_i)$, for $i \in [\![0;N]\!]$, where $N$ is the number of time steps :
$$\mathcal{L}_1 = \sum_{i=1}^N \left[abs(\lambda_i-\lambda_{i-1})+abs(\mu_i-\mu_{i-1})\right]$$
The $\mathcal{L}_1$ norm sums the absolute differences between consecutive speciation and extinction rates. Rate trajectories with low $\mathcal{L}_1$  norms tend to have few and small rate variations.

- The $\mathcal{L}_2$ norm on step-wise variations:
$$\mathcal{L}_2 = \sqrt{\sum_{i=1}^N \left[(\lambda_i-\lambda_{i-1})^2+(\mu_i-\mu_{i-1})^2\right]}$$
The $\mathcal{L}_2$ norm, also known as the Euclidean norm, squares step-wise differences before summing and taking the square root. This has the effect of penalizing larger jumps more severely than smaller ones.

- The second-order $\mathcal{L}_1$ norm on derivative shifts in both speciation and extinction rates:
$$\mathcal{L}_1' = \sum_{i=2}^N \left[abs(\lambda_i-2\lambda_{i-1}+\lambda_{i-2})+abs(\mu_i-2\mu_{i-1}+\mu_{i-2})\right]$$
The second-order $\mathcal{L}_1$ norm ($\mathcal{L}_1'$) focuses on the change in the rate of change, effectively measuring the volatility in the trajectory. This offers an additional layer of smoothness by considering how rate changes are themselves changing.

# Exploring the posterior

Until this point, we have only considered point estimates for the fitted birth-death models. If we only consider the maximum-likelihood estimates, or the posterior median rates, it is easy to forget that there is often substantial variability in the estimated speciation and extinction rates. Previously we asked the question, if patterns in the posterior median rates are robust to the congruence class. What happens if we want to take into account the estimation error? Inspecting a series of samples in the posterior distribution is one way to assess this question.

## Reading the posterior

In `CRABS`, we have included functionality to read RevBayes log outputs.

```{r read-posterior, results = "hide"}
data(primates_ebd_log)

posterior <- read.RevBayes(primates_ebd_log, n_times = 1000,
                           max_t = 65, n_samples = 20)
```
This yields a list of twenty models that are sampled from the posterior. Each sample is equally spaced on the MCMC iterations, meaning that we minimized the autocorrelation among the posterior samples. In `CRABS` we have included functionality to sample the congruence class, while also taking into account the variation in birth-death models in the posterior distribution. 

## Sampling congruent models in the posterior

The function `sample.congruence.class.posterior` iterates over each sample from the posterior, and draws for each sample another 20 models that are congruent with their respective posterior sample. In this example, we draw 20 alternative extinction rate functions, under an autocorrelated HSMRF distribution starting at a log-normal $\mu_0$ with a median of 0.1. The rejection sampling makes sure that any models with a rate higher than 1.0 is discarded. `sample.congruence.class.posterior` is completely linked to `sample.basic.models`, and any additional parameters are directly passed to `sample.basic.models(...)`.

```{r sample-posterior, results = "hide"}
samples <- sample.congruence.class.posterior(posterior, 
                                             num.samples = 20,
                                             rate.type = "extinction",
                                             rate0.median = 0.1,
                                             model = "MRF",
                                             MRF.type = "HSMRF",
                                             max.rate = 1.0)
```
This example yields a list of 20 congruence classes, represented by their 1 reference posterior model, and additional 20 congruent models each, totaling 420 models. First we will plot the rate functions in the posterior:

```{r posterior-rates, figcap = "The speciation rates for twenty MCMC samples in the posterior distribution."}
## We name each posterior sample for plotting purposes
l <- list()
for (i in seq_along(posterior)){
  df <- model2df(posterior[[i]])  
  df$name <- paste0("posterior", i)
  l[[i]] <- df
}
df_posterior <- bind_rows(l)

col_lambda <- colorspace::sequential_hcl(palette = "Blues", n = length(posterior))

p <- bind_rows(df_posterior) %>%
  dplyr::filter(rate == "Speciation") %>%
  ggplot(aes(x = Time, y = value, color = name)) +
  geom_line() +
  scale_x_reverse() +
  theme_classic() +
  scale_color_manual(values = col_lambda) +
  labs(x = "time before present", y = "speciation rate") +
  theme(legend.position = "none")
plot(p)
```
Evidently, the posterior rates exhibit a lot more rate variation than we saw from the posterior median. What happens if we inspect the twenty congruence classes we just explored, in a visual sense? To do so, we may use the `summarize.posterior` function. Similar to `summarize.trends`, it will plot the directional trends in the speciation rate through time. The only difference is that now, the y-axis is sorted, so that each row no longer represents a model.
```{r summarize-posterior, fig.cap = "Summary of directional trends in 20 posterior samples, and for each of those 20 samples, 20 additional congruent models.", fig.width = 6, fig.height = 4}
p <- summarize.posterior(samples, threshold = 0.02)
plot(p)
```
In contrast to the previous summaries (Section \@ref(alt-ext)), the congruence classes in the posterior samples depict a much more ambiguous patterns of the directional trends. There is still strong support for two sharp increases in speciation rate, followed by a decrease toward the near present, but the signal is not equivocal. In the end, with finite amounts of data, it may seem that non-identifiability (the congruence class) in diversification rate models is not as much of a problem as uncertainty in the estimated diversification rates themselves.

# Estimating birth-death models

One of the main features and aims of `CRABS` is to explore the congruence class for the diversification rates that you estimated using another software or R package.
Here we provide you some examples using standard software for estimating diversification rates through time, but if your preferred software is not included then you can either specify the diversification rate functions manually, or try to follow our examples in loading the output.

## Reading the phylogenetic tree

First, we assume that you have a phylogenetic tree which you can load into R.
There are several tutorials on how to do so (e.g., [TESS](https://cran.r-project.org/web/packages/TESS/vignettes/Bayesian_Diversification_Rate_Analysis.pdf) or [APE](https://cran.r-project.org/web/packages/ape/ape.pdf), and we simply refer to these tutorials here).
As an example, we used the *Primates* phylogeny from @Springer2012.
You can download the data from our [RevBayes tutorial](https://revbayes.github.io/tutorials/divrate/ebd.html).
First we need to load the APE package [@Paradis2019].
```{r load-ape}
library(ape)
```
Now we can read our phylogenetic tree
```{r load-primates}
tree <- read.tree(file="data/primates.tre")
```
Next, read the software examples to fit birth-death models to the tree, and estimate diversification rates.

```{r TESS-analysis2, child = "tess_analysis.Rmd"}
```

```{r treepar-analysis, child = "treepar.Rmd"}
```

```{r RevBayes-analysis, child = "RevBayes.Rmd"}
```

# References

<div id="refs"></div>
