3 Features
3.1 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 3.2 |
sample.basic.models.joint |
Jointly samples speciation and extinction trajectories. More information in section 3.3 |
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.
3.2 sample.basic.models
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 3.3.
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 (Magee et al., 2020), 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.
3.3 sample.basic.models.joint
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.