## ----set_options, include = FALSE--------------------------------------------- knitr::opts_chunk$set( eval = FALSE, # Chunks of codes will not be evaluated by default collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5 # Set device size at rendering time (when plots are generated) ) ## ----setup, eval = TRUE, include = FALSE-------------------------------------- library(deepSTRAPP) is_dev_version <- function (pkg = "deepSTRAPP") { # # Check if ran on CRAN # not_cran <- identical(Sys.getenv("NOT_CRAN"), "true") # || interactive() # Version number check version <- tryCatch(as.character(utils::packageVersion(pkg)), error = function(e) "") dev_version <- grepl("\\.9000", version) # not_cran || dev_version return(dev_version) } ## ----adjust_dpi_CRAN, include = FALSE, eval = !is_dev_version()--------------- knitr::opts_chunk$set( dpi = 50 # Lower DPI to save space ) ## ----adjust_dpi_dev, include = FALSE, eval = is_dev_version()----------------- # knitr::opts_chunk$set( # dpi = 72 # Default DPI for the dev version # ) ## ----cut_phylogeny------------------------------------------------------------ # # ------ Example 1: Cut a regular phylogeny ------ # # # # See help of the dedicated function # ?deepSTRAPP::cut_phylo_for_focal_time() # # ## Load eel phylogeny from the R package phytools # # Source: Collar et al., 2014; DOI: 10.1038/ncomms6505 # library(phytools) # data(eel.tree) # # ## Cut phylogeny # # # Cut tree to 30 Mya # cut_tree_with_tip_labels <- cut_phylo_for_focal_time( # tree = eel.tree, # focal_time = 30, # keep_tip_labels = TRUE) # # ## Show tip labels # # # Because we used 'keep_tip_labels = TRUE', we kept tip.label on terminal branches # # with a unique descending tip. # # # Plot internal node labels on initial tree with cut-off # plot(eel.tree) # abline(v = max(phytools::nodeHeights(eel.tree)[,2]) - 30, col = "red", lty = 2, lwd = 2) # nb_tips <- length(eel.tree$tip.label) # nodelabels_in_cut_tree <- (nb_tips + 1):(nb_tips + eel.tree$Nnode) # nodelabels_in_cut_tree[!(nodelabels_in_cut_tree %in% cut_tree_with_tip_labels$initial_nodes_ID)] <- NA # ape::nodelabels(text = nodelabels_in_cut_tree) # title(main = "Current phylogeny - 0 Mya") # # # Plot initial internal node labels on cut tree # plot(cut_tree_with_tip_labels) # ape::nodelabels(text = cut_tree_with_tip_labels$initial_nodes_ID) # title(main = "Past phylogeny - 30 Mya") # # # Plot cut tree without keeping tip.label on terminal branches with a unique descending tip. # # All tip.labels are converted to their descending/tipward node ID # cut_tree_without_tip_labels <- cut_phylo_for_focal_time( # tree = eel.tree, # focal_time = 30, # keep_tip_labels = FALSE) # plot(cut_tree_without_tip_labels) # title(main = "Past phylogeny - 30 Mya") # ## ----cut_phylogeny_eval, fig.width = 15, fig.height = 10, out.width = "100%", eval = TRUE, echo = FALSE---- ## Load eel phylogeny from the R package phytools # Source: Collar et al., 2014; DOI: 10.1038/ncomms6505 suppressWarnings(library(phytools, quietly = TRUE)) data(eel.tree) ## Cut phylogeny # Cut tree to 30 Mya cut_tree_with_tip_labels <- cut_phylo_for_focal_time( tree = eel.tree, focal_time = 30, keep_tip_labels = TRUE) ## Show tip labels old_par <- par(no.readonly = TRUE) par(mfrow = c(1, 2)) # Plot internal node labels on initial tree with cut-off plot(eel.tree, cex = 0.7) abline(v = max(phytools::nodeHeights(eel.tree)[,2]) - 30, col = "red", lty = 2, lwd = 2) nb_tips <- length(eel.tree$tip.label) nodelabels_in_cut_tree <- (nb_tips + 1):(nb_tips + eel.tree$Nnode) nodelabels_in_cut_tree[!(nodelabels_in_cut_tree %in% cut_tree_with_tip_labels$initial_nodes_ID)] <- NA ape::nodelabels(text = nodelabels_in_cut_tree) title(main = "Current phylogeny - 0 Mya") # Plot initial internal node labels on cut tree plot(cut_tree_with_tip_labels) ape::nodelabels(text = cut_tree_with_tip_labels$initial_nodes_ID) title(main = "Past phylogeny - 30 Mya") par(old_par) ## ----cut_phylogeny_2---------------------------------------------------------- # ## Show edge labels # # # Plot edge labels on initial tree with cut-off # plot(eel.tree) # abline(v = max(phytools::nodeHeights(eel.tree)[,2]) - 30, col = "red", lty = 2, lwd = 2) # edgelabels_in_cut_tree <- 1:nrow(eel.tree$edge) # edgelabels_in_cut_tree[!(1:nrow(eel.tree$edge) %in% cut_tree_with_tip_labels$initial_edges_ID)] <- NA # ape::edgelabels(text = edgelabels_in_cut_tree) # title(main = "Current phylogeny - 0 Mya") # # # Plot initial edge labels on cut tree # plot(cut_tree_with_tip_labels) # ape::edgelabels(text = cut_tree_with_tip_labels$initial_edges_ID) # title(main = "Past phylogeny - 30 Mya") ## ----cut_phylogeny_eval_2, fig.width = 15, fig.height = 10, out.width = "100%", eval = TRUE, echo = FALSE---- ## Show edge labels old_par <- par(no.readonly = TRUE) par(mfrow = c(1, 2)) # Plot edge labels on initial tree with cut-off plot(eel.tree, cex = 0.7) abline(v = max(phytools::nodeHeights(eel.tree)[,2]) - 30, col = "red", lty = 2, lwd = 2) edgelabels_in_cut_tree <- 1:nrow(eel.tree$edge) edgelabels_in_cut_tree[!(1:nrow(eel.tree$edge) %in% cut_tree_with_tip_labels$initial_edges_ID)] <- NA ape::edgelabels(text = edgelabels_in_cut_tree) title(main = "Current phylogeny - 0 Mya") # Plot initial edge labels on cut tree plot(cut_tree_with_tip_labels) ape::edgelabels(text = cut_tree_with_tip_labels$initial_edges_ID) title(main = "Past phylogeny - 30 Mya") par(old_par) ## ----cut_contMap-------------------------------------------------------------- # # ------ Example 2: Cut a contMap ------ # # # # See help of the dedicated function # ?deepSTRAPP::cut_contMap_for_focal_time() # # # A contMap is a phylogeny with estimated continuous ancestral trait values mapped along branches. # # It is typically obtained with `[phytools::contMap()]`. # # Within deepSTRAPP, it is part of the output of `[deepSTRAPP::prepare_trait_data()]` # # when used on a continuous trait. # # ## Load mammals phylogeny and data from the R package motmot (data included in deepSTRAPP) # # Initial data source: Slater, 2013; DOI: 10.1111/2041-210X.12084 # data(mammals, package = "deepSTRAPP") # # # Get the phylogeny # mammals_tree <- mammals$mammal.phy # # Get the continuous trait data # mammals_data <- setNames(object = mammals$mammal.mass$mean, # nm = row.names(mammals$mammal.mass))[mammals_tree$tip.label] # # # Run a stochastic mapping based on a Brownian Motion model # # for Ancestral Trait Estimates to obtain a "contMap" object # mammals_contMap <- phytools::contMap(mammals_tree, x = mammals_data, # res = 100, # Number of time steps # plot = FALSE) # # # Set focal time # focal_time <- 80 # # ## Cut contMap to 80 Mya while keeping tip.label # # on terminal branches with a unique descending tip. # updated_contMap <- cut_contMap_for_focal_time( # contMap = mammals_contMap, # focal_time = focal_time, # keep_tip_labels = TRUE) # # # Plot node labels on initial stochastic map with cut-off # plot_contMap(mammals_contMap) # ape::nodelabels() # abline(v = max(phytools::nodeHeights(mammals_contMap$tree)[,2]) - focal_time, # col = "red", lty = 2, lwd = 2) # title(main = "Current contMap - 0 Mya") # # # Plot initial node labels on cut stochastic map # plot_contMap(updated_contMap) # ape::nodelabels(text = updated_contMap$tree$initial_nodes_ID) # title(main = "Past contMap - 80 Mya") # # # ## Cut contMap to 80 Mya while NOT keeping tip.label. # updated_contMap <- cut_contMap_for_focal_time( # contMap = mammals_contMap, # focal_time = focal_time, # keep_tip_labels = FALSE) # # # Plot node labels on initial stochastic map with cut-off # plot_contMap(mammals_contMap) # ape::nodelabels() # abline(v = max(phytools::nodeHeights(mammals_contMap$tree)[,2]) - focal_time, # col = "red", lty = 2, lwd = 2) # title(main = "Current contMap - 0 Mya") # # # Plot initial node labels on cut stochastic map # plot_contMap(updated_contMap) # ape::nodelabels(text = updated_contMap$tree$initial_nodes_ID) # title(main = "Past contMap - 80 Mya") # ## ----cut_contMap_eval, fig.width = 15, fig.height = 12, out.width = "100%", eval = TRUE, echo = FALSE---- ## Load mammals phylogeny and data data(mammals, package = "deepSTRAPP") # Get the phylogeny mammals_tree <- mammals$mammal.phy # Get the continuous trait data mammals_data <- setNames(object = mammals$mammal.mass$mean, nm = row.names(mammals$mammal.mass))[mammals_tree$tip.label] # Run a stochastic mapping based on a Brownian Motion model # for Ancestral Trait Estimates to obtain a "contMap" object mammals_contMap <- phytools::contMap(mammals_tree, x = mammals_data, res = 100, # Number of time steps plot = FALSE) # Set focal time focal_time <- 80 ## Cut contMap to 80 Mya while keeping tip.label # on terminal branches with a unique descending tip. updated_contMap <- cut_contMap_for_focal_time( contMap = mammals_contMap, focal_time = focal_time, keep_tip_labels = TRUE) ## Plot with tip.labels old_par <- par(no.readonly = TRUE) par(mfrow = c(1, 2)) # Plot node labels on initial stochastic map with cut-off plot_contMap(mammals_contMap, fsize = c(0.5, 1), lwd = 0.7) ape::nodelabels(cex = 0.7) abline(v = max(phytools::nodeHeights(mammals_contMap$tree)[,2]) - focal_time, col = "red", lty = 2, lwd = 2) title(main = "\nCurrent contMap - 0 Mya") # Plot initial node labels on cut stochastic map plot_contMap(updated_contMap) ape::nodelabels(text = updated_contMap$tree$initial_nodes_ID) title(main = "\nPast contMap - 80 Mya") par(old_par) ## Cut contMap to 80 Mya while NOT keeping tip.label. updated_contMap <- cut_contMap_for_focal_time( contMap = mammals_contMap, focal_time = focal_time, keep_tip_labels = FALSE) ## Plot without tip.labels old_par <- par(no.readonly = TRUE) par(mfrow = c(1, 2)) # Plot node labels on initial stochastic map with cut-off plot_contMap(mammals_contMap, fsize = c(0.5, 1), lwd = 0.7) ape::nodelabels(cex = 0.7) abline(v = max(phytools::nodeHeights(mammals_contMap$tree)[,2]) - focal_time, col = "red", lty = 2, lwd = 2) title(main = "\nCurrent contMap - 0 Mya") # Plot initial node labels on cut stochastic map plot_contMap(updated_contMap) ape::nodelabels(text = updated_contMap$tree$initial_nodes_ID) title(main = "\nPast contMap - 80 Mya") par(old_par) ## ----cut_densityMap----------------------------------------------------------- # # ------ Example 3: Cut densityMaps ------ # # # # See help of the dedicated function # ?deepSTRAPP::cut_densityMaps_for_focal_time() # # # A densityMap is a phylogeny with posterior probabilities of a given state mapped along branches. # # It is typically obtained with `[phytools::densityMap()]`. # # Within deepSTRAPP, it is part of the output of `[deepSTRAPP::prepare_trait_data()]` # # when used on a categorical trait. # # deepSTRAPP also offers to overlay densityMaps of all states to obtain a single mapping of all states # # on a unique phylogeny with `[deepSTRAPP::plot_densityMaps_overlay()]` # # #### 1/ Prepare data #### # # ## Load mammals phylogeny and data from the R package motmot (data included in deepSTRAPP) # # Initial data source: Slater, 2013; DOI: 10.1111/2041-210X.12084 # data(mammals, package = "deepSTRAPP") # # # Get the phylogeny # mammals_tree <- mammals$mammal.phy # # Get the continuous trait data # mammals_data <- setNames(object = mammals$mammal.mass$mean, # nm = row.names(mammals$mammal.mass))[mammals_tree$tip.label] # # Convert mass data into categories # mammals_mass <- setNames(object = mammals$mammal.mass$mean, # nm = row.names(mammals$mammal.mass))[mammals_tree$tip.label] # mammals_data <- mammals_mass # mammals_data[seq_along(mammals_data)] <- "small" # mammals_data[mammals_mass > 5] <- "medium" # mammals_data[mammals_mass > 10] <- "large" # table(mammals_data) # # ## Select color scheme for states # colors_per_states <- c("darkblue", "dodgerblue", "lightblue") # names(colors_per_states) <- c("small", "medium", "large") # # # Produce densityMaps using stochastic character mapping based on an equal-rates (ER) Mk model # mammals_cat_data <- prepare_trait_data(tip_data = mammals_data, phylo = mammals_tree, # trait_data_type = "categorical", # evolutionary_models = "ER", # nb_simulations = 100, # colors_per_levels = colors_per_states) # # # Set focal time # focal_time <- 80 # # #### 2/ Plot a unique cut densityMap #### # # # Extract the density map for small mammals (state 3 = "small") # mammals_densityMap_small <- mammals_cat_data$densityMaps[[3]] # # ## Cut densityMap to 80 Mya while keeping tip.label # # on terminal branches with a unique descending tip. # updated_mammals_densityMap_small <- cut_densityMap_for_focal_time( # densityMap = mammals_densityMap_small, # focal_time = focal_time, # keep_tip_labels = TRUE) # # # Plot node labels on initial stochastic map with cut-off # phytools::plot.densityMap(mammals_densityMap_small) # ape::nodelabels() # abline(v = max(phytools::nodeHeights(mammals_densityMap_small$tree)[,2]) - focal_time, # col = "red", lty = 2, lwd = 2) # title(main = "Current densityMap - 0 Mya") # # # Plot initial node labels on cut stochastic map # phytools::plot.densityMap(updated_mammals_densityMap_small) # ape::nodelabels(text = updated_mammals_densityMap_small$tree$initial_nodes_ID) # title(main = "Past densityMap - 80 Mya") # ## ----cut_densityMap_eval, fig.width = 15, fig.height = 12, out.width = "100%", eval = TRUE, echo = FALSE---- ## Load mammals phylogeny and data from the R package motmot (data included in deepSTRAPP) # Initial data source: Slater, 2013; DOI: 10.1111/2041-210X.12084 data(mammals, package = "deepSTRAPP") # Get the phylogeny mammals_tree <- mammals$mammal.phy # Get the continuous trait data mammals_data <- setNames(object = mammals$mammal.mass$mean, nm = row.names(mammals$mammal.mass))[mammals_tree$tip.label] # Convert mass data into categories mammals_mass <- setNames(object = mammals$mammal.mass$mean, nm = row.names(mammals$mammal.mass))[mammals_tree$tip.label] mammals_data <- mammals_mass mammals_data[seq_along(mammals_data)] <- "small" mammals_data[mammals_mass > 5] <- "medium" mammals_data[mammals_mass > 10] <- "large" table(mammals_data) # Select color scheme for states colors_per_states <- c("darkblue", "dodgerblue", "lightblue") names(colors_per_states) <- c("small", "medium", "large") # Produce densityMaps using stochastic character mapping based on an equal-rates (ER) Mk model mammals_cat_data <- suppressMessages(prepare_trait_data(tip_data = mammals_data, phylo = mammals_tree, trait_data_type = "categorical", evolutionary_models = "ER", nb_simulations = 100, colors_per_levels = colors_per_states, plot_map = FALSE, verbose = FALSE)) # Set focal time focal_time <- 80 # Extract the density map for small mammals (state 3 = "small") mammals_densityMap_small <- mammals_cat_data$densityMaps[[3]] ## Cut densityMap to 80 Mya while keeping tip.label # on terminal branches with a unique descending tip. updated_mammals_densityMap_small <- cut_densityMap_for_focal_time( densityMap = mammals_densityMap_small, focal_time = focal_time, keep_tip_labels = TRUE) ## Plot with tip.labels old_par <- par(no.readonly = TRUE) par(mfrow = c(1, 2)) # Plot node labels on initial stochastic map with cut-off phytools::plot.densityMap(mammals_densityMap_small, fsize = c(0.5, 1), lwd = 0.7) ape::nodelabels(cex = 0.7) abline(v = max(phytools::nodeHeights(mammals_densityMap_small$tree)[,2]) - focal_time, col = "red", lty = 2, lwd = 2) title(main = "\nCurrent densityMap - 0 Mya") # Plot initial node labels on cut stochastic map phytools::plot.densityMap(updated_mammals_densityMap_small) ape::nodelabels(text = updated_mammals_densityMap_small$tree$initial_nodes_ID) title(main = "\nPast densityMap - 80 Mya") par(old_par) ## ----cut_densityMaps---------------------------------------------------------- # # #### 3/ Plot set of overlaid densityMaps #### # # ## Cut all densityMaps to 80 Mya while keeping tip.label # # on terminal branches with a unique descending tip. # updated_mammals_densityMaps <- cut_densityMaps_for_focal_time( # densityMaps = mammals_cat_data$densityMaps, # focal_time = focal_time, # keep_tip_labels = TRUE) # # # Plot node labels on initial stochastic map with cut-off # plot_densityMaps_overlay(densityMaps = mammals_cat_data$densityMaps, # colors_per_levels = colors_per_states) # ape::nodelabels() # abline(v = max(phytools::nodeHeights(mammals_cat_data$densityMaps[[1]]$tree)[,2]) - focal_time, # col = "red", lty = 2, lwd = 2) # title(main = "Current overlaid densityMaps - 0 Mya") # # # Plot initial node labels on cut stochastic map # plot_densityMaps_overlay(densityMaps = updated_mammals_densityMaps, # colors_per_levels = colors_per_states) # ape::nodelabels(text = updated_mammals_densityMaps[[1]]$tree$initial_nodes_ID) # title(main = "Past overlaid densityMaps - 80 Mya") # ## ----cut_densityMaps_eval, fig.width = 15, fig.height = 12, out.width = "100%", eval = TRUE, echo = FALSE---- ## Cut all densityMaps to 80 Mya while keeping tip.label # on terminal branches with a unique descending tip. updated_mammals_densityMaps <- cut_densityMaps_for_focal_time( densityMaps = mammals_cat_data$densityMaps, focal_time = focal_time, keep_tip_labels = TRUE) ## Plot with tip.labels old_par <- par(no.readonly = TRUE) par(mfrow = c(1, 2)) # Plot node labels on initial stochastic map with cut-off plot_densityMaps_overlay(densityMaps = mammals_cat_data$densityMaps, colors_per_levels = colors_per_states, lwd = 0.7, fsize = 0.5) abline(v = max(phytools::nodeHeights(mammals_cat_data$densityMaps[[1]]$tree)[,2]) - focal_time, col = "red", lty = 2, lwd = 2) title(main = "\nCurrent overlaid densityMaps - 0 Mya") # Plot initial node labels on cut stochastic map plot_densityMaps_overlay(densityMaps = updated_mammals_densityMaps, colors_per_levels = colors_per_states) title(main = "\nPast overlaid densityMaps - 80 Mya") par(old_par) ## ----cut_BAMM_object---------------------------------------------------------- # # ------ Example 3: Cut a BAMM_object ------ # # # # See help of the dedicated function # ?deepSTRAPP::update_rates_and_regimes_for_focal_time() # # # A BAMM_object is the typical output of a Bayesian Analysis of Macroevolutionary Mixtures (BAMM). # # It is a phylogeny with additional elements describing diversification dynamics, # # including estimates of diversification rates over BAMM posterior samples. # # Within deepSTRAPP, it is the output of `[deepSTRAPP::prepare_diversification_data()]`. # # # Similarly to phylogenies, contMap, and densityMaps, it can be updated for a given focal-time # # by cutting branches at focal-time and updating tip rates. # # # ## Load the BAMM_object summarizing 1000 posterior samples of BAMM # data(whale_BAMM_object, package = "deepSTRAPP") # # # Set focal-time to 5 My # focal_time = 5 # # ## Update the BAMM object # whale_BAMM_object_5My <- update_rates_and_regimes_for_focal_time( # BAMM_object = whale_BAMM_object, # focal_time = 5, # update_all_elements = TRUE, # keep_tip_labels = TRUE) # # # ## Explore updated outputs # # # Extract speciation rates for t = 0My # speciation_tip_rates_0My <- whale_BAMM_object$meanTipLambda # names(speciation_tip_rates_0My) <- whale_BAMM_object$tip.label # speciation_tip_rates_0My # # # Extract speciation rates for t = 5My # speciation_tip_rates_5My <- whale_BAMM_object_5My$meanTipLambda # names(speciation_tip_rates_5My) <- whale_BAMM_object_5My$tip.label # speciation_tip_rates_5My # # Speciation rates have been updated so that they reflect values estimated at the focal-time (t = 5My). # # # Extract extinction rates for t = 0My # extinction_tip_rates_0My <- whale_BAMM_object$meanTipMu # names(extinction_tip_rates_0My) <- whale_BAMM_object$tip.label # extinction_tip_rates_0My # # # Extract extinction rates for t = 5My # extinction_tip_rates_5My <- whale_BAMM_object_5My$meanTipMu # names(extinction_tip_rates_5My) <- whale_BAMM_object_5My$tip.label # extinction_tip_rates_5My # # Extinction rates have been updated so that they reflect values estimated at the focal-time (t = 5My). # # # ## Plot BAMM_object # # # Add "phylo" class to be compatible with phytools::nodeHeights() # class(whale_BAMM_object) <- unique(c(class(whale_BAMM_object), "phylo")) # root_age <- max(phytools::nodeHeights(whale_BAMM_object)[,2]) # # Remove temporary "phylo" class # class(whale_BAMM_object) <- setdiff(class(whale_BAMM_object), "phylo") # # # Plot initial BAMM_object for t = 0 My # plot_BAMM_rates(whale_BAMM_object, add_regime_shifts = TRUE, # labels = TRUE, legend = TRUE, # par.reset = FALSE) # Keep plotting parameters in memory to use abline(). # abline(v = root_age - focal_time, # col = "red", lty = 2, lwd = 2) # title(main = "\nPresent BAMM rates - 0 Mya") # # # Plot updated BAMM_object for t = 5 My # plot_BAMM_rates(whale_BAMM_object_5My, add_regime_shifts = TRUE, # labels = TRUE, legend = TRUE) # title(main = "\nPast BAMM rates - 5 Mya") # ## ----cut_BAMM_object_eval, fig.width = 15, fig.height = 10, out.width = "100%", eval = TRUE, echo = FALSE---- ## Load the BAMM_object summarizing 1000 posterior samples of BAMM data(whale_BAMM_object, package = "deepSTRAPP") # Set focal-time to 5 My focal_time = 5 ## Update the BAMM object whale_BAMM_object_5My <- update_rates_and_regimes_for_focal_time( BAMM_object = whale_BAMM_object, focal_time = 5, update_rates = TRUE, update_regimes = TRUE, update_tree = TRUE, update_plot = TRUE, update_all_elements = TRUE, keep_tip_labels = TRUE, verbose = FALSE) ## Plot BAMM_object # Add "phylo" class to be compatible with phytools::nodeHeights() class(whale_BAMM_object) <- unique(c(class(whale_BAMM_object), "phylo")) root_age <- max(phytools::nodeHeights(whale_BAMM_object)[,2]) # Remove temporary "phylo" class class(whale_BAMM_object) <- setdiff(class(whale_BAMM_object), "phylo") old_par <- par(no.readonly = TRUE) par(mfrow = c(1, 2)) # Plot initial BAMM_object for t = 0 My plot_BAMM_rates(whale_BAMM_object, add_regime_shifts = TRUE, cex = 0.5, # Adjust tip.label size regimes_size = 3, labels = TRUE, legend = TRUE, par.reset = FALSE) # Keep plotting parameters in memory to use abline(). abline(v = root_age - focal_time, col = "red", lty = 2, lwd = 2) title(main = "\nPresent BAMM rates - 0 Mya") # Plot updated BAMM_object for t = 5 My plot_BAMM_rates(whale_BAMM_object_5My, add_regime_shifts = TRUE, regimes_size = 3, labels = TRUE, legend = TRUE) title(main = "\nPast BAMM rates - 5 Mya") par(old_par)