In TGuillerme/dads: Trees and Traits Simulations

library(treats)
library(rgl)
knitr::knit_hooks$set(webgl = hook_webgl)

Other functionalities {#others}

`make.bd.params` {#makebdparams}

In most examples above, birth-death parameters were set as fixed values, i.e. a speciation rate of $\lambda$ and an extinction rate of $\mu$ throughout the simulations (with maybe some events modifying these rates). However, it is possible to set these rates as specific or changing distributions. You can do this using the make.bd.params function to provide either a vector of values:

## An example where the speciation is randomly sampled among three values
my_bd_params <- make.bd.params(speciation = c(1/3,42))

## Building a tree using this set of parameters
set.seed(123)
plot(treats(stop.rule = list(max.taxa = 50), bd.params = my_bd_params))
## Note the regions in the tree with short branches
## (that's the speciation being 1/3 while the others are speciation = 42)

Or to directly provide a function from which to sample:

## Another example where speciation is drawn from the interval (0, 1)
make.bd.params(speciation = runif)

In this example, the "bd.params" object passed to the treats function will allow the birth-death process to sample the speciation parameter each time it is called (e.g. during the speciation/extinction step, the branch length step, etc.). If using a function, you can fine tune the arguments to be passed to that function using the speciation.args or the extinction.args arguments (as a named list matching the function's arguments):

## Speciation is drawn from the interval (0.5, 1.5)
make.bd.params(speciation = runif,
               speciation.args = list(min = 0.5, max = 1.5))

When using distributions for both the speciation and extinction parameters, you can run into the (usually) undesired problem of having an extinction rate that is higher than your speciation rate and thus all the taxa in your tree dying out. You can avoid this problem by using the joint distribution argument. This will ensure that the sampling of the extinction rate is always lower or equal to the speciation rate.

## Joint speciation and extinction sampled from uniform distribution 
## with speciation always >= to extinction
make.bd.params(speciation = runif, extinction = runif, joint = TRUE)

Finally to avoid negative sampling values, you can use the argument absolute that will make all the sampled values positive. This argument is set to TRUE by default so you shouldn't have to worry about it most of the time unless you specifically need negative sampled values for your parameters.

## Making the speciation sampling always positive
make.bd.params(speciation = rnorm, absolute = TRUE)

The two other possible arguments for this function are test and update that work the same way as make.traits or make.modifiers

If you are a visual person and your bd.params objects are getting a bit too complicated to remember, you can always quickly plot them. The function will sample from the bd.params object and show the results:

par(mfrow = c(2,3))
plot(make.bd.params(), main = "Default birth-death parameters")
plot(make.bd.params(speciation = runif,
               speciation.args = list(min = 0.5, max = 1.5)),
     main = "Uniform speciation between 0.5 and 1.5\n(no extinction)")
plot(make.bd.params(speciation = runif,
                    speciation.args = list(min = 0.5, max = 1.5)),
     main = "Uniform speciation between 0.5 and 1.5\n(no extinction)")
plot(make.bd.params(speciation = runif, extinction = runif, joint = TRUE),
     main = "Joint uniform speciation and extinction")
plot(make.bd.params(speciation = rnorm, extinction = runif,
                    joint = FALSE, abs = FALSE),
     main = "Disjoint normal speciation and uniform extinction")
par(mfrow = c(1,1))

`drop.things` {#dropthings}

You can use the function drop.things to drop specific elements of the tree and data at the same time by providing the argument what = "fossils" for tips that went extinct, or what = "livings" for tips that where alive at the end of the simulation or what = "singles" to drop internal singleton nodes. Alternatively you can use the function aliases drop.fossils, drop.livings or drop.singles for the exact same results:

## A random tree with fossils and traits and internal nodes every 0.5 times
set.seed(3)
my_data <- treats(stop.rule = list(max.taxa = 20),
                  bd.params = list(speciation = 1, extinction = 1/3),
                  traits    = make.traits(), save.steps = 0.5)

## A tree with 20 tips and 54 nodes
my_data$tree
## And a dataset with 74 rows
dim(my_data$data)

## Removing the fossil species
drop.things(my_data, what = "fossils")$tree
dim(drop.fossils(my_data)$data)

## Removing the living species
drop.things(my_data, what = "livings")$tree
dim(drop.livings(my_data)$data)

## Removing the internal nodes
drop.things(my_data, what = "singles")$tree
dim(drop.singles(my_data)$data)

## Removing the internal nodes AND the fossils
drop.singles(drop.fossils(my_data))

`"treats"` internal utilities

The package also provides utilities for internal functions, namely for designing modifiers or events more easily. These functions don't do anything useful on their own but are optimised to be used internally in treats. For all these functions, you can look at the internal manual for an example (i.e. using ?<function_name>)

So far the package has the following internals:

Function | What it does | Where can it be used ----------|--------------|--------------------- parent.traits | selects the trait values of the current lineage's parents (i.e. direct ancestor) | in make.modifiers taxa.condition | provides a trigger for an event dependent on the number of taxa | in make.events age.condition | provides a trigger for an event dependent on time | in make.events trait.condition | provides a trigger for an event dependent on trait values | in make.events