--- title: "treats workshop" author: "Thomas Guillerme" date: "`r Sys.Date()`" output: html_document: fig_width: 6 fig_height: 6 --- # Pre-intro Most of the information about the workshop comes from the [`treats` manual](https://tguillerme.github.io/treats). You can also always find more documentation about the functions used here using the `R` inbuilt manual by typing `?function.name`. ## `R` level In this workshop I will assume you are already familiar with basic `R`. The basic notions that I'll assume you know are: * What is a package (e.g. `ape` or `dispRity`) * What is an object (e.g. `this_object <- 1`) * What is an object's class (e.g. the class `"matrix"` or `"phylo"`) * What is a function (e.g. the function `mean(c(1,2))`) * How to access function manuals (e.g. `?mean`) Let's get into it. The workshop markdown file: https://github.com/TGuillerme/treats/blob/master/inst/vignettes/treats_workshop.Rmd The same file for direct download (right click + save as): https://raw.githubusercontent.com/TGuillerme/treats/master/inst/vignettes/treats_workshop.Rmd The same file but in html (for visualising the fancy plots): https://cdn.githubraw.com/TGuillerme/treats/master/inst/vignettes/treats_workshop.html First we'll want to download and install the package: ```{r, eval = FALSE} install.packages("treats") library(treats) ``` ```{r, eval = TRUE, echo = FALSE, message = FALSE} library(treats) ``` ## Acknowledgments Note that I am the only author on the paper associated with this software suggesting I work alone and deserve all credits. This is of course not the case and a lot of input on this project was provided by many patient people. Many thanks to Mark Puttick that inspired the idea of the package; Andrew Beckerman, Natalie Cooper and Gavin Thomas for support during the development stage; Ian Brennan, Gustavo Burin, Christopher Cooney, Alex Cranston, William Gearty, Jasmine Hardie, Tom Lansley, Clement Prieul, Tiago Quental, Joe Rees, Louie Rombaut, James Rule and Sophie Ryan for comments and suggestions on late development stages; and Alain Danet, Tom Lansley, Rob McDonald and Ankita Sinha for testing this workshop # Intro What ca we do with `treats`? * Simulate trees and traits for illustration (toy model, teaching); * Simulate trees and traits as a null model for analysis (we'll see an example of that later); * Whatever you want: this is because of the modular architecture of the package allowing you to modify each steps of the simulations easily (we'll see this throughout the workshop); Here are a couple of examples of things that we can do with `treats`: #### Simulating a tree with increased extinction rate for species with negative trait values after some time: ```{r, echo = FALSE, message = FALSE} my_bd_params <- make.bd.params(speciation = 1, extinction = 0.1) ## birth-death rules (default) default_modifiers <- make.modifiers() ## New condition and new modifier (increasing speciation if trait is negative) negative.trait.condition <- function(trait.values, lineage) { return(parent.traits(trait.values, lineage) < 0) } increase.extinction.1 <- function(x, trait.values, lineage) { return(x + 1) } ## Update the modifier change.speciation <- make.events( condition = age.condition(3), target = "modifiers", modification = modifiers.update(speciation = speciation, condition = negative.trait.condition, modify = increase.extinction.1)) ## My favorite seed set.seed(3) ## Some stopping rules my_stop.rule <- list(max.time = 5) ## The simulation my_simulation <- treats(stop.rule = my_stop.rule, bd.params = my_bd_params, traits = make.traits(), modifiers = default_modifiers, events = change.speciation) plot(drop.singles(my_simulation)) ``` #### Simulating a tree with a 2D trait with changing correlation: ```{r, echo = FALSE, message = FALSE} # A 2D variance covariance matrix cor_matrix <- matrix(1, 2, 2) ## A correlated 2D Brownian Motion correlated_2D_BM <- make.traits(n = 2, process = BM.process, process.args = list(Sigma = cor_matrix)) ## New correlation new_cor <- matrix(c(10,3,3,2),2,2) ## Event changing a trait correlation correlation.change <- make.events( condition = trait.condition(3, absolute = TRUE), target = "traits", modification = traits.update(process.args = list(Sigma = new_cor))) ## My favorite seed set.seed(42) ## Some stopping rules my_stop.rule <- list(max.taxa = 100) ## The simulation my_simulation <- treats(stop.rule = my_stop.rule, traits = correlated_2D_BM, events = correlation.change) ## Plot in 2D # plot(my_simulation, trait = c(1:2)) ## 3D plot in 3D plot(drop.singles(my_simulation), trait = c(1:2), use.3D = TRUE) rgl::rglwidget() ``` #### Simulating a tree with sequential extinction events affect species on the edges of the trait distribution each time: ```{r, echo = FALSE, message = FALSE} library(treats) ## The triggering conditions (every year) trigger.year1 <- age.condition(x = 1, condition = `>=`) trigger.year2 <- age.condition(x = 2, condition = `>=`) trigger.year3 <- age.condition(x = 3, condition = `>=`) trigger.year4 <- age.condition(x = 4, condition = `>=`) ## Trimming the tree trim.the.tree <- function(bd.params, lineage, trait.values) { ## Get all the trait values parent_traits <- parent.traits(trait.values, lineage, current = FALSE) ## Measure the spread of these values boundaries <- sort(range(parent_traits[, 1])) spread <- diff(boundaries) ## Scale the boundaries to 12.5% each side boundaries[1] <- boundaries[1] + spread*0.125 boundaries[2] <- boundaries[2] - spread*0.125 ## Select all the nodes that are not the boundary selected_nodes <- as.numeric(names(which(!(parent_traits[, 1] <= boundaries[2] & parent_traits[, 1] >= boundaries[1])))) ## Make them extinct extinct <- which(lineage$parents %in% selected_nodes) ## Update the lineage object lineage$livings <- lineage$livings[!lineage$livings %in% extinct] lineage$n <- length(lineage$livings) return(lineage) } ## Make christmas events christmas <- make.events(target = "taxa", condition = trigger.year1, modification = trim.the.tree, event.name = "year1") christmas <- make.events(target = "taxa", condition = trigger.year2, modification = trim.the.tree, event.name = "year2", add = christmas) christmas <- make.events(target = "taxa", condition = trigger.year3, modification = trim.the.tree, event.name = "year3", add = christmas) christmas <- make.events(target = "taxa", condition = trigger.year4, modification = trim.the.tree, event.name = "year4", add = christmas) ## A beautiful chris set.seed(19) test <- treats(stop.rule = list(max.time = 5), traits = make.traits(), bd.params = make.bd.params(speciation = 1.5, extinction = 0.2), events = christmas, null.error = 100) plot(drop.singles(test), col = c("nodes" = "darkgreen", "fossils" = "gold", "livings" = "darkgreen"), edges = "darkgreen", lwd = 2) points(y = 0, x = 5, pch = 24, col = NULL, bg = "red", cex = 2.6) points(y = 0, x = 5, pch = 25, col = NULL, bg = "red", cex = 2.6) ``` # An example of using some null modeling Say we're looking at an empirical dataset showing trait evolution through time. Here we have (part of) a mammalian dataset from [Beck and Lee 2014](https://royalsocietypublishing.org/doi/full/10.1098/rspb.2014.1278) showing mammalian traits evolution around the K-Pg extinction 66 million years ago (the one that wiped out dinosaurs). We could be interested in the question of whether the mass extinction event freed ecological niches and was responsible for mammalian diversification in terms of traits. To do this we can simply look at the dataset an measure changes in trait diversity through time (disparity). For ecologists: this is very similar to functional diversity or dissimilarity analyses (I would even argue it's identical). ```{r, fig.height = 6, fig.width = 6, echo = FALSE} ## Loading the package and data library(dispRity) data(BeckLee_tree) data(BeckLee_mat99) ## Make the data into treats format for visualising beck_lee_data_3d <- make.treats(tree = BeckLee_tree, data = BeckLee_mat99[, c(1,2)]) plot(beck_lee_data_3d, trait = c(1:2), use.3D = TRUE) rgl::rglwidget() ``` ```{r, fig.height = 6, fig.width = 12, echo = FALSE} ## Creating the time slices time_slices <- chrono.subsets(BeckLee_mat99[, c(1,2)], BeckLee_tree, method = "continuous", model = "proximity", time = seq(from = 96, to = 36, by = -10)) ## Calculating disparity on the two first dimensions only observed_disparity <- dispRity(boot.matrix(time_slices), metric = c(sum, variances), dimensions = c(1,2)) ## Plotting the tree and the disparity through time par(mfrow = c(1,2)) ## Thre tree plot(beck_lee_data_3d, trait = 1, main = "The tree", ylab = "Trait 1") abline(v = BeckLee_tree$root.time - 66, col = "red") ## The disparity plot(observed_disparity, ylab = "Sum of variances", main = "The disparity through time") abline(v = 4, col = "red") ``` Here we can see that, although disparity (trait diversity) changes through time, it doesn't seem that affected by the K-Pg extinction event. This can be due to many causes, namely because of the quality of this dataset (that was not designed for this). But a fair question could be whether mass extinction events in general _CAN_ actually result in changes of disparity through time, regardless of the clade and the data! In other words, could our tools pick up something here? ## Simulating data to test our hypothesis To test our hypothesis **whether mass extinction could liberate niches for trait diversification to occur**, we could simulate some neutral scenarios that would show us what disparity through time (trait diversification) would look like under two neutral scenarios: * What happens to disparity when a random mass extinction occurs? * What happens to disparity when a selective mass extinction occurs? To do so, we would need to simulate a decent amount of datasets under these two specific scenarios and measure disparity from them and then compare it to our observed disparity. This can be done relatively easily in `treats` by design very specific extinction scenarios. For example here we can simulate some similar datasets as the empirical ones with the two extinction scenarios: ```{r, message = FALSE, fig.width = 6, echo = FALSE} ## Setting the stopping rule (stop after 140 time units) stop_rule <- list(max.time = 140) ## Using the birth-death parameters from the observed tree my_bd_params <- make.bd.params(speciation = 0.035, extinction = 0.02) ## A 2D Brownian motion my_traits <- make.traits(process = BM.process, n = 2) ## Creating a random mass extinction random_extinction <- make.events( target = "taxa", condition = age.condition(140-66), modification = random.extinction(0.75)) ## Creating an extinction that removes species with positive trait values positive_extinction <- make.events( target = "taxa", condition = age.condition(140-66), modification = trait.extinction(x = 0, condition = `>=`)) set.seed(123) ## Simulate the tree and traits with a random extinction event sim_rand_extinction <- treats( traits = my_traits, bd.params = my_bd_params, stop.rule = stop_rule, events = random_extinction, null.error = 100, replicates = 50) ## Simulate the tree and traits with a selective extinction event sim_trait_extinction <- treats( traits = my_traits, bd.params = my_bd_params, stop.rule = stop_rule, events = positive_extinction, null.error = 100, replicates = 50) ``` ```{r, echo = FALSE} plot(sim_rand_extinction[[42]], trait = c(1:2), use.3D = TRUE) rgl::rglwidget() ``` ```{r, echo = FALSE} plot(sim_trait_extinction[[50]], trait = c(1:2), use.3D = TRUE) rgl::rglwidget() ``` ```{r, message = FALSE, fig.width = 12, echo = FALSE} par(mfrow = c(1,2)) ## Visualising the difference between both scenarios ## Random extinction plot(sim_rand_extinction[[42]], main = "Random extinction") abline(v = 66, col = "red") abline(h = 0, col = "grey") ## Selective extinction plot(sim_trait_extinction[[50]], main = "Positive extinction") abline(v = 66, col = "red") abline(h = 0, col = "grey") ``` We can then compare the measured disparity in these scenarios to the observed one in the empirical data: ```{r, fig.width = 12, warning = FALSE, message = FALSE, echo = FALSE} ## Remove single nodes simulated at the extinction sim_rand_extinction <- drop.singles(sim_rand_extinction) sim_trait_extinction <- drop.singles(sim_trait_extinction) ## Calculate the dispRity for all the simulations random_extinction_disparity <- dispRitreats(sim_rand_extinction, method = "continuous", model = "proximity", time = seq(from = 96, to = 36, by = -10), metric = c(sum, variances), scale.trees = FALSE) selective_extinction_disparity <- dispRitreats(sim_trait_extinction, method = "continuous", model = "proximity", time = seq(from = 96, to = 36, by = -10), metric = c(sum, variances), scale.trees = FALSE) ## Scale the disparity results (to compare to the observed ones) random_extinction_disparity <- scale.dispRity(random_extinction_disparity) selective_extinction_disparity <- scale.dispRity(selective_extinction_disparity) observed_disparity <- unlist(get.disparity(scale.dispRity(observed_disparity))) ## Plotting the results with the observed disparity par(mfrow = c(1,2)) plot(random_extinction_disparity, main = "Random extinction", ylab = "Scaled sum of variances", ylim = c(0, 1)) abline(v = 4, col = "red") lines(x = 1:7, y = observed_disparity, lty = 2, lwd = 2) plot(selective_extinction_disparity, main = "Selective extinction", ylab = "", ylim = c(0, 1)) abline(v = 4, col = "red") lines(x = 1:7, y = observed_disparity, lty = 2, lwd = 2) ``` So here this example is just to highlight the need for a modular package to simulate trees and traits using very specific simulation scenarios. > Note that here I didn't explain what's going on on the coding level. But you can check out a more detailed version of this case example at the end of this workshop. # The basics: copy-pasting starts here! ```{r} library(treats) ``` ## Growing a tree For the very basics here we're going to see how to grow a tree with the minimum number of options. `treats` is designed to grow trees easily but not with to many default options this is so to encourage users to know what they're doing _at least a minimum_. ### `stop.rule` So what do we need? To grow a tree you need to fix a rule of when the simulation should stop: This is a hard coded rule for stopping the growth of a tree, this can be the number of species (living or fossil) or a certain time of the simulation. Here when I mean time, I mean how long the simulation runs for in arbitrary time units (e.g. million of years, etc...). This is passed to the `stop.rule` argument that is a named list containing integers or numeric values. For example: ```{r} ## Stopping when reaching 200 species stop_200_sp <- list(max.taxa = 200) ## Stopping when reaching 4 time units stop_4_time <- list(max.time = 4) ``` You can then feed these stop rules to the `treats` function and let the magic happen: ```{r} ## Create a tree with 200 taxa tree_200_sp <- treats(stop.rule = stop_200_sp) ## Creating a tree with 4 time units tree_4_time <- treats(stop.rule = stop_4_time) ``` This results in two different trees (`"phylo"` objects) that you can easily visualise with the excellent `"ape"` package: ```{r, fig.width = 12} ## What's in these objects? tree_200_sp tree_4_time ## Displaying the trees op <- par(mfrow = c(1,2)) plot(tree_200_sp, main = "Stop at 200 tips") ; axisPhylo() plot(tree_4_time, main = "Stop at 4 time units") ; axisPhylo() par(op) ``` So the trees here have roughly the same age but a very different shape, we'll get back to that in a moment. One other thing to note for this stopping rule, is that you can stack them: for example, you can ask the simulation to stop when it reaches 200 species _or_ 4 time units. ```{r} ## Multiple rules my_stop_rules <- list(max.taxa = 200, max.time = 4) ## Simulating one tree with these rules my_tree <- treats(stop.rule = my_stop_rules) ## Displaying the resulting tree plot(my_tree) ; axisPhylo() ``` ### `bd.params` As mentioned above, we can have vastly different trees depending on the simulation random seed. You can of course always modify that by modifying the seed. But a more interesting method is to try to control the actual speciation parameter. This is also called the birth or $\lambda$ parameter. You can find details on how the speciation and the extinction parameters are used in the birth-death equations in the [`treats` manual](http://tguillerme.github.io/treats.html) (or anywhere else where there's info about birth-death algorithms) so I'm not gonna go in the details here. Basically, the bigger the speciation parameter, the shorter the tree branches (the bushier the tree). You can set this parameter using the `make.bd.params` function: ```{r, fig.width = 12} set.seed(1) ## Creating a small birth parameter speciation_0.5 <- make.bd.params(speciation = 0.5) ## Creating a big birth parameter (1 is the default by the way) speciation_1.0 <- make.bd.params(speciation = 1.0) ## Simulating the trees my_tree_small <- treats(stop.rule = my_stop_rules, bd.params = speciation_0.5) my_tree_big <- treats(stop.rule = my_stop_rules, bd.params = speciation_1.0) ## Displaying the trees op <- par(mfrow = c(1,2)) plot(my_tree_small, main = "Speciation = 0.5") ; axisPhylo() plot(my_tree_big, main = "Speciation = 1.0") ; axisPhylo() par(op) ``` By the way, throughout the pipeline in the package you can check what's in the object you create. For example, `speciation_0.5` is a `"treats"` and `"bd.params"` object that contains info about the parameters: ```{r} ## What is this? class(speciation_0.5) ## What's in it? speciation_0.5 ``` This allows you to always double check what's in the objects you create! ### The death parameter! Note that this object contained also an extinction parameter (set to 0 by default). This allows to introduce extinction events in your trees! The bigger this parameter, the more fossils in your tree. ```{r, fig.width = 12} ## Creating a small death parameter extinction_0.1 <- make.bd.params(extinction = 0.1, speciation = 1) ## Creating a big death parameter ## (note speciation = 1 is the default) extinction_0.5 <- make.bd.params(extinction = 0.5) ## Simulating the trees set.seed(1) ## Putting the seed here ensures both trees grow the same way few_fossils <- treats(stop.rule = my_stop_rules, bd.params = extinction_0.1) set.seed(1) many_fossils <- treats(stop.rule = my_stop_rules, bd.params = extinction_0.5) ## Displaying the trees op <- par(mfrow = c(1,2)) plot(few_fossils, main = "Extinction = 0.1", show.tip.label = FALSE) ; axisPhylo() plot(many_fossils, main = "Extinction = 0.5", show.tip.label = FALSE) ; axisPhylo() par(op) ``` OK but what's actually happening in the algorithm here? Again there's many ways to illustrate or define the birth-death process (e.g. [wikipedia](https://en.wikipedia.org/wiki/Birth%E2%80%93death_process)) but we can have a look at a simple illustration of each step [here](http://tguillerme.github.io/treats_01.html). > Note that speciation and extinction are fixed values here. But they can also be distributions! For example `make.bd.params(speciation = runif)` creates a `bd.params` object with a speciation always randomly drawn between 0 and 1. You can also link this parameter to the extinction parameter (for example for keep speciation always greater than extinction). Have a look at `?make.bd.params` for more details ### More simulation options for `treats` Just before we get in the more interesting bits, I want to highlight some useful options to deal with stochasticity in `treats`. In fact, sometimes (often?) your trees can "die" before reaching the stop value, especially if the speciation and extinction parameters are close to each other numerically ```{r} ## Creating some dangerous parameters dangerous_params <- make.bd.params(speciation = 1, extinction = 0.9) ``` ```{r, eval = FALSE} ## Simulating a tree that "dies" (returns an error) set.seed(2) dead_tree <- treats(stop.rule = my_stop_rules, bd.params = dangerous_params) ``` This outputs an error message! Aaaah! But fear not, you can handle this error message differently depending on what you need: * you can either set it to null rather than and error (giving you `NULL` rather than an error) using the `null.error = 100` option: ```{r} ## Simulating a tree that "dies" (returns NULL) set.seed(2) dead_tree <- treats(stop.rule = my_stop_rules, bd.params = dangerous_params, null.error = TRUE) dead_tree # is NULL ``` * or you can "brute force" it by keeping trying to generate a tree until it works (or you loose) patience using the `null.error = n` option (where *n* is the number of trials you're happy to wait for: ```{r} ## Simulating a tree that "dies" (try and try again!) set.seed(2) dead_tree <- treats(stop.rule = my_stop_rules, bd.params = dangerous_params, null.error = 100, ## Try up to 100 times! verbose = TRUE) ## Visualise the trials dead_tree # is a tree! ``` Finally, you can also create tree distributions by using the option `replicates = n` that will generate _n_ trees. ## Adding a trait Now this is all nice a well but not very different from other great `R` packages that allow you to grow a birth-death tree. If anything it's maybe even just more cumbersome! But of course, in `treats` you can also do the traits simulation part! To do that we're going to use the `make.traits` function that will create a `"treats"` and `"traits"` object. First we need to choose a trait generation process. This process should be a `function` that intakes the arguments `x0` (the value of the trait in at it's last node) and `edge.length` (the amount of evolutionary time since last node). For example we can create a Brownian Motion (BM) that is process that draws a random value from a normal distribution centered on some value (`x0`) and with a standard deviation that's relative to the evolutionary time spend (`edge.length`). By default it's good practice to set the `edge.length` argument to a specific value (say 1). ```{r} ## Creating a simple BM process my.BM.process <- function(x0, edge.length = 1) { return(rnorm(n = 1, mean = x0, sd = sqrt(edge.length))) } ``` You can then feed this process to `make.traits` and set some specific parameters (e.g. the starting value) or check if it works and visualise it: ```{r} ## Creating the trait object with the BM process ## and a starting value of 0 my_trait <- make.traits(process = my.BM.process, start = 0) ## What is the trait? class(my_trait) ## What's in it? my_trait ## Note the process name is A (you can change that) ## Add a name to the trait object my_trait <- make.traits(process = my.BM.process, start = 0, trait.name = "BM") ## Visualising what the process would look over time plot(my_trait, main = "What does the process look like over time?") ``` > Note that you can of course use much more complex processes (e.g. multidimensional, with correlation, etc...). You can find some inbuilt processes by looking at the documentation for `?trait.process`. Importantly, the function `make.traits` also checks if you have configured your process function correctly and will ping you an error message if not: ```{r, eval = FALSE} ## This doesn't work! badly_formated_trait <- make.traits(process = rnorm) ``` Once you're happy with your process you can simply feed it to `treats` with the other arguments we've covered before: ```{r} ## Generating a tree and a trait tree_and_trait <- treats(stop.rule = my_stop_rules, bd.params = extinction_0.5, traits = my_trait, null.error = 100, verbose = TRUE) ``` The resulting object is now a `"treats"` object (not a `"phylo"` as previously) that contains two main elements: the tree (a `"phylo"` object) and the data (a code `"matrix"`). ```{r} ## What is it? class(tree_and_trait) ## What's in it? tree_and_trait ## The tree: tree_and_trait$tree ## The data head(tree_and_trait$data) ``` You can again also visualise both together using our faithful `plot` function: ```{r} ## Wow! plot(tree_and_trait, main = "A tree with a trait") ``` And that's it! We've simulated a tree and a trait at the same time. The whole point of this package is of course the modularity where you can play around and choose different tree parameters and different traits. In the next part we'll even see how we can complexify all that by making all the birth-death process and the traits interact. Because yes, in this example so far the trait is just "mapped" on the tree: in other words, although both are simulated at the same time, the traits here are not influencing the tree generation. ### Bonus: map.traits This means that you can also just "map" any trait on a tree that you already have. For example, here we're going to map this BM trait on the bird orders example tree from `ape` using the `map.traits` function: ```{r, fig.width = 12} ## Loading the bird.orders tree data(bird.orders) ## map the trait on the tree ## (i.e. simulate the trait using the tree topology) my_bird_BM <- map.traits(traits = my_trait, tree = bird.orders) ## Show the results op <- par(mfrow = c(1,2)) plot(bird.orders, main = "A tree with no traits") plot(my_bird_BM, main = "A tree with a trait") par(op) ``` The resulting object is a `"treats"` one so you can also just quickly generate just the data using `map.traits(traits = my_trait, tree = bird.orders)$data`. ### Bonus: save.steps Note that here we've been generating trait values only at nodes and tips. In some cases, you might be interested to see you trait values generated constantly through time rather than just during speciation and extinction events. For example simulating trait values every _x_ amount of time. You can do this using the `save.steps` in `treats` for example to generate the data every 0.2 time units: ```{r} ## Generating a tree and a trait every 0.2 time steps continuous_trait <- treats(stop.rule = my_stop_rules, bd.params = extinction_0.5, traits = my_trait, save.steps = 0.2, null.error = 100, verbose = TRUE) ## Visualise the tree and trait plot(continuous_trait) ``` Of course though, the bigger the tree and the traits, the slower the algorithm (and the more RAM eaten)! ## Adding modifiers! But let's get back to fun modular stuff: one of the advantages of the `treats` package is that you can modify _all_ steps of the process. The birth-death process is effectively separated into three independent parts: 1. The selection of the lineage of interest. 2. The generation of branch length. 3. The decision whether that lineage goes extinct or speciates. These three steps are encoded by default in the `treats` birth-death algorithm to be: 1. The selection of the lineage of interest: `sample(1, n)` 2. The generation of branch length: `rexp(1, rate = n*(speciation + extinction))` 3. The decision whether that lineage speciates (`TRUE`) or goes extinct (`FALSE`): `runif(1) < speciation/(speciation + extinction)` You can see in more details what's happening by clicking on the modifier at the bottom of the pages [here](http://tguillerme.github.io/treats_01.html). So what's the point of all that and of the modifiers? ### Modifying the speciation process This is the default speciation process that's using the `bd.params` object: ```{r} ## The default speciation process (if TRUE, do speciate) bd.speciating <- function(bd.params) { return(runif(1) <= (bd.params$speciation/(bd.params$speciation + bd.params$extinction))) } ``` Here the speciation process is just depending on the birth-death parameters. But we can be interested in making a simulating a scenario approximating environements limiting capacities: for example, extinction never occurs until there is at least ten species, then it becomes more and more likely the more species in the environement. We can create a modifier that scales the speciation/turnover ratio by the number of species present in the simulation: ```{r} ## The default speciation process (if TRUE, do speciate) density.dependent.sp <- function(bd.params, lineage) { ## Calculate the speciation/turnover ratio (as previously) ratio <- bd.params$speciation/(bd.params$speciation + bd.params$extinction) ## Scale in with the 10s of species return(runif(1) <= ratio*(10/lineage$n)) } ``` > Note here that `lineage` is a new argument introduced here. You can always check what arguments are allowed for internal functions and what internal parts they modify by checking the [manual](http://tguillerme.github.io/treats) "Quick overview" sections. In the case the argument `lineage` points to an internal structure of the `treats` package (no need to worry what's in there for now - but if you're interested, it's all detailed [here](https://github.com/TGuillerme/treats/blob/master/lineage_object_structure.md)). Similarly as for `make.traits` you can create a `"modifiers"` object using `make.modifiers` (again, testing for the validity of your inputs): ```{r} ## A modifier that makes extinction density dependent dd.extinction <- make.modifiers(speciation = density.dependent.sp) ## What is it? class(dd.extinction) ## What's in it? dd.extinction ``` And then feed it to the `treats` function: ```{r, fig.width = 12} set.seed(1) ## Generating a tree and a trait (no modifier) no_modifier <- treats(stop.rule = my_stop_rules, bd.params = extinction_0.5, null.error = 100, verbose = TRUE) set.seed(1) ## Generating a tree and a trait (modified sepciation) modified_speciation <- treats( stop.rule = my_stop_rules, bd.params = extinction_0.5, modifiers = dd.extinction, null.error = 100, verbose = TRUE) ## Plotting the results op <- par(mfrow = c(1,2)) plot(no_modifier, main = "Default birth-death") plot(modified_speciation, main = "Modified speciation") par(op) ``` ## Adding events! You can find a gallery of events [here](https://raw.githack.com/TGuillerme/dads/master/inst/gitbook/_book/specific-examples.html). OK, and finally, the last class of `"treats"` objects are `"events"` created via `make.events`. These work the same way as for `"traits"` and `"modifiers"` but affect the `treats` simulation at a higher level. For example you can use them to simulate mass extinctions after a certain amount of time. Or you can do much more complex things due to the package modular structure: for example, you can use an `"event"` that changes a `"modifier"` depending on some trait values... Let's look at a quick example here. To create an extinction event, you need to set up three aspects: * *What the event is going to affect?* This can be each individual object passed to the simulation (i.e. the `modifiers`, `traits` or `bd.params`) or internal objects like, in our case the we can set to target to `target = "taxa"` (that will affect the internal `"lineage"` object); * *When should the event happen?* This should be a function that is triggered following any specific condition (e.g. some species has a specific trait value, some number of species are present, etc.). In our case we are going to use a pre-implemented function called `age.condition` that generates the function to trigger the event at a specific age: ```{r} ## Making a function that triggers and event after 2 time units age.condition(x = 2, condition = `>=`) ## This is effectively testing whether time is greater or equal than x x = 2 time <- 1 time >= x ## Or, in a more classic function format `>=`(time, x) ``` > Note that for `age.condition`: the age is in units from the past, i.e. since the start of the simulation, not units from the present (i.e. units of time ago). You can always change that by defining it as `age.condition(x = time_in_the_present - time)` * *What the event is going to modify?* This should be also a function that, once the event is triggered, modifies the target of the simulation (e.g. modify the `bd.params` or the `traits` object). In our case, we will use the pre-implemented function `trait.extinction`: ```{r} ## Making a function that removes all species with trait values >= 0 trait.extinction(x = 0, condition = `>=`) ``` We can then create these events the same way as for `make.modifiers` or `make.traits`: ```{r, message = FALSE} ## Creating an extinction that removes species with positive trait values positive_extinction <- make.events( target = "taxa", condition = age.condition(3), modification = trait.extinction(x = 0, condition = `>=`)) ## What's it? class(positive_extinction) ## What's in it? positive_extinction ``` And then simulate a tree and trait with this extinction event: ```{r, message = FALSE} set.seed(7) ## Simulate the tree and traits with a selective extinction event sim_trait_extinction <- treats( traits = my_trait, bd.params = extinction_0.5, stop.rule = my_stop_rules, events = positive_extinction, null.error = 100) ## Selective extinction plot(sim_trait_extinction, main = "Positive extinction") abline(v = 1, col = "red") abline(h = 0, col = "grey") ``` # Bonus: pre-made evolutionary scenarios On the GitHub page of the package you can find [a series of simulation templates](https://github.com/TGuillerme/treats/labels/simulation%20template): these are simplified complex simulation scenarios that you could find useful for your own specific questions. The idea here is to have a library of such templates so if you create one, please share it to the community by following [this template](https://github.com/TGuillerme/treats/issues/new?assignees=&labels=simulation+template&projects=&template=simulation-template.md&title=Simulate+something). # Mega bonus: comparing simulated extinction scenarios to observed ones Here is the detailed script for the simulating extinction scenarios from the start of the workshop. ## Mammalian evolution after the K-Pg extinction One toy example could be to test how a mass extinction event affect species traits evolution. Two potential scenarios could be tested depending on the type of extinction: * The mass extinction is random * The mass extinction targets specific trait values We could simulate both scenarios and compare them to an observed one to see if the trait evolution matches with either of the simulated scenarios. If they do, we could suggest that one of the extinction scenario occurred. Let's first look at the observed example data from the `dispRity` package (from [Beck and Lee 2014](https://royalsocietypublishing.org/doi/full/10.1098/rspb.2014.1278)). We will do a simple disparity through time analyses. Very briefly: we can build a **traitspace** (a multidimensional space that contains all the possible trait combinations for a group of species) and look how species occupied this space through time by looking at traitspace occupancy (**disparity**). This is a common type of analyses in palaebiology but I would argue that it's identical to functional diversity (dissimilarity) analyses in ecology. ```{r, fig.height = 12, fig.width = 12} ## Loading the package and data library(dispRity) data(BeckLee_tree) data(BeckLee_mat99) ## Creating the time slices time_slices <- chrono.subsets(BeckLee_mat99, BeckLee_tree, method = "continuous", model = "proximity", time = seq(from = 96, to = 36, by = -10)) ## Calculating disparity on the two first dimensions only observed_disparity <- dispRity(boot.matrix(time_slices), metric = c(sum, variances), dimensions = c(1,2)) ## Plotting the tree and the disparity through time par(mfrow = c(2,2), bty = "null") ## Thre tree plot(ladderize(BeckLee_tree), show.tip.label = FALSE, main = "The tree") axisPhylo() abline(v = BeckLee_tree$root.time - 66, col = "red") ## The Traitspace plot(time_slices, main = "The traitspace") axisPhylo() abline(v = BeckLee_tree$root.time - 66, col = "red") ## The disparity plot(observed_disparity, ylab = "Sum of variances", main = "The disparity through time") abline(v = 4, col = "red") ``` ## Simulating similar data To simulate the data and the extinction we will need the following: 1. A stopping rule 2. A set of birth death parameters 3. Some traits 4. Some extinction events ### 1. The stopping rule Stopping when reaching 140 time units: ```{r} ## Setting the stopping rule (stop after 140 time units) stop_rule <- list(max.time = 140) ``` ### 2. The birth death parameters ```{r} ## Using the birth-death parameters from the observed tree my_bd_params <- make.bd.params(speciation = 0.035, extinction = 0.02) ``` ### 3. The traits ```{r} ## A 2D Brownian motion my_traits <- make.traits(process = BM.process, n = 2) ``` ### 4. The extinction events ```{r} ## Creating a random mass extinction random_extinction <- make.events( target = "taxa", condition = age.condition(140-66), modification = random.extinction(0.75)) ## Creating an extinction that removes species with positive trait values positive_extinction <- make.events( target = "taxa", condition = age.condition(140-66), modification = trait.extinction(x = 0, condition = `>=`)) ``` ## Simulating the data We can then feed all that to our two simulations: ```{r, message = FALSE, fig.width = 12} set.seed(123) ## Simulate the tree and traits with a random extinction event sim_rand_extinction <- treats( traits = my_traits, bd.params = my_bd_params, stop.rule = stop_rule, events = random_extinction, null.error = 100, replicates = 50) ## Simulate the tree and traits with a selective extinction event sim_trait_extinction <- treats( traits = my_traits, bd.params = my_bd_params, stop.rule = stop_rule, events = positive_extinction, null.error = 100, replicates = 50) par(mfrow = c(1,2)) ## Visualising the difference between both scenarios ## Random extinction plot(sim_rand_extinction[[42]], main = "Random extinction") abline(v = 66, col = "red") abline(h = 0, col = "grey") legend("topleft", pch = NULL, legend = "A", bty = "n", cex = 2) ## Selective extinction plot(sim_trait_extinction[[50]], main = "Positive extinction") abline(v = 66, col = "red") abline(h = 0, col = "grey") legend("topleft", pch = NULL, legend = "B", bty = "n", cex = 2) ``` Tada! > Note that because of the stochastic nature of these simulations we will simulate a series of trees and traits, not just one. This can easily be done in `treats` using the `replicates` option. Here we will do 50 replicates and focus on just 2 traits. ## Comparing both simualted and empirical scenarios And from there, we can compare the simulated data to the empirical on to see if they match or not. In our specific question we can use the `dispRitreats` function that's a handy interface between the `dispRity` package and `treats` and use that to measure disparity in our simulated scenarios. Again, I'm not going into the details there because this kind of comparisons requires much more work and thoughts but it's here just used as a dummy example on how you can use `treats` for your empirical questions! ```{r, fig.width = 12, warning = FALSE, message = FALSE} ## Remove single nodes simulated at the extinction sim_rand_extinction <- drop.singles(sim_rand_extinction) sim_trait_extinction <- drop.singles(sim_trait_extinction) ## Calculate the dispRity for all the simulations random_extinction_disparity <- dispRitreats(sim_rand_extinction, method = "continuous", model = "proximity", time = seq(from = 96, to = 36, by = -10), metric = c(sum, variances), scale.trees = FALSE) selective_extinction_disparity <- dispRitreats(sim_trait_extinction, method = "continuous", model = "proximity", time = seq(from = 96, to = 36, by = -10), metric = c(sum, variances), scale.trees = FALSE) ## Scale the disparity results (to compare to the observed ones) random_extinction_disparity <- scale.dispRity(random_extinction_disparity) selective_extinction_disparity <- scale.dispRity(selective_extinction_disparity) observed_disparity <- unlist(get.disparity(scale.dispRity(observed_disparity))) ## Plotting the results with the observed disparity par(mfrow = c(1,2)) plot(random_extinction_disparity, main = "Random extinction", ylab = "Scaled sum of variances", ylim = c(0, 1)) abline(v = 4, col = "red") lines(x = 1:7, y = observed_disparity, lty = 2, lwd = 2) plot(selective_extinction_disparity, main = "Selective extinction", ylab = "", ylim = c(0, 1)) abline(v = 4, col = "red") lines(x = 1:7, y = observed_disparity, lty = 2, lwd = 2) ```