In TGuillerme/dads: Trees and Traits Simulations

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

Thorough example

Here we are going to look at an example on how to simulate complex evolutionary histories with treats.

Let say we want to simulate the following story:

1. A lineage evolves in a normal birth-death fashion (with extinction) from a single origin and with two traits called:
   • "Distance": the physical distance between species in an arbitrary 1D space (evolving in a Brownian motion).
   • "Variable": a 2D uncorrelated Ornstein-Uhlenbeck trait that corresponds to some traits the lineages have. During the whole process, lineages that are more further away from the other species in terms of distance have a higher change of speciating (proportional to that distance).

Three events then happen during the evolutionary history of this lineage:

1. "Correlation event": after reaching time 2 an event happens making the "Variables" trait dimensions become correlated.
2. "Selective extinction event": after reaching a certain point in time, lineages with their "Variables" values that are negative go extinct
3. "Founding event": after reaching a value for the "Distance" trait, a founding event happens with a new birth-death process with no extinction and where the "Variable" is now a 2D correlated Brownian Motion.
4. "Mass extinction event": after reaching a certain number of taxa in the "Founding event", 80% of taxa go extinct randomly.

We are going to run each step one by one in a cumulative manner to see how they pile up:

Step 1:

The birth death parameters

First let's set the overall birth death parameters (a speciation of 1.2 and an extinction of 0.3):

## The birth death parameters
bd_params <- list(speciation = 1.5, extinction = 0.3)
## The stopping rule
stop_rule <- list(max.time = 4)

And let's quickly test if the simulation works by running a birth death tree without traits and anything:

set.seed(0) # 123 also works
## Running a birth death tree 
results <- treats(bd.params = bd_params,
                stop.rule = stop_rule)

## And displaying the results
plot(results)
axisPhylo()

The traits

Then let's create the traits:

## The trait object
two_traits <- make.traits(process = c(BM.process, OU.process), n = c(1, 2),
                          trait.names = c("Distance", "Variables"))

## Visualise the 1D trait (BM.process)
plot(two_traits, trait = 1)

## Visualise the 2D trait (uncorrelated OU.process)
plot(two_traits, trait = 2, use.3D = TRUE)

Note that specifying the correlation matrix is not necessary here per se (since no correlation is the default) but we will modify it through an event later on so it's good practice to have it properly stated before hand (to avoid any confusion).

We can then run the tree with the traits:

set.seed(0)
## Running a birth death tree 
results <- treats(bd.params  = bd_params,
                stop.rule  = stop_rule,
                traits     = two_traits)

## And displaying the results
par(mfrow = c(1,2))
plot(results, trait = 1, main = "Trait 1: Distance")
plot(results, trait = c(2,3), main = "Variables correlation")

## Visualising the tree in 3D
plot(results, trait = c(2,3), use.3D = TRUE)

The modifiers

## Speciation event is more likely if lineage's ancestor is further away from the mean trait value
distance.modify <- function(x, trait.values, lineage) return(x + x * abs(parent.traits(trait.values, lineage)[1] - mean(trait.values[,1])))

## Make a distance modifier
distance.speciation <- make.modifiers(speciation = speciation,
                                      modify = distance.modify)

set.seed(13)
## Running a birth death tree 
results <- treats(bd.params  = bd_params,
                stop.rule  = stop_rule,
                traits     = two_traits,
                modifiers  = distance.speciation)

## And displaying the results
par(mfrow = c(1,2))
plot(results, trait = 1, main = "Trait 1: Distance")
plot(results, trait = c(2,3), main = "Variables correlation")

## Visualising the tree in 3D
plot(results, trait = c(2,3), use.3D = TRUE)

Step 2:

1. "Correlation event": after reaching time 2 an event happens making the "Variable" trait dimensions become correlated.

## The correlation change event
update.correlation <- traits.update(process.args =list(Sigma = matrix(c(10,3,3,2),2,2)))

events_list <- make.events(event.name   = "Correlation",
                           target       = "traits",
                           condition    = age.condition(2),
                           modification = traits.update(
                    process.args =list(Sigma = matrix(c(10,3,3,2),2,2)),
                    trait.names  = "Variables"))

results <- NULL
set.seed(1)
## Running a birth death tree 
results <- treats(bd.params  = bd_params,
                stop.rule  = stop_rule,
                traits     = two_traits,
                modifiers  = distance.speciation,
                events     = events_list)

## And displaying the results
par(mfrow = c(1,2))
plot(results, trait = 1, main = "Trait 1: Distance")
plot(results, trait = c(2,3), main = "Variables correlation")

## Visualising the tree in 3D
plot(results, trait = c(2,3), use.3D = TRUE)

Step 3:

1. "Selective extinction event": after reaching a certain point in time, lineages with their "Variables" values that are negative go extinct

BUGGED FROM HERE!

## Adding the selective event
events_list <- make.events(event.name   = "Selection",
                           add          = events_list,
                           target       = "taxa",
                           condition    = age.condition(4),
                           modification = trait.extinction(x = 0, condition = `<`, trait = 2))

Problem with modification argument here

## Running a birth death tree 
results <- treats(bd.params  = bd_params,
                stop.rule  = stop_rule,
                traits     = two_traits,
                modifiers  = distance.speciation,
                events     = events_list)

## And displaying the results
par(mfrow = c(1,2))
plot(results, trait = 1, main = "Trait 1: Distance")
plot(results, trait = c(2,3), main = "Variables correlation")

Step 4 and 5:

1. "Founding event": after reaching a value for the "Distance" trait, a founding event happens with a new birth-death process with no extinction and where the "Variable" is now a 2D correlated Brownian Motion.
2. "Mass extinction event": after reaching a certain number of taxa in the "Founding event", 80% of taxa go extinct randomly.

## The new trait for the founding event
new_traits <- make.traits(process = c(BM.process, BM.process), n = c(1, 2),
                          trait.names = c("Distance", "Variables"))

## The mass extinction event after reaching 50 taxa
mass_extinction <- make.events(target       = "taxa",
                               condition    = taxa.condition(50),
                               modification = random.extinction(0.8))

## The selective event
events_list <- make.events(event.name   = "Founding", 
                           add          = events_list,
                           target       = "founding",
                           condition    = trait.condition(3),
                           modification = founding.event(traits    = new_traits,
                                                         events    = mass_extinction,
                                                         bd.params = list(
                                                            speciation = 2,
                                                            extinction = 0)),
                           additional.args = list(prefix = "found_"))

## Running a birth death tree 
results <- treats(bd.params  = bd_params,
                stop.rule  = stop_rule,
                traits     = two_traits,
                modifiers  = distance.speciation,
                events     = events_list,
                null.error = 100)

## And displaying the results
par(mfrow = c(1,2))
plot(results, trait = 1, main = "Trait 1: Distance")
plot(results, trait = c(2,3), main = "Variables correlation")

TGuillerme/dads documentation built on July 16, 2025, 9:14 p.m.