Package 'deepSTRAPP'

Aggregate a list of contMaps into a unique mean/median contMap

Description

Aggregate a list of contMaps into a unique mean/median contMap as produced by the phytools::contMap() function.

Usage

aggregate_contMaps(
  contMaps,
  fun = "mean",
  display_plot = TRUE,
  ...,
  color_scale = NULL,
  verbose = TRUE
)

Arguments

contMaps	List of objects of class "contMap" that represent independent simulations of the evolution of a continuous trait (i.e., continuous stochastic maps).
fun	Character string. Select the aggregating function. Available options are "mean" and "median". Default = "mean".
display_plot	Logical. Whether to plot the resulting aggregated contMap. Default = TRUE.
...	List of named arguments. Additional arguments to be passed down to deepSTRAPP::plot_contMap().
color_scale	Character vector. List of colors to use to build the color scale with grDevices::colorRampPalette() showing the evolution of the continuous trait. From lowest values to highest values. Default (color_scale = NULL) is using the color palette recorded in the contMap$cols item. If none was provided, the rainbow() palette is used. Provide the full argument name to avoid ambiguity with the col graphical argument.
verbose	Logical. Whether to display progress every 100 edges.

Details

The function is primarily designed to average trait values across multiple continuous stochastic maps produced from the same simulation, typically with contsimmap::make.contsimmap() and then converted to a list of contMaps (phytools format) with the convert_contsimmap_to_contMaps() function.

The result is a unique contMap which should be consistent with the contMap produced by phytools::contMap() when interpolating maximum likelihood estimates of ancestral trait values between nodes.

Value

Returns a unique "contMap" object representing the average trait evolutionary history recorded across all continuous stochastic maps provided as input in contMaps.

The resulting aggregated "contMap" is a list with three items:

• ⁠$tree⁠ A list of classes "simmap" and "phylo". The mapped phylogeny including:

  • ⁠$maps⁠ A list of named numeric vectors. Provides the mapping of aggregated trait values along each edge (scaled between 0 and 1000).

  • ⁠$mapped.edge⁠ A numeric matrix. Provides the evolutionary time spent across trait values (columns) along the edges (rows).

• ⁠$cols⁠ A named character vector mapping values to colors used to display trait evolution along the branches.

• ⁠$lims⁠ A numeric vector. Minimum and maximum aggregated trait values recorded.

Author(s)

Maël Doré

References

For continuous stochastic mapping: Martin, B. S., & Weber, M. G. (2026). Stochastic character mapping of continuous traits on phylogenies. Systematic Biology, syag031. doi:10.1093/sysbio/syag031.

See Also

contsimmap::make.contsimmap() phytools::contMap() plot_contMap()

Examples

if (deepSTRAPP::is_dev_version())
{
 ## The R package 'contsimmap' is needed for this example to work.
 # Please install it manually from: https://github.com/bstaggmartin/contsimmap.

 # ----- Prepare data ----- #

 # Load eel phylogeny and tip data from the R package phytools
 # Source: Collar et al., 2014; DOI: 10.1038/ncomms6505

 data("eel.tree", package = "phytools")
 data("eel.data", package = "phytools")

 # Extract body size
 eel_tip_data <- stats::setNames(eel.data$Max_TL_cm,
                                 rownames(eel.data))

 # ----- Run continuous stochastic mapping ----- #

  # (May take several seconds to run)
 eel_contsimmap <- contsimmap::make.contsimmap(
    tree = eel.tree,
    trait.data = eel_tip_data ,
    nsim = 100, res = 100,
    verbose = TRUE)

 dim(eel_contsimmap) # 121 edges, 1 trait, 100 simulations

 # ----- Convert to a list of contMaps ----- #

 eel_contMaps <- convert_contsimmap_to_contMaps(
    contsimmap = eel_contsimmap, verbose = TRUE)

 # ----- Explore contMaps ----- #

 # Plot contMap from simulation n°1
 plot_contMap(eel_contMaps[[1]])
 # Plot contMap from simulation n°50
 plot_contMap(eel_contMaps[[50]])

 # ----- Aggregate all simulations ----- #

 eel_aggregated_contMap <- aggregate_contMaps(
    contMaps = eel_contMaps,
    verbose = TRUE, display_plot = FALSE)

 plot_contMap(eel_aggregated_contMap)

 # ----- Compare to interpolated contMap from phytools ----- #

 # Produce interpolated contMap with phytools
 eel_contMap <- phytools::contMap(tree = eel.tree, x = eel_tip_data,
                                  res = 100, # Number of time steps
                                  plot = FALSE)

 # Plot the interpolated contMap mapping Maximum Likelihood estimates
 # of ancestral trait values interpolated between nodes
 plot_contMap(eel_contMap)
 # Plot the aggregated contMap representing average trait values
 # from true continuous stochastic mapping simulations
 plot_contMap(eel_aggregated_contMap)
 # Both should be visually similar
 
}

---

Template file for BAMM diversification analyses

Description

Template file for BAMM diversification analyses provided as character strings.

Source: bamm-2.5.0 References: http://bamm-project.org/; https://github.com/macroevolution/bamm

Usage

data(BAMM_template_diversification)

Format

A vector of 260 character strings.

Details

A vector of 260 character strings that can be displayed with print(BAMM_template_diversification). This is the template used to generate the 'config_file.txt' controlling settings used for a BAMM run. It provides detailed explanations of the additional_BAMM_settings that can be used in prepare_diversification_data() to customize the BAMM run. It is called internally by prepare_diversification_data() to produce the custom 'config_file' used in the subsequent BAMM run.

References

Authors: Daniel Rabosky (BAMM) & Pascal Title ({BAMMtools}). Modified by Maël Doré for deepSTRAPP.

See Also

BAMM software website: http://bamm-project.org/

---

Build a BAMM object for a deepSTRAPP run

Description

Build a BAMM object of class bammdata based on the output file of a BAMM run that contains a phylogenetic tree and associated diversification rates mapped along branches across BAMM posterior samples.

The BAMM_object output is typically used as input to run deepSTRAPP with run_deepSTRAPP_for_focal_time() or run_deepSTRAPP_over_time().

This is a wrapper of the original BAMMtools::getEventData() function that additionally provides information on:

• the Marginal Shift Probability (MSP) = the probability of a regime shift to occur along each branch.

• the Maximum A Posteriori probability (MAP) configurations among the posterior samples = the configurations of regime shifts that were sampled most frequently (See BAMMtools::getBestShiftConfiguration()).

• the Maximum Shift Credibility (MSC) configurations among the posterior samples = the configurations of regime shift location with the highest product of marginal probabilities across branches (See BAMMtools::maximumShiftCredibility()).

Those additional elements are used by plot_BAMM_rates() to display regime shift probabilities and locations.

This function enables users to integrate the results of their own BAMM analyses into the deepSTRAPP framework. Alternatively, a full BAMM analysis starting from a time-calibrated phylogeny alone can be carried out within deepSTRAPP with prepare_diversification_data().

Usage

build_BAMM_object(
  phylo,
  eventdata,
  burn_in = 0.25,
  nb_posterior_samples = NULL,
  seed = NULL,
  expectedNumberOfShifts,
  MAP_odds_ratio_threshold = 5,
  verbose = FALSE
)

Arguments

phylo	Object of class "phylo" as defined in R package {ape}. Time-calibrated phylogeny that was used to produce the BAMM run. The phylogeny must be rooted and fully resolved.
eventdata	Character string specifying the path to a BAMM event-data file. Alternatively, an object of class data.frame that includes the event data from a BAMM run.
burn_in	Numeric. Proportion of posterior samples removed from the BAMM output to ensure that the remaining samples were drawn once the equilibrium distribution was reached. Default is 0.25
nb_posterior_samples	Integer. Number of posterior samples to extract, after removing the burn-in, in the final BAMM_object to use for downstream analyses. If set to NULL (default), all samples remaining after removing the burn-in will be kept.
seed	Integer. Set the seed to ensure reproducibility when drawing random posterior samples. Default is NULL (a random seed is used).
expectedNumberOfShifts	Integer. The expected number of regime shifts set during the BAMM run as a hyperparameter controlling the exponential prior distribution used to modulate reversible jumps across model configurations in the rjMCMC run. This is needed to compute priors for regime shift along branches.
MAP_odds_ratio_threshold	Numeric. Controls the definition of 'core-shifts' used to distinguish across configurations when fetching the MAP samples. Shifts that have an odds ratio of marginal posterior probability / prior lower than MAP_odds_ratio_threshold are ignored. See BAMMtools::getBestShiftConfiguration(). Default = 5.
verbose	Logical. Whether to display progress in the console. Default = FALSE.

Value

The function returns a BAMM_object of class "bammdata" which is a list with at least 23 elements.

Phylogeny-related elements used to plot a phylogeny with ape::plot.phylo():

• ⁠$edge⁠ Integer matrix. Defines the tree topology by providing rootward and tipward node ID of each edge.

• ⁠$Nnode⁠ Integer. Number of internal nodes.

• ⁠$tip.label⁠ Character vector. Labels of all tips.

• ⁠$edge.length⁠ Numeric vector. Length of edges/branches.

• ⁠$node.label⁠ Character vector. Labels of all internal nodes. (Present only if present in the initial phylo)

BAMM internal elements used for tree exploration:

• ⁠$begin⁠ Numeric vector. Absolute time since root of edge/branch start (rootward).

• ⁠$end⁠ Numeric vector. Absolute time since root of edge/branch end (tipward).

• ⁠$downseq⁠ Integer vector. Order of node visits when using a pre-order tree traversal.

• ⁠$lastvisit⁠ ID of the last node visited when starting from the node in the corresponding position in ⁠$downseq⁠.

BAMM elements summarizing diversification data:

• ⁠$numberEvents⁠ Integer vector. Number of events/macroevolutionary regimes (k+1) recorded in each posterior configuration. k = number of shifts.

• ⁠$eventData⁠ List of data.frames. One per posterior sample. Records shift events and macroevolutionary regimes parameters. 1st line = Background root regime.

• ⁠$eventVectors⁠ List of integer vectors. One per posterior sample. Record regime ID per branch.

• ⁠$tipStates⁠ List of named integer vectors. One per posterior sample. Record regime ID per tip.

• ⁠$tipLambda⁠ List of named numeric vectors. One per posterior sample. Record speciation rates per tip.

• ⁠$tipMu⁠ List of named numeric vectors. One per posterior sample. Record extinction rates per tip.

• ⁠$eventBranchSegs⁠ List of numeric matrices. One per posterior sample. Record regime ID per segment of branches.

• ⁠$meanTipLambda⁠ Named numeric vector. Mean tip speciation rates across all posterior configurations of tips.

• ⁠$meanTipMu⁠ Named numeric vector. Mean tip extinction rates across all posterior configurations of tips.

• ⁠$type⁠ Character string. Set the type of data modeled with BAMM. Should be "diversification".

Additional elements providing key information for downstream analyses:

• ⁠$expectedNumberOfShifts⁠ Integer. The expected number of regime shifts used to set the prior in BAMM.

• ⁠$MSP_tree⁠ Object of class phylo. List of 4 elements duplicating information from the Phylogeny-related elements above, except ⁠$MSP_tree$edge.length⁠ is recording the Marginal Shift Probability of each branch (i.e., the probability of a regime shift to occur along each branch)

• ⁠$MAP_indices⁠ Integer vector. The indices of the Maximum A Posteriori probability (MAP) configurations among the posterior samples.

• ⁠$MAP_BAMM_object⁠ List of 18 elements of class "bammdata" recording the mean rates and regime shift locations found across the Maximum A Posteriori probability (MAP) configurations. All BAMM elements summarizing diversification data hold a single entry describing this mean diversification history.

• ⁠$MSC_indices⁠ Integer vector. The indices of the Maximum Shift Credibility (MSC) configurations among the posterior samples.

• ⁠$MSC_BAMM_object⁠ List of 18 elements of class "bammdata" recording the mean rates and regime shift locations found across the Maximum Shift Credibility (MSC) configurations. All BAMM elements summarizing diversification data hold a single entry describing this mean diversification history.

Note on Bayesian Analysis of Macroevolutionary Mixtures (BAMM)

BAMM is a model of diversification for time-calibrated phylogenies that explores complex diversification dynamics by allowing multiple regime shifts across clades without a priori hypotheses on the location of such shifts. It uses reversible jump Markov chain Monte Carlo (rjMCMC) to automatically explore a vast range of models with different speciation and extinction rates, and different number and location of regime shifts.

BAMM is one option among others for modeling diversification on phylogenies. You may wish to explore alternative models such as LSBDS model in RevBayes (Höhna et al., 2016), the MTBD model (Barido-Sottani et al., 2020), or the ClaDS2 model (Maliet et al., 2019) for your own data. However, you will need Bayesian models that infer regime shifts to be able to perform STRAPP tests (Rabosky & Huang, 2016). Additionally, you need to format the model output such as in BAMM_object, so it can be used in a deepSTRAPP workflow.

Author(s)

Maël Doré

References

For BAMM: Rabosky, D. L. (2014). Automatic detection of key innovations, rate shifts, and diversity-dependence on phylogenetic trees. PloS one, 9(2), e89543. doi:10.1371/journal.pone.0089543. Website: http://bamm-project.org/.

For {BAMMtools}: Rabosky, D. L., Grundler, M., Anderson, C., Title, P., Shi, J. J., Brown, J. W., ... & Larson, J. G. (2014). BAMM tools: an R package for the analysis of evolutionary dynamics on phylogenetic trees. Methods in Ecology and Evolution, 5(7), 701-707. doi:10.1111/2041-210X.12199

See Also

Initial functions in BAMMtools: BAMMtools::getEventData() BAMMtools::getBestShiftConfiguration() BAMMtools::maximumShiftCredibility()

Associated functions in deepSTRAPP: prepare_diversification_data() plot_BAMM_rates() run_deepSTRAPP_for_focal_time() run_deepSTRAPP_over_time()

For a guided tutorial, see this vignette: vignette("model_diversification_dynamics", package = "deepSTRAPP")

Examples

# The key output from a BAMM is the 'event_data.txt' file
# It can be loaded in R to use as input for deepSTRAPP

# The 'whale_event_data.txt' file used for example here is provided
# within deepSTRAPP in the 'extdata' directory

library(phytools)
data(whale.tree)

BAMM_object <- build_BAMM_object(
   phylo = whale.tree,
   eventdata = system.file("extdata", "whale_event_data.txt", package = "deepSTRAPP"),
   burn_in = 0.25, # Remove 25% as burn-in
   nb_posterior_samples = 1000, # Retain 1000 samples
   expectedNumberOfShifts = 1,
   verbose = TRUE)
str(BAMM_object, 1)

---

Compute STRAPP to test for a relationship between diversification rates and trait data

Description

Carries out the appropriate statistical method to test for a relationship between diversification rates and trait data for a given point in the past (i.e. the focal_time). Tests are based on block-permutations: rates data are randomized across tips following blocks defined by the diversification regimes identified on each tip (typically from a BAMM).

Such tests are called STructured RAte Permutations on Phylogenies (STRAPP) as described in Rabosky, D. L., & Huang, H. (2016). A robust semi-parametric test for detecting trait-dependent diversification. Systematic biology, 65(2), 181-193. doi:10.1093/sysbio/syv066.

The function is an extension of the original BAMMtools::traitDependentBAMM() function used to carry out STRAPP test on extant time-calibrated phylogenies.

Tests can be carried out on speciation, extinction and net diversification rates.

deepSTRAPP::compute_STRAPP_test_for_focal_time() can handle three types of statistical tests depending on the type of trait data provided:

Continuous trait data

Tests for correlations between trait and rates carried out with deepSTRAPP::compute_STRAPP_test_for_continuous_data(). The associated test is the Spearman's rank correlation test (See stats::cor.test).

Binary trait data

For categorical and biogeographic trait data that have only two states (ex: 'Nearctic' vs. 'Neotropics'). Tests for differences in rates between states are carried out with deepSTRAPP::compute_STRAPP_test_for_binary_data(). The associated test is the Mann-Whitney-Wilcoxon rank-sum test (See stats::wilcox.test).

Multinominal trait data

For categorical and biogeographic trait data with more than two states (ex: 'No leg' vs. 'Two legs' vs. 'Four legs'). Tests for differences in rates between states are carried out with deepSTRAPP::compute_STRAPP_test_for_multinomial_data(). The associated test for all states is the Kruskal-Wallis H test (See stats::kruskal.test). If posthoc_pairwise_tests = TRUE, post hoc pairwise tests between pairs of states will be carried out too. The associated test for post hoc pairwise tests is Dunn's post hoc pairwise rank-sum test (See dunn.test::dunn.test).

Usage

compute_STRAPP_test_for_focal_time(
  BAMM_object,
  trait_data_list,
  rate_type = "net_diversification",
  uncertainty_strategy = "paired",
  trait_maps_vs_BAMM_samples_list = NULL,
  seed = NULL,
  nb_permutations = NULL,
  alpha = 0.05,
  two_tailed = TRUE,
  one_tailed_hypothesis = NULL,
  posthoc_pairwise_tests = FALSE,
  p.adjust_method = "none",
  trait_data_type_for_stats = NULL,
  return_perm_data = FALSE,
  nthreads = 1,
  print_hypothesis = TRUE
)

Arguments

BAMM_object	Object of class "bammdata", typically generated with update_rates_and_regimes_for_focal_time(), that contains a phylogenetic tree and associated diversification rates across selected posterior samples updated to a specific time in the past (i.e. the focal_time).
trait_data_list	List obtained from extract_most_likely_trait_values_for_focal_time() or extract_all_trait_values_for_focal_time() that contains at least a ⁠$trait_data⁠ element, a ⁠$focal_time⁠ element, and a ⁠$trait_data_type⁠. ⁠$trait_data⁠ is a named vector with the trait data found on the phylogeny at focal_time. ⁠$focal_time⁠ informs on the time in the past at which the trait and rates data will be tested. ⁠$trait_data_type⁠ informs on the type of trait data: 'continuous', 'categorical', or 'biogeographic'.
rate_type	A character string specifying the type of diversification rates to use. Must be one of 'speciation', 'extinction' or 'net_diversification' (default).
uncertainty_strategy	Character string. To select the strategy used to account for uncertainty in estimates. "rates_only": Only accounts for diversification-rate uncertainty across BAMM posterior samples. Uses ML estimates for continuous traits and the most frequent state/range observed across stochastic maps for categorical and biogeographic data. "paired": Default option. Accounts for both diversification-rate and ancestral trait/range reconstruction uncertainty by pairing BAMM posterior samples with stochastic maps. When the number of BAMM samples and stochastic maps differ, random pairing with replacement from the smaller set is used so that all posterior samples and stochastic maps contribute to the analysis. "full": Exhaustive option that accounts for trait/range- and rate- uncertainty by crossing all BAMM posterior samples with all stochastic maps. Accounts for both diversification-rate and ancestral reconstruction uncertainty by evaluating every combination of BAMM posterior sample and stochastic map. WARNING: This exhaustive approach can substantially increase computation time and memory requirements and is therefore recommended only for moderate-sized analyses.
trait_maps_vs_BAMM_samples_list	(Optional) Data.frame of two variables manually providing the names to associate stochastic maps (⁠$trait_map_ID⁠) with BAMM samples (⁠$BAMM_posterior_sample_ID⁠). This is typically used to ensure the same stochastic maps and BAMM samples are used to test across multiple time-steps. Values are the names of the objects such as "Map_X" and "BAMM_X". Default = NULL. * For uncertainty_strategy == 'rates_only', the ⁠$trait_map_ID⁠ must be "Map_ML" as only the ML estimates of trait values/states/ranges are used. * For uncertainty_strategy == 'paired', each pair of stochastic map and BAMM sample will be used once. * For uncertainty_strategy == 'full', all stochastic maps will be matched with all BAMM samples. Those may partly differ from the actual maps and BAMM samples used for the test, as recorded in ⁠$perm_data_df⁠, because invalid maps with not enough states/ranges are discarded.
seed	Integer. Set the seed to ensure reproducibility. Default is NULL (a random seed is used).
nb_permutations	Integer. To select the number of random permutations to perform during the tests. If NULL (default), all BAMM posterior samples will be used once.
alpha	Numeric. Significance level to use to compute the estimate corresponding to the values of the test statistic used to assess significance of the test. This does NOT affect p-values. Default is 0.05.
two_tailed	Logical. To define the type of tests. If TRUE (default), tests for correlations/differences in rates will be carried out with a null hypothesis that rates are not correlated with trait values (continuous data) or equal between trait states (categorical and biogeographic data). If FALSE, one-tailed tests are carried out. For continuous data, it involves defining a one_tailed_hypothesis testing for either a "positive" or "negative" correlation under the alternative hypothesis. For binary data (two states), it involves defining a one_tailed_hypothesis indicating which states have higher rates under the alternative hypothesis. For multinomial data (more than two states), it defines the type of post hoc pairwise tests to carry out between pairs of states. If posthoc_pairwise_tests = TRUE, all two-tailed (if two_tailed = TRUE) or one-tailed (if two_tailed = FALSE) tests are automatically carried out.
one_tailed_hypothesis	A character string specifying the alternative hypothesis in the one-tailed test. For continuous data, it is either "negative" or "positive" correlation. For binary data, it lists the trait states with states ordered in increasing rates under the alternative hypothesis, separated by a greater-than such as c('A > B').
posthoc_pairwise_tests	Logical. Only for multinomial data (with more than two states). If TRUE, all possible post hoc pairwise (Dunn) tests will be computed across all pairs of states. This is a way to detect which pairs of states have significant differences in rates if the overall test (Kruskal-Wallis) is significant. Default is FALSE.
p.adjust_method	A character string. Only for multinomial data (with more than two states). It specifies the type of correction to apply to the p-values in the post hoc pairwise tests to account for multiple comparisons. See stats::p.adjust() for the available methods. Default is none.
trait_data_type_for_stats	(Optional) Character string. One of "continuous", "binary", or "multinomial". Imposes the statistical method to use, instead of deriving it from the type of traits and the number of states/ranges observed at this focal_time. This is set automatically by run_deepSTRAPP_over_time(), which fixes the method once for the whole run from the complete trait mapping, so that the same test is used at every time step and the p-values along the trajectory remain comparable. Default is NULL, in which case the method is selected from the number of states/ranges in trait_data_list$trait_data observed at this focal_time.
return_perm_data	Logical. Whether to return the stats data computed from the posterior samples for observed and permuted data in the output. This is needed to plot the histogram of the null distribution used to assess significance of the test with plot_histogram_STRAPP_test_for_focal_time(). Default is FALSE.
nthreads	Integer. Number of threads to use for parallel computing of the tests across the permutations. The R package parallel must be loaded for nthreads > 1. Default is 1.
print_hypothesis	Logical. Whether to print information on what test is carried out, detailing the null and alternative hypotheses, what significance level is used to reject or not the null hypothesis, and how uncertainty in trait estimates is handled. Default is TRUE.

Details

This set of functions carries out the STructured RAte Permutations on Phylogenies (STRAPP) test as defined in Rabosky, D. L., & Huang, H. (2016). A robust semi-parametric test for detecting trait-dependent diversification. Systematic biology, 65(2), 181-193.

It is an extension of the original BAMMtools::traitDependentBAMM() function used to carry out STRAPP test on extant time-calibrated phylogenies, but allowing here to test for differences/correlations at any point in the past (i.e. the focal_time).

It takes an object of class "bammdata" (BAMM_object) that was updated such that its diversification rates (⁠$tipLambda⁠ and ⁠$tipMu⁠) and regimes (⁠$tipStates⁠) are reflecting values observed at a specific time in the past (i.e. the ⁠$focal_time⁠). Similarly, it takes a list (trait_data_list) that provides ⁠$trait_data⁠ as observed on branches at the same focal_time as the diversification rates and regimes.

A STRAPP test is carried out by drawing a random set of posterior samples from the BAMM_object, then randomly permuting rates across blocks of tips defined by the macroevolutionary regimes. Test statistics are then computed across the initial observed data and the permuted data for each sample. In a two-tailed test, the p-value is the proportion of posterior samples in which the test stats is more extreme in the permuted than in the observed data. In a one-tailed test, the p-value is the proportion of posterior samples in which the test stats is higher in the permuted than in the observed data.

———- Major changes compared to BAMMtools::traitDependentBAMM() ———-

• Allow to account for uncertainty in trait estimates by pairing/mapping multiple trait data extracted from stochastic maps across the BAMM samples. See the uncertainty_strategy argument for details.

• Add post hoc pairwise tests (Dunn test) for multinomial data. Use posthoc_pairwise_tests = TRUE.

• Provide outputs tailored for histogram plots plot_histogram_STRAPP_test_for_focal_time() and p-value time-series plots plot_STRAPP_pvalues_over_time().

• Add prints detailing what test is carried out, what are the null and alternative hypotheses, and what significance level is used to reject or not the null hypothesis. (Enabled with print_hypothesis = TRUE).

• Split the function in multiple sub-functions according to the type of data (⁠$trait_data_type⁠).

• Prevent using Pearson's correlation tests and applying log-transformation for continuous data. The rationale is that there is no reason to assume that tip rates are distributed normally or log-normally. Thus, a Spearman's rank correlation test is favored.

Value

The function returns a list with at least eleven elements.

Summary elements for the main test:

• ⁠$estimate⁠ Named numeric. Value of the test statistic used to assess significance of the test according to the significance level provided (alpha). The test is significant if ⁠$estimate⁠ is higher than zero.

• ⁠$stats_median⁠ Numeric. Median value of the distribution of test statistics.

• ⁠$nb_test_stats⁠ Integer. Number of test stats in the distribution.

• ⁠$p-value⁠ Numeric. P-value of the test. The test is considered significant if ⁠$p-value⁠ is lower than alpha.

• ⁠$method⁠ Character string. The statistical method used to carry out the test.

• ⁠$rate_type⁠ Character string. The type of diversification rates tested. One of 'speciation', 'extinction' or 'net_diversification'.

• ⁠$trait_data_type⁠ Character string. The type of trait data as found in 'trait_data_list$trait_data_type'. One of 'continuous', 'categorical', or 'biogeographic'.

• ⁠$trait_data_type_for_stats⁠ Character string. The type of trait data used to select statistical method. One of 'continuous', 'binary', or 'multinomial'.

• ⁠$focal_time⁠ The time in the past at which the trait and rates data were tested.

• ⁠$uncertainty_strategy⁠ Character string. The strategy used to account for uncertainty in estimates.

• ⁠$trait_maps_vs_BAMM_samples_list⁠ List of two elements recording the stochastic maps (⁠$trait_map_ID⁠) and BAMM samples (⁠$BAMM_posterior_sample_ID⁠) chosen for testing. Those may partly differ from the actual maps and BAMM samples used for the test as recorded in ⁠$perm_data_df⁠ because invalid maps with not enough states/ranges are discarded.

• ⁠$states_observed⁠ (Only for categorical and biogeographic data) Character string vector of the states/ranges actually found in trait_data_list$trait_data at this focal_time. This can be a subset of the states/ranges described in the complete trait mapping, since a state/range may be absent from the deeper time steps.

• ⁠$nb_states_observed⁠ (Only for categorical and biogeographic data) Integer. Number of states/ranges actually found in trait_data_list$trait_data at this focal_time.

If using continuous or binary data:

• ⁠$two-tailed⁠ Logical. Record the type of test used: two-tailed if TRUE, one-tailed if FALSE. If one_tailed_hypothesis is provided (only for continuous and binary trait data):

• ⁠$one_tailed_hypothesis⁠ Character string. Record of the alternative hypothesis used for the one-tailed tests.

If posthoc_pairwise_tests = TRUE (only for multinomial trait data):

• ⁠$posthoc_pairwise_tests⁠ List of at least 3 sub-elements:

  • ⁠$summary_df⁠ Data.frame of six variables providing the summary results of post hoc pairwise tests

  • ⁠$method⁠ Character string. The statistical method used to carry out the test. Here, "Dunn".

  • ⁠$two-tailed⁠ Logical. Record the type of post hoc pairwise tests used: two-tailed if TRUE, one-tailed if FALSE.

If return_perm_data = TRUE, the stats data computed from the posterior samples for observed and permuted data are provided. This is needed to plot the histogram of the null distribution used to assess significance of the test with plot_histogram_STRAPP_test_for_focal_time().

• ⁠$perm_data_df⁠ A data.frame with four variables summarizing the data generated during the STRAPP test:

  • ⁠$trait_map_ID⁠ Integer. ID of the stochastic map from which trait data were extracted and used for the STRAPP test.

  • ⁠$BAMM_posterior_sample_ID⁠ Integer. ID of the posterior samples randomly drawn and used for the STRAPP test.

  • ⁠$*_obs⁠ Numeric. Test stats computed from the observed data in the posterior samples. Name depends on the test used.

  • ⁠$*_perm⁠ Numeric. Test stats computed from the permuted data in the posterior samples. Name depends on the test used.

  • ⁠$delta_*⁠ OR ⁠$abs_delta_*⁠ Numeric. Test stats computed for the STRAPP test comparing observed stats and permuted stats. Name depends on the test used and the type of tests (two-tailed compare absolute values; one-tailed compare raw values). Combined with posthoc_pairwise_tests = TRUE, the stats data are also provided for the post hoc pairwise tests:

• ⁠$posthoc_pairwise_tests$perm_data_array⁠ A 3D array containing stats data for all post hoc pairwise tests in a similar format to ⁠$perm_data_df⁠.

If no STRAPP test was performed in the case of categorical/biogeographic data with a single state/range at focal_time, only the ⁠$trait_data_type⁠, ⁠$trait_data_type_for_stats⁠ = "none", ⁠$focal_time⁠, ⁠$trait_maps_vs_BAMM_samples_list⁠, ⁠$states_observed⁠, and ⁠$nb_states_observed⁠ are returned.

Author(s)

Maël Doré

References

For STRAPP: Rabosky, D. L., & Huang, H. (2016). A robust semi-parametric test for detecting trait-dependent diversification. Systematic biology, 65(2), 181-193. doi:10.1093/sysbio/syv066.

For STRAPP in deep times: Doré, M., Borowiec, M. L., Branstetter, M. G., Camacho, G. P., Fisher, B. L., Longino, J. T., Ward, P. S., Blaimer, B. B., (2025), Evolutionary history of ponerine ants highlights how the timing of dispersal events shapes modern biodiversity, Nature Communications. doi:10.1038/s41467-025-63709-3

See Also

Associated functions in deepSTRAPP: extract_most_likely_trait_values_for_focal_time() update_rates_and_regimes_for_focal_time()

Original function in BAMMtools: BAMMtools::traitDependentBAMM()

Statistical tests: stats::cor.test() stats::wilcox.test() stats::kruskal.test() dunn.test::dunn.test()

For a guided tutorial, see this vignette: vignette("explore_STRAPP_test_types", package = "deepSTRAPP")

Examples

if (deepSTRAPP::is_dev_version())
{
 # ------ Prepare data ------ #

 ## Load the BAMM_object summarizing 1000 posterior samples of BAMM with diversification rates
 # for ponerine ants extracted for 10My ago.
 data(Ponerinae_BAMM_object_10My, package = "deepSTRAPP")
 ## This dataset is only available in development versions installed from GitHub.
 # It is not available in CRAN versions.
# Use remotes::install_github(repo = "MaelDore/deepSTRAPP") to get the latest development version.

 # Plot the associated phylogeny with mapped rates
 plot_BAMM_rates(Ponerinae_BAMM_object_10My)

 ## Load the object containing head width trait data for ponerine ants extracted for 10My ago.
 data(Ponerinae_trait_cont_tip_data_10My, package = "deepSTRAPP")

 # Plot the associated contMap (continuous trait stochastic map)
 plot_contMap(Ponerinae_trait_cont_tip_data_10My$contMap)

 # Check that objects are ordered in the same fashion
 identical(names(Ponerinae_BAMM_object_10My$tipStates[[1]]),
           names(Ponerinae_trait_cont_tip_data_10My$trait_data))

 # Save continuous data
 trait_data_continuous <- Ponerinae_trait_cont_tip_data_10My

 ## Transform trait data into binary and multinomial data

 # Binarize data into two states
 trait_data_binary <- trait_data_continuous
 trait_data_binary$trait_data[trait_data_continuous$trait_data < 0.5] <- "state_A"
 trait_data_binary$trait_data[trait_data_continuous$trait_data >= 0.5] <- "state_B"
 trait_data_binary$trait_data_type <- "categorical"

 table(trait_data_binary$trait_data)

 # Categorize data into three states
 trait_data_multinomial <- trait_data_continuous
 trait_data_multinomial$trait_data[trait_data_continuous$trait_data < 0.6] <- "state_B"
 trait_data_multinomial$trait_data[trait_data_continuous$trait_data < 0.4] <- "state_A"
 trait_data_multinomial$trait_data[trait_data_continuous$trait_data >= 0.6] <- "state_C"
 trait_data_multinomial$trait_data_type <- "categorical"

 table(trait_data_multinomial$trait_data)

 ## Duplicate trait data as if extracted from multiple stochastic maps

 trait_data_continuous_multimaps <- Ponerinae_trait_cont_tip_data_10My
 trait_data_initial <- trait_data_continuous_multimaps$trait_data

 trait_data_multimaps <- list()
 for (i in 1:10)
 {
   # Add a bit of randomness across all stochastic maps data
   trait_data_multimaps[[i]] <- trait_data_initial +
       rnorm(n = length(trait_data_initial), mean = 0, sd = sd(trait_data_initial) / 10)
 }
 names(trait_data_multimaps) <- paste0("Map_", 1:10)
 trait_data_continuous_multimaps$trait_data <- trait_data_multimaps


  # (May take several minutes to run)
 # ------ Compute STRAPP test for continuous data ------ #

 plot(x = trait_data_continuous$trait_data, y = Ponerinae_BAMM_object_10My$tipLambda[[1]])

 # Compute STRAPP test under the alternative hypothesis of a "negative" correlation
 # between "net_diversification" rates and trait data
 STRAPP_results <- compute_STRAPP_test_for_focal_time(
    BAMM_object = Ponerinae_BAMM_object_10My,
    trait_data_list = trait_data_continuous,
    uncertainty_strategy = "rates_only",
    two_tailed = FALSE,
    one_tailed_hypothesis = "negative",
    return_perm_data = TRUE)
 str(STRAPP_results, max.level = 2)
 # Data from the posterior samples is available in STRAPP_results$perm_data_df
 head(STRAPP_results$perm_data_df)

 # ------ Compute STRAPP test for binary data ------ #

 # Compute STRAPP test under the alternative hypothesis that "state_A" is associated
 # with higher "net_diversification" that "state_B"
 STRAPP_results <- compute_STRAPP_test_for_focal_time(
    BAMM_object = Ponerinae_BAMM_object_10My,
    trait_data_list = trait_data_binary,
    uncertainty_strategy = "rates_only",
    two_tailed = FALSE,
    one_tailed_hypothesis = c("state_A > state_B"))
 str(STRAPP_results, max.level = 1)

 # Compute STRAPP test under the alternative hypothesis that "state_B" is associated
 # with higher "net_diversification" that "state_A"
 STRAPP_results <- compute_STRAPP_test_for_focal_time(BAMM_object = Ponerinae_BAMM_object_10My,
    trait_data_list = trait_data_binary,
    uncertainty_strategy = "rates_only",
    two_tailed = FALSE,
    one_tailed_hypothesis = c("state_B > state_A"))
 str(STRAPP_results, max.level = 1)

 # ------ Compute STRAPP test for multinomial data ------ #

 # Compute STRAPP test between all three states, and compute post hoc tests
 # for differences in rates between all possible pairs of states
 # with a p-value adjusted for multiple comparison using Bonferroni's correction
 STRAPP_results <- compute_STRAPP_test_for_focal_time(
    BAMM_object = Ponerinae_BAMM_object_10My,
    trait_data_list = trait_data_multinomial,
    uncertainty_strategy = "rates_only",
    posthoc_pairwise_tests = TRUE,
    two_tailed = TRUE,
    p.adjust_method = "bonferroni",
    return_perm_data = TRUE)
 str(STRAPP_results, max.level = 2)
 # All post hoc pairwise test summaries are available in $summary_df
 STRAPP_results$posthoc_pairwise_tests$summary_df

 # ------ Compute STRAPP test with the 'paired' strategy ------ #

 # Account for uncertainty in trait estimates
 # by pairing stochastic maps with BAMM posterior samples

 # Compute STRAPP test under the alternative hypothesis of a "negative" correlation
 # between "net_diversification" rates and trait data
 STRAPP_results <- compute_STRAPP_test_for_focal_time(
   BAMM_object = Ponerinae_BAMM_object_10My,
   trait_data_list = trait_data_continuous_multimaps,
   uncertainty_strategy = "paired",
   two_tailed = FALSE,
   one_tailed_hypothesis = "negative",
   return_perm_data = TRUE)
 str(STRAPP_results, max.level = 2)
 # Data from the paired stochastic maps and BAMM posterior samples
 # is available in STRAPP_results$perm_data_df
 head(STRAPP_results$perm_data_df)
 # Tests were performed across 1000 iterations since the stochastic maps
 # have been paired with the 1000 BAMM samples
 nrow(STRAPP_results$perm_data_df)
 # Each of the 10 maps was paired ca. 100 times with a unique BAMM sample
 table(STRAPP_results$perm_data_df$trait_map_ID)

 # ------ Compute STRAPP test with the 'full' strategy ------ #

 # Account for uncertainty in trait estimates
 # by crossing stochastic maps with BAMM posterior samples

 # Compute STRAPP test under the alternative hypothesis of a "negative" correlation
 # between "net_diversification" rates and trait data
 STRAPP_results <- compute_STRAPP_test_for_focal_time(
   BAMM_object = Ponerinae_BAMM_object_10My,
   trait_data_list = trait_data_continuous_multimaps,
   uncertainty_strategy = "full",
   two_tailed = FALSE,
   one_tailed_hypothesis = "negative",
   return_perm_data = TRUE)
 str(STRAPP_results, max.level = 2)
 # Data from the combined stochastic maps X BAMM posterior samples
 # is available in STRAPP_results$perm_data_df
 head(STRAPP_results$perm_data_df)
 # Tests were performed across 10000 iterations since the 10 stochastic maps
 # have been combined with the 1000 BAMM samples
 nrow(STRAPP_results$perm_data_df)
 # Each of the 10 maps was paired with all 1000 BAMM samples
 table(STRAPP_results$perm_data_df$trait_map_ID)
 
}

---

Convert Biogeographic Stochastic Map (BSM) to phytools SIMMAP stochastic map (SM) format

Description

These functions convert a Biogeographic Stochastic Map (BSM) output from BioGeoBEARS into a simmap object from R package {phytools} (See phytools::make.simmap()).

They require a model fit with BioGeoBEARS::bears_optim_run() and the output of a Biogeographic Stochastic Mapping performed with BioGeoBEARS::runBSM() to produce simmap objects as phylogenies with the associated mapping of range evolution along branches across simulations.

• convert_BSM_to_simmap(): Produce one simmap for the required simulation (index of the simulation provided with sim_index).

• convert_BSMs_to_simmaps(): Produce all simmap objects for all simulations stored in a unique multiSimmap object.

Initial functions in R package BioGeoBEARS by Nicholas J. Matzke:

• BioGeoBEARS::BSM_to_phytools_SM()

• BioGeoBEARS::BSMs_to_phytools_SMs()

Usage

convert_BSM_to_simmap(model_fit, phylo, BSM_output, sim_index)

convert_BSMs_to_simmaps(model_fit, phylo, BSM_output)

Arguments

model_fit	A BioGeoBEARS results object, produced by ML inference via BioGeoBEARS::bears_optim_run().
phylo	Time-calibrated phylogeny used in the BioGeoBEARS analyses to produce the historical biogeographic inference and run the Biogeographic Stochastic Mapping. Object of class "phylo" as defined in {ape}.
BSM_output	A list with two objects, a cladogenetic events table and an anagenetic events table, as the result of Biogeographic Stochastic Mapping conducted with BioGeoBEARS::runBSM().
sim_index	Integer. Index of the biogeographic simulation targeted to produce the simmap with convert_BSM_to_simmap().

Details

• BioGeoBEARS::BSM_to_phytools_SM()

• BioGeoBEARS::BSMs_to_phytools_SMs()

These functions are slight adaptations of original functions from the R Package BioGeoBEARS by N. Matzke.

Initial functions: BioGeoBEARS::BSM_to_phytools_SM() BioGeoBEARS::BSMs_to_phytools_SMs()

Changes:

• Solves issue with differences in ranges allowed across time-strata.

• Requires directly the output of BioGeoBEARS::runBSM() instead of separated cladogenetic and anagenetic event tables.

• Resolve paths to local tree and tip range files.

• Restore the S4 class 'BioGeoBEARS_model' if needed.

• Update the documentation.

Value

The convert_BSM_to_simmap() function returns a list with two elements:

• ⁠$simmap⁠ A unique simmap for a given biogeographic simulation as an object of classes c("simmap", "phylo"). This is a modified {ape} tree with additional elements to report range mapping, model parameters and likelihood.

  • ⁠$maps⁠ A list of named numeric vectors. Provides the mapping of ranges along each remaining edge. Names are the ranges. Values are residence times in each state across segments

  • ⁠$mapped.edge⁠ A numeric matrix. Provides the evolutionary time spent across ranges (columns) along the edges (rows). row.names() are the node ID at the rootward and tipward ends of each edge.

  • ⁠$Q⁠ Numerical matrix. The transition rates across ranges calculated from the ML parameter estimates of the model.

  • ⁠$logL⁠ Numeric. The log-likelihood of the data under the ML model.

• ⁠$residence_times⁠ Data.frame with two rows. Summarizes the residence time spent in each range along all branches, in (raw) evolutionary time (i.e., branch lengths), and in percentage (perc).

The convert_BSMs_to_simmaps() function loops around the convert_BSM_to_simmap() function to aggregate all simmaps from all biogeographic simulations in a unique list of classes c("multiSimmap", "multiPhylo").

• Each element is the ⁠$simmap⁠ of a biogeographic simulation obtained with convert_BSM_to_simmap().

• ⁠$residence_times⁠ summary data.frames are not preserved.

Notes on BioGeoBEARS

The R package BioGeoBEARS is needed for this function to work with biogeographic data. Please install it manually from: https://github.com/nmatzke/BioGeoBEARS.

Notes on using the resulting simmap object in phytools (adapted from Nicholas J. Matzke)

The phytools functions, like phytools::countSimmap(), will only count the anagenetic events (range transitions occurring along branches) as they were written assuming purely anagenetic models.

It remains possible to extract cladogenetic events (range transitions occurring at speciation) by comparing the last-state-below-a-node with the descendant-pairs-above-a-node. However, it is recommended to use the built-in functions from BioGeoBEARS to summarize the biogeographic history based on the tables of cladogenetic and anagenetic events obtained from BioGeoBEARS::runBSM(). simmap objects should primarily be considered as a tool for visualization.

Associated functions in R package BioGeoBEARS:

• BioGeoBEARS::simulate_source_areas_ana_clado(): To select randomly a unique area source for transition from a multi-area state to a single area.

• BioGeoBEARS::get_dmat_times_from_res(): To generate matrices of range expansion from source area to destination area.

• BioGeoBEARS::count_ana_clado_events(): To count the number and type of events from BSM tables.

• BioGeoBEARS::hist_event_counts(): To plot histograms of event counts across BSM tables.

Please note carefully that area-to-area dispersal events are not identical with the state transitions. For example, a state can be a geographic range with multiple areas, but under the logic of DEC-type models, a range-expansion event like ABC->ABCD actually means that a dispersal happened from some specific area (A, B, or C) to the new area. BSMs track this area-to-area sourcing in their cladogenetic and anagenetic event tables, at least if BioGeoBEARS::simulate_source_areas_ana_clado() has been run on the output of BioGeoBEARS::runBSM().

Author(s)

Nicholas J. Matzke. Contact: [email protected]

Changes by Maël Doré (see Details)

References

For BioGeoBEARS: Matzke, Nicholas J. (2018). BioGeoBEARS: BioGeography with Bayesian (and likelihood) Evolutionary Analysis with R Scripts. version 1.1.1, published on GitHub on November 6, 2018. doi:10.5281/zenodo.1478250. Website: http://phylo.wikidot.com/biogeobears.

See Also

phytools::countSimmap() phytools::make.simmap() BioGeoBEARS::simulate_source_areas_ana_clado() BioGeoBEARS::get_dmat_times_from_res() BioGeoBEARS::count_ana_clado_events() BioGeoBEARS::hist_event_counts()

Examples

if (deepSTRAPP::is_dev_version())
{
 ## Run only if you have R package 'BioGeoBEARS' installed.
 # Please install it manually from: https://github.com/nmatzke/BioGeoBEARS")

 ## Load phylogeny and tip data
 library(phytools)
 data(eel.tree)

  # (May take several minutes to run)
 ## Load directly output of prepare_trait_data() run on biogeographic data
 data(eel_biogeo_data, package = "deepSTRAPP")

 ## Convert BSM output into a unique simmap, including residence times
 simmap_1 <- convert_BSM_to_simmap(model_fit = eel_biogeo_data$best_model_fit,
                                    phylo = eel.tree,
                                    BSM_output = eel_biogeo_data$BSM_output,
                                    sim_index = 1)
 # Explore output
 str(simmap_1, max.level = 1)
 # Print residence times in each range
 simmap_1$residence_times
 # Plot simmap
 plot(simmap_1$simmap)

 ## Convert BSM output into all simmaps in a multiSimmap/multiPhylo object
 all_simmaps <- convert_BSMs_to_simmaps(model_fit = eel_biogeo_data$best_model_fit,
                                        phylo = eel.tree,
                                        BSM_output = eel_biogeo_data$BSM_output)
 # Explore output
 str(all_simmaps, max.level = 1)
 # Plot simmap n°1
 plot(all_simmaps[[1]]) 
}

---

Convert a contsimmap object into a list of contMaps

Description

Convert a contsimmap object generated with contsimmap::make.contsimmap() into a list of contMaps as produced by the phytools::contMap() function.

Usage

convert_contsimmap_to_contMaps(contsimmap, verbose = FALSE)

Arguments

contsimmap	List of class "contsimmap", typically generated with contsimmap::make.contsimmap(), that summarizes a set of continuous stochastic maps. Each map represents a simulation of the evolution of a continuous trait compatible with the associated evolutionary model fit on observed trait data.
verbose	Logical. Whether to display conversion progress every 10 maps.

Details

A contsimmap object is typically produced by the contsimmap::make.contsimmap() function. It stores results of a continuous stochastic mapping simulation that represent independent evolutionary history of a continuous trait mapped along a time-calibrated phylogeny, and compatible with the associated evolutionary model fit on observed trait data. The method is described in Martin & Weber, 2026.

A typical "contsimmap" object stores simulated trait values across maps in a 3D array, recording values for all time-points used to represent the continuous evolution. While "contsimmap" objects can record evolution of multivariate traits, deepSTRAPP currently works only with simulation of univariate trait evolution. See the contsimmap::make.contsimmap() help section for a detailed description of the structure of the object.

A "contMap" object typically produced by the phytools::contMap() function also records the evolution of a continuous trait alongside branches of a time-calibrated phylogeny. Trait values are stored in ⁠$maps⁠ as list of named vectors recording trait values (as names scaled between 0 and 1000) and length of the associated branch segment between two time-points (as values). Therefore, a unique "contsimmap" object is converted into a list of "contMaps", with each item being a "contMap" that represents a unique evolutionary history simulation conditioned on the observed trait data and model fit.

Value

Returns a list of "contMap" objects. Each "contMap" represents a unique evolutionary history simulation conditioned on the observed trait data and model fit (i.e., a continuous stochastic map).

A "contMap" typically contains:

• ⁠$tree⁠ A list of classes "simmap" and "phylo". The mapped phylogeny including:

  • ⁠$maps⁠ A list of named numeric vectors. Provides the mapping of trait values along each edge (scaled between 0 and 1000).

  • ⁠$mapped.edge⁠ A numeric matrix. Provides the evolutionary time spent across trait values (columns) along the edges (rows).

• ⁠$cols⁠ A named vector mapping values to colors used to display trait evolution along the branches.

• ⁠$lims⁠ A numeric vector. Minimum and maximum trait values recorded during the simulation.

Author(s)

Maël Doré

References

For continuous stochastic mapping: Martin, B. S., & Weber, M. G. (2026). Stochastic character mapping of continuous traits on phylogenies. Systematic Biology, syag031. doi:10.1093/sysbio/syag031.

See Also

contsimmap::make.contsimmap() phytools::contMap() prepare_trait_data()

Examples

if (deepSTRAPP::is_dev_version())
{
 ## The R package 'contsimmap' is needed for this example to work.
 # Please install it manually from: https://github.com/bstaggmartin/contsimmap.

 # ----- Prepare data ----- #

 # Load eel phylogeny and tip data from the R package phytools
 # Source: Collar et al., 2014; DOI: 10.1038/ncomms6505

 data("eel.tree", package = "phytools")
 data("eel.data", package = "phytools")

 # Extract body size
 eel_tip_data <- stats::setNames(eel.data$Max_TL_cm,
                                 rownames(eel.data))

 # ----- Infer evolutionary rate ----- #

 # Use phytools::anc.ML() to extract the evolutionary rate under a BM
 BM_fit <- phytools::anc.ML(tree = eel.tree,
                            x = eel_tip_data,
                            model = "BM")
 # Extract evolutionary rate = sigma²
 sigma2 <- BM_fit$sig2

 # ----- Run continuous stochastic mapping ----- #

  # (May take several seconds to run)
 eel_contsimmap <- contsimmap::make.contsimmap(
    tree = eel.tree,
    trait.data = eel_tip_data,
    Xsig2 = list(sigma2), # Evolutionary rate
    nsim = 100, res = 100,
    verbose = TRUE)

 dim(eel_contsimmap) # 121 edges, 1 trait, 100 simulations

 # ----- Convert to a list of contMaps ----- #

 eel_contMaps <- convert_contsimmap_to_contMaps(
    contsimmap = eel_contsimmap, verbose = TRUE)

 # ----- Explore contMaps ----- #

 # Plot contMap from simulation n°1
 plot_contMap(eel_contMaps[[1]])
 # Plot contMap from simulation n°50
 plot_contMap(eel_contMaps[[50]])
 
}

---

Convert simmaps into densityMaps suitable for deepSTRAPP

Description

Convert simmaps representing independent simulations of evolutionary history into a densityMaps object suitable to use as input for a deepSTRAPP run.

simmaps are mapping discrete character/geographic evolutionary history as transitions in character states/geographic ranges across nodes and branches of a phylogeny.

densityMaps are recording the posterior probability/frequencies of being in a given state/range along branches. In deepSTRAPP format, densityMaps is a list of objects of class "densityMap" where each object corresponds to the frequency of presence/absence of a state/range.

Usage

convert_simmaps_to_densityMaps(
  simmaps,
  colors_per_levels = NULL,
  tol = 1e-05,
  verbose = TRUE
)

Arguments

simmaps	List of objects of class "simmap", typically generated with prepare_trait_data() or phytools::make.simmap(), that represent discrete character/geographic evolutionary history (i.e., transitions in character states/geographic ranges) mapped along branches.
colors_per_levels	Named character string. To set the colors to use to map each state/range posterior probabilities. Names = states/ranges; values = colors. If NULL (default), the rainbow() color scale will be used.
tol	Positive numeric. To set the tolerance used to match node ages and time steps (i.e., consider them equal). Default = 1e-5.
verbose	Logical. Whether to display progress every 100 edges. Default = TRUE.

Details

The function is a wrapper of the original phytools::densityMap() by Liam Revell. However, it does not produce a single densityMap for binary states, but rather can handle simmaps mapping any number of states/ranges, and produce a densityMaps list that summarizes the frequency of presence/absence of each state/range in subsequent densityMap objects.

Both simmaps and densityMaps objects can be used as inputs for a deepSTRAPP run with run_deepSTRAPP_for_focal_time() or run_deepSTRAPP_over_time().

simmaps retain the identity of each simulated history, allowing deepSTRAPP to keep track of which simulation generated each set of trait values (i.e., states/ranges). However, retaining all simulated histories can require substantial RAM, particularly for large phylogenies or many simulations.

Alternatively, densityMaps summarize the frequency of states/ranges across simulations. They require substantially less RAM and can be used to visualize the overall uncertainty in trait evolution with plot_densityMaps_overlay().

The trade-off is that densityMaps discard the identity of individual simulations and therefore cannot be used to track which simulated history generated a given set of trait values.

Value

The function returns a densityMaps as a list of objects of class "densityMap", where each object is mapping the frequency of presence/absence of a state/range.

The number of objects depends on the number of states/ranges observed across the simmaps provided for conversion. Each densityMap is named as "Density_map_X" with X being each of the state/range recorded across the simmaps.

Each densityMap is a list of three elements:

• ⁠$tree⁠ List of classes "simmap" and "phylo" that contains the phylogeny in ape format and the mapping of states/ranges frequencies along branches in ⁠$tree$maps⁠.

• ⁠$cols⁠ Named character strings. Colors mapped to the 0 to 1000 scale used to record frequencies in ⁠$tree$maps⁠.

• ⁠$states⁠ Character string with two values. First entry is the absence of the state/range recorded as "Not X". Second entry is the presence of the state/range recorded as "X", X being the state/range name.

Author(s)

Maël Doré

See Also

phytools::densityMap() plot_densityMaps_overlay() prepare_trait_data()

Examples

## Load data

# Load trait df
data(Ponerinae_trait_tip_data, package = "deepSTRAPP")
# Load phylogeny
data(Ponerinae_tree_old_calib, package = "deepSTRAPP")

# Extract categorical data with 3-levels
Ponerinae_cat_3lvl_tip_data <- setNames(object = Ponerinae_trait_tip_data$fake_cat_3lvl_tip_data,
                                         nm = Ponerinae_trait_tip_data$Taxa)
table(Ponerinae_cat_3lvl_tip_data)

# Select color scheme for states
colors_per_states <- c("forestgreen", "sienna", "goldenrod")
names(colors_per_states) <- c("arboreal", "subterranean", "terricolous")

 # (May take several minutes to run)
## Produce densityMaps using stochastic character mapping based on an ER Mk model
Ponerinae_cat_3lvl_data_old_calib <- prepare_trait_data(
    tip_data = Ponerinae_cat_3lvl_tip_data,
    phylo = Ponerinae_tree_old_calib,
    trait_data_type = "categorical",
    colors_per_levels = colors_per_states,
    evolutionary_models = "ER", # Use ER model
    run_stochastic_maps = TRUE,
    nb_simulations = 100, # Reduce number of simulations to save time
    seed = 1234, # Set seed for reproducibility
    return_simmaps = TRUE, # Return simmaps in the output
    plot_map = FALSE)

# Note that densityMaps are already produced by the [deepSTRAPP::prepare_trait_data()] function,
# but for the sake of example, we can convert the simmaps stored in the output into densityMaps,
# using [deepSTRAPP::convert_simmaps_to_densityMaps()].

## Convert simmaps to densityMaps as input format for deepSTRAPP
Ponerinae_densityMaps <- convert_simmaps_to_densityMaps(
   simmaps = Ponerinae_cat_3lvl_data_old_calib$simmaps)

# Plot densityMaps one by one
plot(Ponerinae_densityMaps[[1]]) # densityMap for state n°1 ("arboreal")
plot(Ponerinae_densityMaps[[2]]) # densityMap for state n°1 ("subterranean")
plot(Ponerinae_densityMaps[[3]]) # densityMap for state n°1 ("terricolous")

# Plot overlay of all densityMaps
plot_densityMaps_overlay(densityMaps = Ponerinae_densityMaps)

---

Cut the phylogeny and continuous trait mapping for a given focal time in the past

Description

Cuts off all the branches of the phylogeny which are younger than a specific time in the past (i.e. the focal_time). Branches overlapping the focal_time are shortened to the focal_time. Likewise, remove continuous trait mapping for the cut off branches by updating the ⁠$tree$maps⁠ and ⁠$tree$mapped.edge⁠ elements.

Usage

cut_contMap_for_focal_time(contMap, focal_time, keep_tip_labels = TRUE)

Arguments

contMap	Object of class "contMap", typically generated with phytools::contMap(), that contains a phylogenetic tree and associated continuous trait mapping. The phylogenetic tree must be rooted and fully resolved/dichotomous, but it does not need to be ultrametric (it can include fossils).
focal_time	Numeric. The time, in terms of time distance from the present, for which the tree and mapping must be cut. It must be smaller than the root age of the phylogeny.
keep_tip_labels	Logical. Specify whether terminal branches with a single descendant tip must retain their initial tip.label. Default is TRUE.

Details

The phylogenetic tree is cut for a specific time in the past (i.e. the focal_time).

When a branch with a single descendant tip is cut and keep_tip_labels = TRUE, the leaf left is labeled with the tip.label of the unique descendant tip.

When a branch with a single descendant tip is cut and keep_tip_labels = FALSE, the leaf left is labeled with the node ID of the unique descendant tip.

In all cases, when a branch with multiple descendant tips (i.e., a clade) is cut, the leaf left is labeled with the node ID of the MRCA of the cut-off clade.

The continuous trait mapping is updated accordingly by removing mapping associated with the cut off branches.

Value

The function returns the cut contMap as an object of class "contMap" as a list of three elements:

• ⁠$tree⁠. A list of classes "simmap" and "phylo". This function updates and adds multiple useful sub-elements to the ⁠$tree⁠ element.

  • ⁠$maps⁠ An updated list of named numeric vectors. Provides the mapping of trait values along each remaining edge.

  • ⁠$mapped.edge⁠ An updated matrix. Provides the evolutionary time spent across trait values (columns) along the remaining edges (rows).

  • ⁠$root_age⁠ Integer. Stores the age of the root of the tree.

  • ⁠$nodes_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_node_ID to the initial_node_ID. Each row is a node.

  • ⁠$initial_nodes_ID⁠ Character vector. Provides the initial ID of internal nodes. Used to plot internal node IDs as labels with ape::nodelabels().

  • ⁠$edges_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_edge_ID to the initial_edge_ID. Each row is an edge/branch.

  • ⁠$initial_edges_ID⁠ Character vector. Provides the initial ID of edges/branches. Used to plot edge/branch IDs as labels with ape::edgelabels().

• ⁠$cols⁠ A named character vector mapping trait values with colors to display on the plot.

• ⁠$lims⁠ A numeric vector of 2 storing the initial min and max values of the trait (trait values are scaled between 0 and 1000 in ⁠$tree$maps⁠)

Author(s)

Maël Doré

See Also

cut_phylo_for_focal_time() extract_most_likely_trait_values_for_focal_time()

For a guided tutorial, see this vignette: vignette("cut_phylogenies", package = "deepSTRAPP")

Examples

# ----- 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")

mammals_tree <- mammals$mammal.phy
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

# ----- Example 1: keep_tip_labels = TRUE ----- #

# 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, lwd = 2, fsize = c(0.5, 1))
ape::nodelabels(cex = 0.5)
abline(v = max(phytools::nodeHeights(mammals_contMap$tree)[,2]) - focal_time,
       col = "red", lty = 2, lwd = 2)

# Plot initial node labels on cut stochastic map
plot_contMap(updated_contMap, fsize = c(0.8, 1))
ape::nodelabels(cex = 0.8, text = updated_contMap$tree$initial_nodes_ID)

# ----- Example 2: keep_tip_labels = FALSE ----- #

# 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, fsize = c(0.5, 1))
ape::nodelabels(cex = 0.5)
abline(v = max(phytools::nodeHeights(mammals_contMap$tree)[,2]) - focal_time,
       col = "red", lty = 2, lwd = 2)

# Plot initial node labels on cut stochastic map
plot_contMap(updated_contMap, fsize = c(0.8, 1))
ape::nodelabels(cex = 0.8, text = updated_contMap$tree$initial_nodes_ID)

---

Cut the phylogenies and continuous trait mappings of a list of contMaps for a given focal time in the past

Description

Cuts off all the branches of the phylogeny which are younger than a specific time in the past (i.e. the focal_time). Branches overlapping the focal_time are shortened to the focal_time. Likewise, remove continuous trait mapping for the cut off branches by updating the ⁠$tree$maps⁠ and ⁠$tree$mapped.edge⁠ elements.

Usage

cut_contMaps_for_focal_time(contMaps, focal_time, keep_tip_labels = TRUE)

Arguments

contMaps	List of objects of class "contMap", typically generated with phytools::contMap(), that contains a phylogenetic tree and associated continuous trait mapping. The phylogenetic tree must be rooted and fully resolved/dichotomous, but it does not need to be ultrametric (it can include fossils).
focal_time	Numeric. The time, in terms of time distance from the present, for which the tree and mapping must be cut. It must be smaller than the root age of the phylogeny.
keep_tip_labels	Logical. Specify whether terminal branches with a single descendant tip must retain their initial tip.label. Default is TRUE.

Details

The phylogenetic trees are cut for a specific time in the past (i.e. the focal_time).

When a branch with a single descendant tip is cut and keep_tip_labels = TRUE, the leaf left is labeled with the tip.label of the unique descendant tip.

When a branch with a single descendant tip is cut and keep_tip_labels = FALSE, the leaf left is labeled with the node ID of the unique descendant tip.

In all cases, when a branch with multiple descendant tips (i.e., a clade) is cut, the leaf left is labeled with the node ID of the MRCA of the cut-off clade.

The continuous trait mappings are updated accordingly by removing mapping associated with the cut off branches.

Value

The function returns an updated list of objects as cut contMaps of class "contMap".

Each contMap object contains three elements:

• ⁠$tree⁠. A list of classes "simmap" and "phylo". This function updates and adds multiple useful sub-elements to the ⁠$tree⁠ element.

  • ⁠$maps⁠ An updated list of named numeric vectors. Provides the mapping of trait values along each remaining edge.

  • ⁠$mapped.edge⁠ An updated matrix. Provides the evolutionary time spent across trait values (columns) along the remaining edges (rows).

  • ⁠$root_age⁠ Integer. Stores the age of the root of the tree.

  • ⁠$nodes_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_node_ID to the initial_node_ID. Each row is a node.

  • ⁠$initial_nodes_ID⁠ Character vector. Provides the initial ID of internal nodes. Used to plot internal node IDs as labels with ape::nodelabels().

  • ⁠$edges_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_edge_ID to the initial_edge_ID. Each row is an edge/branch.

  • ⁠$initial_edges_ID⁠ Character vector. Provides the initial ID of edges/branches. Used to plot edge/branch IDs as labels with ape::edgelabels().

• ⁠$cols⁠ A named character vector mapping trait values with colors to display on the plot.

• ⁠$lims⁠ A numeric vector of length 2 storing the initial min and max values of the trait (trait values are scaled between 0 and 1000 in ⁠$tree$maps⁠)

Author(s)

Maël Doré

See Also

cut_phylo_for_focal_time() extract_most_likely_trait_values_for_focal_time() extract_all_trait_values_for_focal_time()

For a guided tutorial, see this vignette: vignette("cut_phylogenies", package = "deepSTRAPP")

Examples

if (deepSTRAPP::is_dev_version())
{
 ## The R package 'contsimmap' is needed for this example to work.
 # Please install it manually from: https://github.com/bstaggmartin/contsimmap.

 # ----- Prepare data ----- #

 # Load phylogeny and tip data
 # Load eel phylogeny and tip data from the R package phytools
 # Source: Collar et al., 2014; DOI: 10.1038/ncomms6505
 data("eel.tree", package = "phytools")
 data("eel.data", package = "phytools")

 # Extract body size
 eel_data <- stats::setNames(eel.data$Max_TL_cm,
                             rownames(eel.data))

  # (May take several seconds to run)
 ## Map multiple trait evolution histories on the phylogeny using continuous stochastic mapping
 eel_mapped_trait_data <- prepare_trait_data(
   tip_data = eel_data,
   trait_data_type = "continuous",
   phylo = eel.tree,
   evolutionary_models = "BM",
   run_stochastic_maps = TRUE,
   nb_simulations = 10,
   plot_map = FALSE,
   verbose = TRUE)

 # Extract list of stochastic maps
 eel_contMaps <- eel_mapped_trait_data$contMaps

 # ----- Example 1: keep_tip_labels = TRUE ----- #

 # Cut contMaps to 20 Mya while keeping tip.label
 # on terminal branches with a unique descending tip.
 updated_contMaps <- cut_contMaps_for_focal_time(contMaps = eel_contMaps,
                                                 focal_time = 20,
                                                 keep_tip_labels = TRUE)

 # Plot node labels on stochastic map n°1 with cut-off
 plot_contMap(eel_contMaps[[1]], lwd = 2, fsize = c(0.5, 1))
 ape::nodelabels(cex = 0.5)
 abline(v = max(phytools::nodeHeights(eel_contMaps[[1]]$tree)[,2]) - 20,
        col = "red", lty = 2, lwd = 2)

 # Plot initial node labels on cut stochastic map n°1
 plot_contMap(updated_contMaps[[1]], fsize = c(0.8, 1))
 ape::nodelabels(cex = 0.8, text = updated_contMaps[[1]]$tree$initial_nodes_ID)

 # ----- Example 2: keep_tip_labels = FALSE ----- #

 # Cut contMap to 20 Mya while NOT keeping tip.label.
 updated_contMaps <- cut_contMaps_for_focal_time(contMaps = eel_contMaps,
                                                 focal_time = 20,
                                                 keep_tip_labels = FALSE)

 # Plot node labels on stochastic map n°1 with cut-off
 plot_contMap(eel_contMaps[[1]], fsize = c(0.5, 1))
 ape::nodelabels(cex = 0.5)
 abline(v = max(phytools::nodeHeights(eel_contMaps[[1]]$tree)[,2]) - 20,
        col = "red", lty = 2, lwd = 2)

 # Plot initial node labels on cut stochastic map
 plot_contMap(updated_contMaps[[1]], fsize = c(0.8, 1))
 ape::nodelabels(cex = 0.8, text = updated_contMaps[[1]]$tree$initial_nodes_ID)
 
}

---

Cut the phylogeny and posterior probability mapping of a categorical trait for a given focal time in the past

Description

Cuts off all the branches of the phylogeny which are younger than a specific time in the past (i.e. the focal_time). Branches overlapping the focal_time are shortened to the focal_time. Likewise, remove posterior probability mapping of the categorical trait for the cut off branches by updating the ⁠$tree$maps⁠ and ⁠$tree$mapped.edge⁠ elements.

Usage

cut_densityMap_for_focal_time(densityMap, focal_time, keep_tip_labels = TRUE)

Arguments

densityMap	Object of class "densityMap", typically generated with phytools::densityMap(), that contains a phylogenetic tree and associated posterior probability mapping of a categorical trait. The phylogenetic tree must be rooted and fully resolved/dichotomous, but it does not need to be ultrametric (it can include fossils).
focal_time	Numeric. The time, in terms of time distance from the present, for which the tree and mapping must be cut. It must be smaller than the root age of the phylogeny.
keep_tip_labels	Logical. Specify whether terminal branches with a single descendant tip must retain their initial tip.label. Default is TRUE.

Details

The phylogenetic tree is cut for a specific time in the past (i.e. the focal_time).

When a branch with a single descendant tip is cut and keep_tip_labels = TRUE, the leaf left is labeled with the tip.label of the unique descendant tip.

When a branch with a single descendant tip is cut and keep_tip_labels = FALSE, the leaf left is labeled with the node ID of the unique descendant tip.

In all cases, when a branch with multiple descendant tips (i.e., a clade) is cut, the leaf left is labeled with the node ID of the MRCA of the cut-off clade.

The posterior probability mapping of a categorical trait is updated accordingly by removing mapping associated with the cut off branches.

Value

The function returns the cut densityMap as an object of class "densityMap" with three elements.

It contains a ⁠$tree⁠ element of classes "simmap" and "phylo". This function updates and adds multiple useful sub-elements to the ⁠$tree⁠ element. * ⁠$maps⁠ An updated list of named numeric vectors. Provides the mapping of posterior probability of the state along each remaining edge. * ⁠$mapped.edge⁠ An updated matrix. Provides the evolutionary time spent across posterior probabilities (columns) along the remaining edges (rows). * ⁠$root_age⁠ Integer. Stores the age of the root of the tree. * ⁠$nodes_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_node_ID to the initial_node_ID. Each row is a node. * ⁠$initial_nodes_ID⁠ Character vector. Provides the initial ID of internal nodes. Used to plot internal node IDs as labels with ape::nodelabels(). * ⁠$edges_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_edge_ID to the initial_edge_ID. Each row is an edge/branch. * ⁠$initial_edges_ID⁠ Character vector. Provides the initial ID of edges/branches. Used to plot edge/branch IDs as labels with ape::edgelabels().

The ⁠$col⁠ element describes the colors used to map each possible posterior probability value from 0 to 1000.

The ⁠$states⁠ element provides the name of the states. Here, the first value is the absence of the state labeled as "Not X" with X being the state. The second value is the name of the state.

High posterior probability reflects high likelihood to harbor the state. Low probability reflects high likelihood to NOT harbor the state.

Author(s)

Maël Doré

See Also

cut_phylo_for_focal_time() extract_most_likely_trait_values_for_focal_time() extract_all_trait_values_for_focal_time()

For a guided tutorial, see this vignette: vignette("cut_phylogenies", package = "deepSTRAPP")

Examples

# ----- Prepare data ----- #

# Load mammals phylogeny and data from the R package motmot included within deepSTRAPP
# Data source: Slater, 2013; DOI: 10.1111/2041-210X.12084
data("mammals", package = "deepSTRAPP")

# Obtain mammal tree
mammals_tree <- mammals$mammal.phy
# 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)

 # (May take several minutes to run)
# 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,
                                       plot_map = 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]]

# ----- Example 1: keep_tip_labels = TRUE ----- #

# 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, fsize = 0.7, lwd = 2)
ape::nodelabels(cex = 0.7)
abline(v = max(phytools::nodeHeights(mammals_densityMap_small$tree)[,2]) - focal_time,
       col = "red", lty = 2, lwd = 2)

# Plot initial node labels on cut stochastic map
phytools::plot.densityMap(updated_mammals_densityMap_small, fsize = 0.8)
ape::nodelabels(cex = 0.8, text = updated_mammals_densityMap_small$tree$initial_nodes_ID)

# ----- Example 2: keep_tip_labels = FALSE ----- #

# Cut densityMap to 80 Mya while NOT keeping tip.label
updated_mammals_densityMap_small <- cut_densityMap_for_focal_time(
     densityMap = mammals_densityMap_small,
     focal_time = focal_time,
     keep_tip_labels = FALSE)

# Plot node labels on initial stochastic map with cut-off
phytools::plot.densityMap(mammals_densityMap_small, fsize = 0.7, lwd = 2)
ape::nodelabels(cex = 0.7)
abline(v = max(phytools::nodeHeights(mammals_densityMap_small$tree)[,2]) - focal_time,
       col = "red", lty = 2, lwd = 2)

# Plot initial node labels on cut stochastic map
phytools::plot.densityMap(updated_mammals_densityMap_small, fsize = 0.8)
ape::nodelabels(cex = 0.8, text = updated_mammals_densityMap_small$tree$initial_nodes_ID)

---

Cut phylogenies and posterior probability mapping of each state for a given focal time in the past

Description

Cuts off all the branches of the phylogeny which are younger than a specific time in the past (i.e. the focal_time). Branches overlapping the focal_time are shortened to the focal_time. Likewise, remove posterior probability mapping of the categorical trait for the cut off branches by updating the ⁠$tree$maps⁠ and ⁠$tree$mapped.edge⁠ elements.

Usage

cut_densityMaps_for_focal_time(densityMaps, focal_time, keep_tip_labels = TRUE)

Arguments

densityMaps	List of objects of class "densityMap", typically generated with prepare_trait_data(). Each densityMap (see phytools::densityMap()) contains a phylogenetic tree and associated posterior probability mapping of a categorical trait. The phylogenetic tree must be rooted and fully resolved/dichotomous, but it does not need to be ultrametric (it can include fossils).
focal_time	Numeric. The time, in terms of time distance from the present, for which the tree and mapping must be cut. It must be smaller than the root age of the phylogeny.
keep_tip_labels	Logical. Specify whether terminal branches with a single descendant tip must retain their initial tip.label. Default is TRUE.

Details

The phylogenetic tree is cut for a specific time in the past (i.e. the focal_time).

When a branch with a single descendant tip is cut and keep_tip_labels = TRUE, the leaf left is labeled with the tip.label of the unique descendant tip.

When a branch with a single descendant tip is cut and keep_tip_labels = FALSE, the leaf left is labeled with the node ID of the unique descendant tip.

In all cases, when a branch with multiple descendant tips (i.e., a clade) is cut, the leaf left is labeled with the node ID of the MRCA of the cut-off clade.

The continuous trait mapping is updated accordingly by removing mapping associated with the cut off branches.

Value

The function returns an updated list of objects as cut densityMaps of class "densityMap".

Each densityMap object contains three elements:

• ⁠$tree⁠ element of classes "simmap" and "phylo". This function updates and adds multiple useful sub-elements to the ⁠$tree⁠ element.

  • ⁠$maps⁠ An updated list of named numeric vectors. Provides the mapping of posterior probability of the state along each remaining edge.

  • ⁠$mapped.edge⁠ An updated matrix. Provides the evolutionary time spent across posterior probabilities (columns) along the remaining edges (rows).

  • ⁠$root_age⁠ Integer. Stores the age of the root of the tree.

  • ⁠$nodes_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_node_ID to the initial_node_ID. Each row is a node.

  • ⁠$initial_nodes_ID⁠ Character vector. Provides the initial ID of internal nodes. Used to plot internal node IDs as labels with ape::nodelabels().

  • ⁠$edges_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_edge_ID to the initial_edge_ID. Each row is an edge/branch.

  • ⁠$initial_edges_ID⁠ Character vector. Provides the initial ID of edges/branches. Used to plot edge/branch IDs as labels with ape::edgelabels().

• ⁠$col⁠ element describes the colors used to map each possible posterior probability value from 0 to 1000.

• ⁠$states⁠ element provide the name of the states. Here, the first value is the absence of the state labeled as "Not X" with X being the state. The second value is the name of the state.

High posterior probability reflects high likelihood to harbor the state. Low probability reflects high likelihood to NOT harbor the state.

Author(s)

Maël Doré

See Also

cut_phylo_for_focal_time() cut_densityMap_for_focal_time() extract_most_likely_trait_values_for_focal_time() extract_all_trait_values_for_focal_time()

Examples

# ----- Prepare data ----- #

# Load mammals phylogeny and data from the R package motmot, and implemented in deepSTRAPP
# Data source: Slater, 2013; DOI: 10.1111/2041-210X.12084
data("mammals", package = "deepSTRAPP")

# Obtain mammal tree
mammals_tree <- mammals$mammal.phy
# 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)

 # (May take several minutes to run)
# 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,
                                       plot_map = FALSE)

# Set focal time
focal_time <- 80

# Extract the density maps
mammals_densityMaps <- mammals_cat_data$densityMaps

# ----- Example 1: keep_tip_labels = TRUE ----- #

# Cut 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_densityMaps,
    focal_time = focal_time,
    keep_tip_labels = TRUE)

## Plot density maps as overlay of all state posterior probabilities
# ?plot_densityMaps_overlay

# Plot initial density maps
plot_densityMaps_overlay(densityMaps = mammals_densityMaps, fsize = 0.5)
abline(v = max(phytools::nodeHeights(mammals_densityMaps[[1]]$tree)[,2]) - focal_time,
       col = "red", lty = 2, lwd = 2)

# Plot updated/cut density maps
plot_densityMaps_overlay(densityMaps = updated_mammals_densityMaps, fsize = 0.8)

# ----- Example 2: keep_tip_labels = FALSE ----- #

# Cut densityMap to 80 Mya while NOT keeping tip.label
updated_mammals_densityMaps <- cut_densityMaps_for_focal_time(
    densityMaps = mammals_densityMaps,
    focal_time = focal_time,
    keep_tip_labels = FALSE)

# Plot initial density maps
plot_densityMaps_overlay(densityMaps = mammals_densityMaps, fsize = 0.5)
abline(v = max(phytools::nodeHeights(mammals_densityMaps[[1]]$tree)[,2]) - focal_time,
       col = "red", lty = 2, lwd = 2)

# Plot updated/cut density maps
plot_densityMaps_overlay(densityMaps = updated_mammals_densityMaps, fsize = 0.8)

---

Cut the phylogeny for a given time in the past

Description

Cuts off all the branches of the phylogeny which are younger than a specific time in the past (i.e. the focal_time). Branches overlapping the focal_time are shortened to the focal_time.

Usage

cut_phylo_for_focal_time(tree, focal_time, keep_tip_labels = TRUE)

Arguments

tree	Object of class "phylo". The phylogenetic tree must be rooted and fully resolved/dichotomous, but it does not need to be ultrametric (it can include fossils).
focal_time	Numeric. The time, in terms of time distance from the present, for which the tree must be cut. It must be smaller than the root age of the phylogeny.
keep_tip_labels	Logical. Specify whether terminal branches with a single descendant tip must retain their initial tip.label. Default is TRUE.

Details

When a branch with a single descendant tip is cut and keep_tip_labels = TRUE, the leaf left is labeled with the tip.label of the unique descendant tip.

When a branch with a single descendant tip is cut and keep_tip_labels = FALSE, the leaf left is labeled with the node ID of the unique descendant tip.

In all cases, when a branch with multiple descendant tips (i.e., a clade) is cut, the leaf left is labeled with the node ID of the MRCA of the cut-off clade.

Value

The function returns the cut phylogeny as an object of class "phylo". It adds multiple useful elements to the object.

• ⁠$root_age⁠ Integer. Stores the age of the root of the tree.

• ⁠$nodes_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_node_ID to the initial_node_ID. Each row is a node.

• ⁠$initial_nodes_ID⁠ Character vector. Provides the initial ID of internal nodes. Used to plot internal node IDs as labels with ape::nodelabels().

• ⁠$edges_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_edge_ID to the initial_edge_ID. Each row is an edge/branch.

• ⁠$initial_edges_ID⁠ Character vector. Provides the initial ID of edges/branches. Used to plot edge/branch IDs as labels with ape::edgelabels().

Author(s)

Maël Doré

See Also

cut_contMap_for_focal_time() cut_densityMaps_for_focal_time()

For a guided tutorial, see this vignette: vignette("cut_phylogenies", package = "deepSTRAPP")

Examples

# Load eel phylogeny from the R package phytools
# Source: Collar et al., 2014; DOI: 10.1038/ncomms6505
data("eel.tree", package = "phytools")

# ----- Example 1: keep_tip_labels = TRUE ----- #

# Cut tree to 30 Mya while keeping tip.label on terminal branches with a unique descending tip.
cut_eel.tree <- cut_phylo_for_focal_time(tree = eel.tree, focal_time = 30, keep_tip_labels = TRUE)

# Plot internal node labels on initial tree with cut-off
plot(eel.tree, cex = 0.5)
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_eel.tree$initial_nodes_ID)] <- NA
ape::nodelabels(text = nodelabels_in_cut_tree)

# Plot initial internal node labels on cut tree
plot(cut_eel.tree, cex = 0.8)
ape::nodelabels(text = cut_eel.tree$initial_nodes_ID)

# Plot edge labels on initial tree with cut-off
plot(eel.tree, cex = 0.5)
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_eel.tree$initial_edges_ID)] <- NA
ape::edgelabels(text = edgelabels_in_cut_tree)

# Plot initial edge labels on cut tree
plot(cut_eel.tree, cex = 0.8)
ape::edgelabels(text = cut_eel.tree$initial_edges_ID)

# ----- Example 2: keep_tip_labels = FALSE ----- #

# Cut tree to 30 Mya without keeping tip.label on terminal branches with a unique descending tip.
# All tip.labels are converted to their descending/tipward node ID
cut_eel.tree <- cut_phylo_for_focal_time(tree = eel.tree, focal_time = 30, keep_tip_labels = FALSE)
plot(cut_eel.tree, cex = 0.8)

---

Cut the phylogeny and categorical trait/range mapping for a given focal time in the past

Description

Cuts off all the branches of the phylogeny which are younger than a specific time in the past (i.e. the focal_time). Branches overlapping the focal_time are shortened to the focal_time. Likewise, remove mapping for the cut off branches by updating the ⁠$maps⁠ and ⁠$mapped.edge⁠ elements.

Usage

cut_simmap_for_focal_time(simmap, focal_time, keep_tip_labels = TRUE)

Arguments

simmap	Object with the classes "phylo" and "simmap", typically generated with phytools::make.simmap(), that contains a phylogenetic tree and associated categorical trait/range mapping. The phylogenetic tree must be rooted and fully resolved/dichotomous, but it does not need to be ultrametric (it can include fossils).
focal_time	Numeric. The time, in terms of time distance from the present, for which the tree and mapping must be cut. It must be smaller than the root age of the phylogeny.
keep_tip_labels	Logical. Specify whether terminal branches with a single descendant tip must retain their initial tip.label. Default is TRUE.

Details

The mapped phylogenetic tree is cut for a specific time in the past (i.e. the focal_time).

When a branch with a single descendant tip is cut and keep_tip_labels = TRUE, the leaf left is labeled with the tip.label of the unique descendant tip.

When a branch with a single descendant tip is cut and keep_tip_labels = FALSE, the leaf left is labeled with the node ID of the unique descendant tip.

In all cases, when a branch with multiple descendant tips (i.e., a clade) is cut, the leaf left is labeled with the node ID of the MRCA of the cut-off clade.

The stochastic trait mapping is updated accordingly by removing mapping associated with the cut off branches.

Value

The function returns the cut/updated simmap as a list of at least six elements with the classes "phylo" and "simmap":

Initial updated elements:

• ⁠$edge⁠ A numeric matrix with two columns listing the updated ID of rootward and tipward nodes of the remaining edges

• ⁠$edge.length⁠ A numeric vector. Updated length of remaining edges.

• ⁠$Nnode⁠ Integer. Number of remaining internal nodes.

• ⁠$tip.label⁠ Character vector. Labels of the branches cut at the focal_time.

  • If keep_tip_labels = TRUE, the cut terminal branches are labeled with the tip.label of their unique descendant tip.

  • If keep_tip_labels = FALSE, the cut terminal branches are labeled with the node ID of their unique descendant tip.

• ⁠$maps⁠ An updated list of named numeric vectors. Provides the mapping of trait values along each remaining edge.

• ⁠$mapped.edge⁠ An updated matrix. Provides the evolutionary time spent across trait values (columns) along the remaining edges (rows).

Additional elements:

• ⁠$root_age⁠ Integer. Stores the age of the root of the tree.

• ⁠$nodes_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_node_ID to the initial_node_ID. Each row is a node.

• ⁠$initial_nodes_ID⁠ Character vector. Provides the initial ID of internal nodes. Used to plot internal node IDs as labels with ape::nodelabels().

• ⁠$edges_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_edge_ID to the initial_edge_ID. Each row is an edge/branch.

• ⁠$initial_edges_ID⁠ Character vector. Provides the initial ID of edges/branches. Used to plot edge/branch IDs as labels with ape::edgelabels().

Author(s)

Maël Doré

See Also

cut_phylo_for_focal_time() extract_most_likely_trait_values_for_focal_time() extract_all_trait_values_for_focal_time()

For a guided tutorial, see this vignette: vignette("cut_phylogenies", package = "deepSTRAPP")

Examples

# ----- Prepare data ----- #

## Load mammals phylogeny and data from the R package motmot, and implemented in deepSTRAPP
# Data source: Slater, 2013; DOI: 10.1111/2041-210X.12084
data("mammals", package = "deepSTRAPP")

# Obtain mammal tree
mammals_tree <- mammals$mammal.phy
# 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)
 # (May take several seconds to run)
# 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,
                                       return_simmaps = TRUE,
                                       plot_map = FALSE)

# Extract the simulated evolutionary histories (simmaps)
mammals_simmaps <- mammals_cat_data$simmaps
# Extract simulation n°1
mammals_simmap_1 <- mammals_simmaps[[1]]

# Set focal time to 80Mya
focal_time <- 80

# ----- Example 1: keep_tip_labels = TRUE ----- #

# Cut the simmap to 80 Mya while keeping tip.label
# on terminal branches with a unique descending tip.
updated_simmap_1 <- cut_simmap_for_focal_time(simmap = mammals_simmap_1,
                                              focal_time = focal_time,
                                              keep_tip_labels = TRUE)

# Plot node labels on initial stochastic map with cut-off
plot(mammals_simmap_1, lwd = 2, fsize = 0.5,
     colors = c(large = "black", medium = "blue", small = "dodgerblue"))
ape::nodelabels(cex = 0.5)
abline(v = max(phytools::nodeHeights(mammals_simmap_1)[,2]) - focal_time,
       col = "red", lty = 2, lwd = 2)

# Plot initial node labels on cut stochastic map
plot(updated_simmap_1, fsize = 0.8,
     colors = c(large = "black", medium = "blue", small = "dodgerblue"))
ape::nodelabels(cex = 0.8, text = updated_simmap_1$initial_nodes_ID)

# ----- Example 2: keep_tip_labels = FALSE ----- #

# Cut the simmap to 80 Mya while NOT keeping tip.label.
updated_simmap_1 <- cut_simmap_for_focal_time(simmap = mammals_simmap_1,
                                              focal_time = focal_time,
                                              keep_tip_labels = FALSE)

# Plot node labels on initial stochastic map with cut-off
plot(mammals_simmap_1, lwd = 2, fsize = 0.5,
     colors = c(large = "black", medium = "blue", small = "dodgerblue"))
ape::nodelabels(cex = 0.5)
abline(v = max(phytools::nodeHeights(mammals_simmap_1)[,2]) - focal_time,
       col = "red", lty = 2, lwd = 2)

# Plot initial node labels on cut stochastic map
plot(updated_simmap_1, fsize = 0.8,
     colors = c(large = "black", medium = "blue", small = "dodgerblue"))
ape::nodelabels(cex = 0.8, text = updated_simmap_1$initial_nodes_ID)

---

Cut the phylogenies and categorical trait/range mappings of a list of simmaps for a given focal time in the past

Description

Cuts off all the branches of the phylogeny which are younger than a specific time in the past (i.e. the focal_time). Branches overlapping the focal_time are shortened to the focal_time. Likewise, remove continuous trait mapping for the cut off branches by updating the ⁠$maps⁠ and ⁠$mapped.edge⁠ elements.

Usage

cut_simmaps_for_focal_time(simmaps, focal_time, keep_tip_labels = TRUE)

Arguments

simmaps	List of objects with the classes "phylo" and "simmap", typically generated with phytools::make.simmap(), that contains a phylogenetic tree and associated categorical trait/range mapping. The phylogenetic tree must be rooted and fully resolved/dichotomous, but it does not need to be ultrametric (it can include fossils).
focal_time	Numeric. The time, in terms of time distance from the present, for which the tree and mapping must be cut. It must be smaller than the root age of the phylogeny.
keep_tip_labels	Logical. Specify whether terminal branches with a single descendant tip must retain their initial tip.label. Default is TRUE.

Details

The phylogenetic trees are cut for a specific time in the past (i.e. the focal_time).

When a branch with a single descendant tip is cut and keep_tip_labels = TRUE, the leaf left is labeled with the tip.label of the unique descendant tip.

When a branch with a single descendant tip is cut and keep_tip_labels = FALSE, the leaf left is labeled with the node ID of the unique descendant tip.

In all cases, when a branch with multiple descendant tips (i.e., a clade) is cut, the leaf left is labeled with the node ID of the MRCA of the cut-off clade.

The categorical trait/range mappings are updated accordingly by removing mapping associated with the cut off branches.

Value

The function returns an updated list of objects as cut/updated simmaps of classes "phylo" and "simmap".

Each simmap object represents an updated stochastic map cut/updated simmap as a list of at least six elements:

Initial updated elements:

• ⁠$edge⁠ A numeric matrix with two columns listing the updated ID of rootward and tipward nodes of the remaining edges

• ⁠$edge.length⁠ A numeric vector. Updated length of remaining edges.

• ⁠$Nnode⁠ Integer. Number of remaining internal nodes.

• ⁠$tip.label⁠ Character vector. Labels of the branches cut at the focal_time.

  • If keep_tip_labels = TRUE, the cut terminal branches are labeled with the tip.label of their unique descendant tip.

  • If keep_tip_labels = FALSE, the cut terminal branches are labeled with the node ID of their unique descendant tip.

• ⁠$maps⁠ An updated list of named numeric vectors. Provides the mapping of trait values along each remaining edge.

• ⁠$mapped.edge⁠ An updated matrix. Provides the evolutionary time spent across trait values (columns) along the remaining edges (rows).

Additional elements:

• ⁠$root_age⁠ Integer. Stores the age of the root of the tree.

• ⁠$nodes_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_node_ID to the initial_node_ID. Each row is a node.

• ⁠$initial_nodes_ID⁠ Character vector. Provides the initial ID of internal nodes. Used to plot internal node IDs as labels with ape::nodelabels().

• ⁠$edges_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_edge_ID to the initial_edge_ID. Each row is an edge/branch.

• ⁠$initial_edges_ID⁠ Character vector. Provides the initial ID of edges/branches. Used to plot edge/branch IDs as labels with ape::edgelabels().

Author(s)

Maël Doré

See Also

cut_phylo_for_focal_time() extract_most_likely_trait_values_for_focal_time() extract_all_trait_values_for_focal_time()

For a guided tutorial, see this vignette: vignette("cut_phylogenies", package = "deepSTRAPP")

Examples

# ----- Prepare data ----- #

## Load mammals phylogeny and data from the R package motmot, and implemented in deepSTRAPP
# Data source: Slater, 2013; DOI: 10.1111/2041-210X.12084
data("mammals", package = "deepSTRAPP")

# Obtain mammal tree
mammals_tree <- mammals$mammal.phy
# 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)

 # (May take several seconds to run)
# 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,
                                       return_simmaps = TRUE,
                                       plot_map = FALSE)

# Extract the simulated evolutionary histories (simmaps)
mammals_simmaps <- mammals_cat_data$simmaps

# Set focal time to 80Mya
focal_time <- 80

# ----- Example 1: keep_tip_labels = TRUE ----- #

# Cut the simmap to 80 Mya while keeping tip.label
# on terminal branches with a unique descending tip.
updated_simmaps <- cut_simmaps_for_focal_time(simmaps = mammals_simmaps,
                                              focal_time = focal_time,
                                              keep_tip_labels = TRUE)

# Plot node labels on initial stochastic map n°1 with cut-off
plot(mammals_simmaps[[1]], lwd = 2, fsize = 0.5,
     colors = c(large = "black", medium = "blue", small = "dodgerblue"))
ape::nodelabels(cex = 0.5)
abline(v = max(phytools::nodeHeights(mammals_simmaps[[1]])[,2]) - focal_time,
       col = "red", lty = 2, lwd = 2)

# Plot initial node labels on cut stochastic map n°1
plot(updated_simmaps[[1]], fsize = 0.8,
     colors = c(large = "black", medium = "blue", small = "dodgerblue"))
ape::nodelabels(cex = 0.8, text = updated_simmaps[[1]]$initial_nodes_ID)

# ----- Example 2: keep_tip_labels = FALSE ----- #

# Cut the simmap to 80 Mya while NOT keeping tip.label.
updated_simmaps <- cut_simmaps_for_focal_time(simmaps = mammals_simmaps,
                                              focal_time = focal_time,
                                              keep_tip_labels = FALSE)

# Plot node labels on initial stochastic map n°1 with cut-off
plot(mammals_simmaps[[1]], lwd = 2, fsize = 0.5,
     colors = c(large = "black", medium = "blue", small = "dodgerblue"))
ape::nodelabels(cex = 0.5)
abline(v = max(phytools::nodeHeights(mammals_simmaps[[1]])[,2]) - focal_time,
       col = "red", lty = 2, lwd = 2)

# Plot initial node labels on cut stochastic map n°1
plot(updated_simmaps[[1]], fsize = 0.8,
     colors = c(large = "black", medium = "blue", small = "dodgerblue"))
ape::nodelabels(cex = 0.8, text = updated_simmaps[[1]]$initial_nodes_ID)

---

Data summarizing the evolution of geographic ranges in eels

Description

A list containing (fake) geographic ranges data of eels mapped on the phylogeny, modeled with R package BioGeoBEARS. This object was obtained with prepare_trait_data(). Initial data based on feeding habits was altered to be transformed into ranges "A" and "B", and then adding arbitrarily multi-area "AB" ranges. This is NOT real biogeographic data. Please refer to the initial article for real data.

Original data source: Collar, D. C., P. C. Wainwright, M. E. Alfaro, L. J. Revell, and R. S. Mehta (2014) Biting disrupts integration to spur skull evolution in eels. Nature Communications, 5, 5505. doi:10.1038/ncomms6505

Usage

data(eel_biogeo_data)

Format

A list with 9 elements.

Details

A list of 9 elements containing information on the evolution of geographic ranges in eels. This object was obtained with prepare_trait_data().

• ⁠$densityMaps⁠ List of objects of class "densityMap" that contains a phylogenetic tree and associated mapping of probability to harbor a given range along branches. The list contains only a "densityMap" per unique area because split_multi_area_ranges was set to TRUE.

• ⁠$densityMaps_all_ranges⁠ List of objects of class "densityMap" that contains a phylogenetic tree and associated mapping of probability to harbor a given range along branches. The list contains one "densityMap" per range found along branches during the simulated biogeographic histories.

• ⁠$trait_data_type⁠ Character string. Record the type of trait data. Here: "biogeographic".

• ⁠$ace⁠ Numeric matrix. Record the posterior probabilities of ancestral ranges estimated at internal nodes. Only unique areas are considered among the ranges. Multi-area ranges have been split among unique ranges. Rows are internal nodes. Columns are ranges. Values are posterior probabilities of each range per node.

• ⁠$ace_all_ranges⁠ Numeric matrix. Record the posterior probabilities of ancestral ranges estimated at internal nodes. All ranges observed along branches during the simulated biogeographic histories are present. Rows are internal nodes. Columns are ranges. Values are posterior probabilities of each range per node.

• ⁠$BSM_output⁠ List of two lists that contain summary information on cladogenetic (⁠$RES_caldo_events_tables⁠) and anagenetic (⁠$RES_ana_events_tables⁠) events recorded across the 1000 simulations of biogeographic histories performed during Biogeographic Stochastic Mapping (BSM). Each element of the list is a data.frame recording events occurring during one simulation.

• ⁠$simmaps⁠ List of 1000 objects of class "simmap". Each simmap object is a phylogeny with one simulated biogeographic history (i.e., transitions in geographic ranges) mapped along branches.

• ⁠$best_model_fit⁠ List that provides the output of the best fitting model (Here: DEC+J model).

• ⁠$model_selection_df⁠ Data.frame that summarizes model comparisons used to select the best fitting model.

References

Collar, D. C., P. C. Wainwright, M. E. Alfaro, L. J. Revell, and R. S. Mehta (2014) Biting disrupts integration to spur skull evolution in eels. Nature Communications, 5, 5505. doi:10.1038/ncomms6505

See Also

prepare_trait_data()

---

Data summarizing the evolution of feeding habits in eels using a 3-level factor as categorical trait

Description

A list containing feeding habits data of eels mapped on the phylogeny, modeled with geiger::fitDiscrete. This object was obtained with prepare_trait_data(). Initial data was altered arbitrarily to create three categories, adding a "kiss" feeding habit to the initial "bite" and "suction" data. This is NOT real biological data. Please refer to the initial article for real data.

Original data source: Collar, D. C., P. C. Wainwright, M. E. Alfaro, L. J. Revell, and R. S. Mehta (2014) Biting disrupts integration to spur skull evolution in eels. Nature Communications, 5, 5505. doi:10.1038/ncomms6505

Usage

data(eel_cat_3lvl_data)

Format

A list with 6 elements.

Details

A list of six elements containing information on the evolution of feeding habits in eels. This object was obtained with prepare_trait_data().

• ⁠$densityMaps⁠ List of objects of class "densityMap" that contains a phylogenetic tree and associated mapping of probability to harbor a given state/range along branches. The list contains one "densityMap" per state/range found in the tip_data.

• ⁠$trait_data_type⁠ Character string. Record the type of trait data. Here: "categorical".

• ⁠$simmaps⁠ List of 100 stochastic mapping simulations for trait evolution. Each element is a "simmap" object (see phytools::make.simmap) representing a possible evolutionary history that fits states observed on tips, inferred ACE at internal nodes, and transition rates as estimated from the best fit model.

• ⁠$ace⁠ Numeric matrix. Record the posterior probabilities of ancestral states/ranges (characters) estimates (ACE) at internal nodes. Rows are internal nodes. Columns are states/ranges. Values are posterior probabilities of each state per node.

• ⁠$best_model_fit⁠ List that provides the output of the best fitting model (Here: ER model).

• ⁠$model_selection_df⁠ Data.frame that summarizes model comparisons used to select the best fitting model.

References

Collar, D. C., P. C. Wainwright, M. E. Alfaro, L. J. Revell, and R. S. Mehta (2014) Biting disrupts integration to spur skull evolution in eels. Nature Communications, 5, 5505. doi:10.1038/ncomms6505

See Also

prepare_trait_data()

---

Extract all trait data from stochastic maps at a given time in the past

Description

Extracts all trait values/states/ranges found along branches of multiple stochastic maps at a specific time in the past (i.e. the focal_time). Optionally, the function can update the mapped phylogenies (contMaps or densityMaps) such as branches overlapping the focal_time are shortened to the focal_time, and the trait mapping for the cut off branches are removed by updating the ⁠$maps⁠ and ⁠$mapped.edge⁠ elements.

Usage

extract_all_trait_values_for_focal_time(
  contMaps = NULL,
  densityMaps = NULL,
  simmaps = NULL,
  nb_simulations = NULL,
  tip_data = NULL,
  trait_data_type,
  focal_time,
  update_Map = FALSE,
  keep_tip_labels = TRUE
)

Arguments

contMaps	For continuous trait data. List of objects of class "contMap", typically generated with prepare_trait_data(), each containing a phylogenetic tree and associated continuous trait mapping that represents an independent evolutionary history of ancestral trait evolution, conditioned on observed trait data and model fit (i.e., stochastic maps). The phylogenetic tree must be rooted and fully resolved/dichotomous, but it does not need to be ultrametric (it can include fossils).
densityMaps	For categorical trait or biogeographic data. List of objects of class "densityMap", typically generated with prepare_trait_data(), that contains a phylogenetic tree and associated posterior probability of being in a given state/range along branches. Each object (i.e., densityMap) corresponds to a state/range. The phylogenetic tree must be rooted and fully resolved/dichotomous, but it does not need to be ultrametric (it can include fossils).
simmaps	For categorical trait or biogeographic data. List of objects of classes "phylo" and "simmap", typically generated with prepare_trait_data() or phytools::make.simmap(), that represent discrete character/geographic evolutionary history (i.e., transitions in character states/geographic ranges) mapped along branches. This is needed to be able to track which simulated history provided which trait data in downstream analyses employing the "paired" or "full" strategies to account for uncertainty in ancestral trait estimates.
nb_simulations	For categorical trait or biogeographic data. Integer. The number of stochastic maps used to simulate trait evolution. This is needed for the "paired" and "full" strategies to account for trait estimate uncertainty, if only densityMaps summarizing posterior state/range density are provided, but not the simmaps representing all evolutionary histories.
tip_data	(Optional) Named vector of tip values of the trait. For continuous trait data: Named numeric vector of trait values. For categorical trait or biogeographic data: Character string vector of states/ranges Names are nodes_ID of the internal nodes. Ensure accurate tip values are used.
trait_data_type	Character string. Specify the type of trait data. Must be one of "continuous", "categorical", "biogeographic".
focal_time	Integer. The time, in terms of time distance from the present, at which the tree and mapping must be cut. It must be smaller than the root age of the phylogeny.
update_Map	Logical. Specify whether the mapped phylogeny (contMap, densityMaps, and/or simmaps) provided as input should be updated for visualization and returned among the outputs. Default is FALSE. The update consists of cutting off branches and mapping that are younger than the focal_time.
keep_tip_labels	Logical. Specify whether terminal branches with a single descendant tip must retain their initial tip.label on the updated contMap. Default is TRUE. Used only if update_Map = TRUE.

Details

The mapped phylogenies (contMaps, densityMaps, or simmaps) are cut at a specific time in the past (i.e. the focal_time) and the associated trait values of the overlapping edges/branches are extracted.

—– Extract trait_data —–

For continuous trait data:

Simulated trait data are extracted from contMaps (i.e. continuous stochastic maps).

For categorical trait and biogeographic data:

If simmaps are provided as input, all states/ranges are extracted directly from the stochastic maps.

If only densityMaps are provided, the posterior probabilities of states/ranges are extracted. Posterior probabilities are multiplied by the nb_simulations to record the distribution of states/ranges across dummy stochastic maps and assigned to each tip and cut branches at focal_time.

True tip states/ranges will be used if tip_data are provided as optional inputs. In practice the discrepancy is negligible.

—– Update the contMap/densityMaps/simmaps —–

To obtain an updated contMap/densityMaps/simmaps alongside the trait data, set update_Map = TRUE. The update consists of cutting off branches and mapping that are younger than the focal_time.

• When a branch with a single descendant tip is cut and keep_tip_labels = TRUE, the leaf left is labeled with the tip.label of the unique descendant tip.

• When a branch with a single descendant tip is cut and keep_tip_labels = FALSE, the leaf left is labeled with the node ID of the unique descendant tip.

• In all cases, when a branch with multiple descendant tips (i.e., a clade) is cut, the leaf left is labeled with the node ID of the MRCA of the cut-off clade.

The mapping in contMap/densityMaps/simmaps (⁠$maps⁠ and ⁠$mapped.edge⁠) is updated accordingly by removing mapping associated with the cut off branches.

To extract only the most likely trait value/state/range, see this associated function: extract_most_likely_trait_values_for_focal_time()

Value

By default, the function returns a list with three elements.

• ⁠$trait_data⁠ A list of named numeric vectors with simulated trait values found along branches overlapping the focal_time across each of the stochastic maps. Names are the tip.label/tipward node ID. ID of the associated stochastic maps associated with each item of the list can only be provided if simmaps are provided as inputs. If only densityMaps are provided, the distribution of states/ranges is associated with dummy maps_ID.

• ⁠$focal_time⁠ Integer. The time, in terms of time distance from the present, at which the trait data were extracted.

• ⁠$trait_data_type⁠ Character string. Define the type of trait data as "continuous", "categorical", or "biogeographic". Used in downstream analyses to select appropriate statistical processing.

If update_Map = TRUE, the output is a list with four to five elements: ⁠$trait_data⁠, ⁠$focal_time⁠, ⁠$trait_data_type⁠, ⁠$contMaps⁠ or ⁠$densityMaps⁠, and/or ⁠$simmaps⁠.

For continuous trait data:

• ⁠$contMaps⁠ A list of objects of class "contMap" that contains the updated contMaps with branches and mapping that are younger than the focal_time cut off. The function also adds multiple useful sub-elements to the ⁠$contMap$tree⁠ elements.

  • ⁠$root_age⁠ Integer. Stores the age of the root of the tree.

  • ⁠$nodes_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_node_ID to the initial_node_ID. Each row is a node.

  • ⁠$initial_nodes_ID⁠ Character vector. Provides the initial ID of internal nodes. Used to plot internal node IDs as labels with ape::nodelabels().

  • ⁠$edges_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_edge_ID to the initial_edge_ID. Each row is an edge/branch.

  • ⁠$initial_edges_ID⁠ Character vector. Provides the initial ID of edges/branches. Used to plot edge/branch IDs as labels with ape::edgelabels().

For categorical trait and biogeographic data:

• ⁠$densityMaps⁠ A list of objects of class "densityMap" that contains the updated densityMap of each state/range, with branches and mapping that are younger than the focal_time cut off. The function also adds multiple useful sub-elements to the ⁠$densityMaps$tree⁠ elements.

  • ⁠$root_age⁠ Integer. Stores the age of the root of the tree.

  • ⁠$nodes_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_node_ID to the initial_node_ID. Each row is a node.

  • ⁠$initial_nodes_ID⁠ Character vector. Provides the initial ID of internal nodes. Used to plot internal node IDs as labels with ape::nodelabels().

  • ⁠$edges_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_edge_ID to the initial_edge_ID. Each row is an edge/branch.

  • ⁠$initial_edges_ID⁠ Character vector. Provides the initial ID of edges/branches. Used to plot edge/branch IDs as labels with ape::edgelabels().

• ⁠$simmaps⁠ A list of objects with the class "simmap" that contains the updated simmap = stochastic maps representing discrete character/geographic evolutionary history, with branches and mapping that are younger than the focal_time cut off. The function also adds multiple useful sub-elements to each simmap.

  • ⁠$root_age⁠ Integer. Stores the age of the root of the tree.

  • ⁠$nodes_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_node_ID to the initial_node_ID. Each row is a node.

  • ⁠$initial_nodes_ID⁠ Character vector. Provides the initial ID of internal nodes. Used to plot internal node IDs as labels with ape::nodelabels().

  • ⁠$edges_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_edge_ID to the initial_edge_ID. Each row is an edge/branch.

  • ⁠$initial_edges_ID⁠ Character vector. Provides the initial ID of edges/branches. Used to plot edge/branch IDs as labels with ape::edgelabels().

Author(s)

Maël Doré

See Also

cut_phylo_for_focal_time() cut_contMaps_for_focal_time() cut_densityMaps_for_focal_time() cut_simmaps_for_focal_time()

Equivalent function to extract the most likely trait value/state/range from mapped phylogenies: extract_most_likely_trait_values_for_focal_time()

Examples

# ----- Example 1: Continuous trait ----- #

if (deepSTRAPP::is_dev_version())
{
 ## The R package 'contsimmap' is needed for this example to work.
 # Please install it manually from: https://github.com/bstaggmartin/contsimmap.

 ## Prepare data

 # Load eel data from the R package phytools
 # Source: Collar et al., 2014; DOI: 10.1038/ncomms6505

 library(phytools)
 data(eel.tree)
 data(eel.data)

 # Extract body size
 eel_data <- setNames(eel.data$Max_TL_cm,
                      rownames(eel.data))

  # (May take several minutes to run)
 ## Map continuous trait evolution on the phylogeny
 eel_cont_data  <- prepare_trait_data(
    tip_data = eel_data,
    trait_data_type = "continuous",
    phylo = eel.tree,
    # evolutionary_models = c("BM", "OU", "lambda", "kappa"),
    evolutionary_models = c("BM"),
    # Perform stochastic mapping to obtain multiple evolutionary histories
    run_stochastic_maps = TRUE,
    nb_simulations = 100,
    verbose = TRUE)

 # Extract continuous stochastic maps (contMaps)
 eel_contMaps <- eel_cont_data$contMaps

 # Plot the interpolated map of ML estimates
 plot_contMap(contMap = eel_cont_data$contMap)

 # Plot several continuous stochastic maps (contMaps)
 plot_contMap(contMap = eel_contMaps[[1]])
 plot_contMap(contMap = eel_contMaps[[10]])
 plot_contMap(contMap = eel_contMaps[[100]])

 # Set focal time to 10 Mya
 focal_time <- 10

 ## Extract all trait values for focal time = 10 Mya
 extract_trait_data_10My <- extract_all_trait_values_for_focal_time(
    contMaps = eel_contMaps, trait_data_type = "continuous",
    focal_time = 10, update_Map = TRUE)

 # Convert in data.frame
  # Rows = Stochastic maps
  # Columns = Cut branches at 10 Mya
 trait_data_df <- as.data.frame(do.call(rbind, extract_trait_data_10My$trait_data))
 head(trait_data_df)

 ## Plot updated contMaps

 updated_contMaps_10My <- extract_trait_data_10My$contMaps
 plot_contMap(contMap = updated_contMaps_10My[[1]])
 plot_contMap(contMap = updated_contMaps_10My[[10]])
 plot_contMap(contMap = updated_contMaps_10My[[100]])
 
}


# ----- Example 2: Categorical trait ----- #

## Load categorical trait data mapped on a phylogeny
data(eel_cat_3lvl_data, package = "deepSTRAPP")

# Explore data
str(eel_cat_3lvl_data, 1)
eel_cat_3lvl_data$densityMaps # Three density maps: one per state

# Set focal time to 10 Mya
focal_time <- 10

 # (May take several minutes to run)
## Extract trait data and update densityMaps for the given focal_time

# Extract from the densityMaps
eel_cat_3lvl_data_10My <- extract_all_trait_values_for_focal_time(
   densityMaps = eel_cat_3lvl_data$densityMaps,
   trait_data_type = "categorical",
   nb_simulations = 100,
   focal_time = focal_time,
   update_Map = TRUE)

## Print trait data
str(eel_cat_3lvl_data_10My, 1)
eel_cat_3lvl_data_10My$trait_data[1:2]

# Convert in data.frame
 # Rows = Dummy stochastic maps
 # Columns = Cut branches at 10 Mya
trait_data_df <- as.data.frame(do.call(rbind, eel_cat_3lvl_data_10My$trait_data))
trait_data_df[50:60, ]

# Distributions of states across dummy stochastic maps reflect posterior probabilities
# recorded in densityMaps, but they are not true evolutionary histories.
# If you want to relate trait states with a simulated evolutionary history,
# you need to provide simmaps as input.



# ----- Example 3: Biogeographic ranges ----- #

## Load biogeographic range data mapped on a phylogeny
data(eel_biogeo_data, package = "deepSTRAPP")

# Explore data
str(eel_biogeo_data, 1)
length(eel_biogeo_data$simmaps) # 100 stochastic maps: one per simulated biogeographic history

# Set focal time to 10 Mya
focal_time <- 10

 # (May take several minutes to run)
## Extract trait data and update simmaps for the given focal_time

# Extract from the simmaps
eel_biogeo_data_10My <- extract_all_trait_values_for_focal_time(
   simmaps = eel_biogeo_data$simmaps,
   trait_data_type = "biogeographic",
   focal_time = focal_time,
   update_Map = TRUE)

## Print trait data
str(eel_biogeo_data_10My, 1)
eel_biogeo_data_10My$trait_data[1:2]

# Convert in data.frame
 # Rows = Stochastic maps
 # Columns = Cut branches at 10 Mya
trait_data_df <- as.data.frame(do.call(rbind, eel_biogeo_data_10My$trait_data))
trait_data_df[1:10, ]

# Distributions of ranges are recorded across true stochastic maps
# as you provided simmaps as input.

## Plot updated stochastic maps

# Plot initial stochastic map n°1
plot(eel_biogeo_data$simmaps[[1]], fsize = 0.5)
abline(v = max(phytools::nodeHeights(eel_biogeo_data$simmaps[[1]])[,2]) - focal_time,
       col = "red", lty = 2, lwd = 2)

# Plot updated stochastic map n°1, cut at 10 Mya
plot(eel_biogeo_data_10My$simmaps[[1]], fsize = 0.7)

---

Extract diversification data from a BAMM_object

Description

Extracts regime IDs and tip rates from a BAMM_object that has been updated to provide diversification data for a specific time in the past (i.e. the focal_time). Use update_rates_and_regimes_for_focal_time() to obtain a BAMM_object updated for a given focal_time.

Usage

extract_diversification_data_melted_df_for_focal_time(
  BAMM_object,
  verbose = TRUE
)

Arguments

BAMM_object	Object of class "bammdata", typically generated with update_rates_and_regimes_for_focal_time(), that contains a phylogenetic tree and associated diversification rate mapping across selected posterior samples.
verbose	Logical. Should progression be displayed? A message will be printed for every batch of 100 BAMM posterior samples extracted. Default is TRUE.

Value

Returns a data.frame with six columns.

• ⁠$focal_time⁠ Integer. The time, in terms of time distance from the present, at which the trait data were extracted. Should be equal for all rows as a unique BAMM_object updated for a unique focal_time is being extracted.

• ⁠$BAMM_sample_ID⁠ Character string. ID of the posterior samples from which the diversification data are extracted named as 'BAMM_X' where X is the numerical ID.

• ⁠$tip_ID⁠ Character string. Tip labels of the branches cut-off at focal_time.

  • If keep_tip_labels = TRUE was used in update_rates_and_regimes_for_focal_time(), cut-off branches with a single descendant tip retain their initial tip.label.

  • If keep_tip_labels = FALSE was used in update_rates_and_regimes_for_focal_time(), all cut-off branches are labeled using their tipward node ID.

• ⁠$regime_ID⁠ Integer. The regime ID on tips at focal_time. IDs are integer. The root process equals "1", then they are incremented from oldest to youngest. Regime IDs are independent across posterior samples.

• ⁠$rate_type⁠ Character string. Type of rates: "lambda" for speciation rates, "mu" for extinction rates, and "net_diversification" for net diversification rates (lambda - mu).

• ⁠$rates⁠ Numeric. Rates in \[number of events / branch / evolutionary time\].

Author(s)

Maël Doré

Examples

if (deepSTRAPP::is_dev_version())
{
 # Load the BAMM_object summarizing 1000 posterior samples of BAMM updated for a focal_time of 10 My
 data(Ponerinae_BAMM_object_10My)
 ## This dataset is only available in development versions installed from GitHub.
 # It is not available in CRAN versions.
 # Use remotes::install_github(repo = "MaelDore/deepSTRAPP") to get the latest development version.

  # (May take several minutes to run)
 # Extract diversification data
 diversification_data_df <- extract_diversification_data_melted_df_for_focal_time(
     BAMM_object = Ponerinae_BAMM_object_10My,
     verbose = TRUE)
 # Print output
 head(diversification_data_df) 
}

---

Extract most likely trait data mapped on a phylogeny at a given time in the past

Description

Extracts the most likely trait values/states/ranges found along branches at a specific time in the past (i.e. the focal_time). Optionally, the function can update the mapped phylogeny (contMap or densityMaps) such as branches overlapping the focal_time are shortened to the focal_time, and the trait mapping for the cut off branches are removed by updating the ⁠$tree$maps⁠ and ⁠$tree$mapped.edge⁠ elements.

Usage

extract_most_likely_trait_values_for_focal_time(
  contMap = NULL,
  contMaps = NULL,
  densityMaps = NULL,
  simmaps = NULL,
  ace = NULL,
  tip_data = NULL,
  trait_data_type,
  focal_time,
  update_Map = FALSE,
  keep_tip_labels = TRUE
)

Arguments

contMap	For continuous trait data. Object of class "contMap", typically generated with prepare_trait_data() or phytools::contMap(), that contains a phylogenetic tree and associated continuous trait mapping, representing interpolated ML estimates of ancestral trait values. The phylogenetic tree must be rooted and fully resolved/dichotomous, but it does not need to be ultrametric (it can include fossils).
contMaps	For continuous trait data. List of objects of class "contMap", typically generated with prepare_trait_data(), each containing a phylogenetic tree and associated continuous trait mapping that represents an independent evolutionary history of ancestral trait evolution, conditioned on observed trait data and model fit (i.e., stochastic maps).
densityMaps	For categorical trait or biogeographic data. List of objects of class "densityMap", typically generated with prepare_trait_data(), that contains a phylogenetic tree and associated posterior probability of being in a given state/range along branches. Each object (i.e., densityMap) corresponds to a state/range. The phylogenetic tree must be rooted and fully resolved/dichotomous, but it does not need to be ultrametric (it can include fossils).
simmaps	For categorical trait or biogeographic data. List of objects of classes "phylo" and "simmap", typically generated with prepare_trait_data() or phytools::make.simmap(), that represent discrete character/geographic evolutionary history (i.e., transitions in character states/geographic ranges) mapped along branches.
ace	(Optional) Ancestral Character Estimates (ACE) at the internal nodes. Obtained with prepare_trait_data() as output in the ⁠$ace⁠ slot. For continuous trait data: Named numeric vector typically generated with phytools::fastAnc(), phytools::anc.ML(), or ape::ace(). Names are nodes_ID of the internal nodes. Values are ACE of the trait. For categorical trait or biogeographic data: Matrix that records the posterior probabilities of ancestral states/ranges. Rows are internal nodes_ID. Columns are states/ranges. Values are posterior probabilities of each state per node. Needed in all cases to provide accurate estimates of trait values.
tip_data	(Optional) Named vector of tip values of the trait. For continuous trait data: Named numeric vector of trait values. For categorical trait or biogeographic data: Character string vector of states/ranges Names are nodes_ID of the internal nodes. Needed to provide accurate tip values.
trait_data_type	Character string. Specify the type of trait data. Must be one of "continuous", "categorical", "biogeographic".
focal_time	Integer. The time, in terms of time distance from the present, at which the tree and mapping must be cut. It must be smaller than the root age of the phylogeny.
update_Map	Logical. Specify whether the mapped phylogeny (contMap, densityMaps, simmaps) provided as input should be updated for visualization and returned among the outputs. Default is FALSE. The update consists of cutting off branches and mapping that are younger than the focal_time.
keep_tip_labels	Logical. Specify whether terminal branches with a single descendant tip must retain their initial tip.label on the updated contMap. Default is TRUE. Used only if update_Map = TRUE.

Details

The mapped phylogeny (contMap or densityMaps) is cut at a specific time in the past (i.e. the focal_time) and the current trait values of the overlapping edges/branches are extracted.

—– Extract trait_data —–

For continuous trait data:

contMap is the default input representing the interpolated ML estimates of trait values along the phylogeny.

If providing only the contMap trait values at tips and internal nodes will be extracted from the mapping of the contMap leading to a slight discrepancy with the actual tip data and estimated ancestral character values.

True ML trait estimates will be used if tip_data and/or ace are provided as optional inputs. In practice the discrepancy is negligible.

Alternatively to a unique contMap, the user can provide a full set of continuous stochastic maps as contMaps. Each map then represents an independent evolutionary history of ancestral trait evolution, conditioned on observed trait data and model fit (i.e., stochastic maps). The 'most likely' trait values will be extracted as the mean of values observed across the stochastic maps. This quantity closely approximates the ML estimates as typically provided with a contMap.

For categorical trait and biogeographic data:

densityMaps are the default input. Most likely states/ranges are directly extracted from the posterior probabilities displayed in the densityMaps. The state/range with the highest probability is assigned to each tip and cut branches at focal_time.

True ML states/ranges will be used if tip_data and/or ace are provided as optional inputs. In practice the discrepancy is negligible.

Alternatively to densityMaps, the user can provide a full set of stochastic maps as simmaps. Each map then represents an independent evolutionary history of ancestral trait evolution, conditioned on observed trait data and model fit (i.e., stochastic maps). The 'most likely' states / ranges will be extracted as the most frequent state/range observed across the stochastic maps. This quantity is fully equivalent to what is derived from the densityMaps.

—– Update the contMap(s)/densityMaps/simmaps —–

To obtain updated contMap(s)/densityMaps/simmaps alongside the trait data, set update_Map = TRUE. The update consists of cutting off branches and mapping that are younger than the focal_time.

• When a branch with a single descendant tip is cut and keep_tip_labels = TRUE, the leaf left is labeled with the tip.label of the unique descendant tip.

• When a branch with a single descendant tip is cut and keep_tip_labels = FALSE, the leaf left is labeled with the node ID of the unique descendant tip.

• In all cases, when a branch with multiple descendant tips (i.e., a clade) is cut, the leaf left is labeled with the node ID of the MRCA of the cut-off clade.

The mapping in contMap(s)/densityMaps (⁠$tree$maps⁠ and ⁠$tree$mapped.edge⁠) or simmaps (⁠$maps⁠ and ⁠$mapped.edge⁠) is updated accordingly by removing mapping associated with the cut off branches.

To extract all trait data across multiple stochastic maps, and not just the most likely trait value/state/range, see this associated function: extract_all_trait_values_for_focal_time()

Value

By default, the function returns a list with three elements.

• ⁠$trait_data⁠ A named numeric vector with ML trait values found along branches overlapping the focal_time. Names are the tip.label/tipward node ID.

• ⁠$focal_time⁠ Integer. The time, in terms of time distance from the present, at which the trait data were extracted.

• ⁠$trait_data_type⁠ Character string. Define the type of trait data as "continuous", "categorical", or "biogeographic". Used in downstream analyses to select appropriate statistical processing.

If update_Map = TRUE, the output also contains updated mapped phylogenies as: ⁠$contMap⁠ and/or ⁠$contMaps⁠ or ⁠$densityMaps⁠ and/or ⁠$simmaps⁠.

For continuous trait data:

• ⁠$contMap⁠ An object of class "contMap" that contains the updated contMap with branches and mapping that are younger than the focal_time cut off. The function also adds multiple useful sub-elements to the ⁠$contMap$tree⁠ element.

  • ⁠$root_age⁠ Integer. Stores the age of the root of the tree.

  • ⁠$nodes_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_node_ID to the initial_node_ID. Each row is a node.

  • ⁠$initial_nodes_ID⁠ Character vector. Provides the initial ID of internal nodes. Used to plot internal node IDs as labels with ape::nodelabels().

  • ⁠$edges_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_edge_ID to the initial_edge_ID. Each row is an edge/branch.

  • ⁠$initial_edges_ID⁠ Character vector. Provides the initial ID of edges/branches. Used to plot edge/branch IDs as labels with ape::edgelabels().

• ⁠$contMaps⁠ List of "contMap" objects as described above.

For categorical trait and biogeographic data:

• ⁠$densityMaps⁠ A list of objects of class "densityMap" that contains the updated densityMap of each state/range, with branches and mapping that are younger than the focal_time cut off. The function also adds multiple useful sub-elements to the ⁠$densityMaps$tree⁠ elements.

  • ⁠$root_age⁠ Integer. Stores the age of the root of the tree.

  • ⁠$nodes_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_node_ID to the initial_node_ID. Each row is a node.

  • ⁠$initial_nodes_ID⁠ Character vector. Provides the initial ID of internal nodes. Used to plot internal node IDs as labels with ape::nodelabels().

  • ⁠$edges_ID_df⁠ Data.frame with two columns. Provides the conversion from the new_edge_ID to the initial_edge_ID. Each row is an edge/branch.

  • ⁠$initial_edges_ID⁠ Character vector. Provides the initial ID of edges/branches. Used to plot edge/branch IDs as labels with ape::edgelabels().

• ⁠$simmaps⁠ A list of objects of classes "phylo" and "simmap", that contains the updated stochastic maps simmap for each simulation of trait/range evolution. The same useful sub-elements as described above are also added to each simmap object.

Author(s)

Maël Doré

See Also

cut_phylo_for_focal_time() cut_contMap_for_focal_time() cut_densityMaps_for_focal_time()

Equivalent function to extract all data for multiple stochastic maps extract_all_trait_values_for_focal_time()

Examples

# ----- Example 1: Continuous trait ----- #

## Prepare data

# Load eel data from the R package phytools
# Source: Collar et al., 2014; DOI: 10.1038/ncomms6505

library(phytools)
data(eel.tree)
data(eel.data)

# Extract body size
eel_data <- setNames(eel.data$Max_TL_cm,
                     rownames(eel.data))

 # (May take several minutes to run)
## Get Ancestral Character Estimates based on a Brownian Motion model
# To obtain values at internal nodes
eel_ACE <- phytools::fastAnc(tree = eel.tree, x = eel_data)

## Interpolate ML estimates of trait values based on a Brownian Motion model
# to interpolate values along branches and obtain a "contMap" object
eel_contMap <- phytools::contMap(eel.tree, x = eel_data,
                                 res = 100, # Number of time steps
                                 plot = FALSE)

# Set focal time to 50 Mya
focal_time <- 50

## Extract trait data and update contMap for the given focal_time

# Extract from the contMap (values are not exact ML estimates)
eel_cont_50 <- extract_most_likely_trait_values_for_focal_time(
   contMap = eel_contMap,
   trait_data_type = "continuous",
   focal_time = focal_time,
   update_Map = TRUE)
# Extract from tip data and ML estimates of ancestral characters (values are true ML estimates)
eel_cont_50 <- extract_most_likely_trait_values_for_focal_time(
   contMap = eel_contMap,
   ace = eel_ACE, tip_data = eel_data,
   trait_data_type = "continuous",
   focal_time = focal_time,
   update_Map = TRUE)

## Visualize outputs

# Print trait data
eel_cont_50$trait_data

# Plot node labels on initial stochastic map with cut-off
plot(eel_contMap, fsize = c(0.5, 1))
ape::nodelabels()
abline(v = max(phytools::nodeHeights(eel_contMap$tree)[,2]) - focal_time,
       col = "red", lty = 2, lwd = 2)

# Plot updated contMap with initial node labels
plot(eel_cont_50$contMap)
ape::nodelabels(text = eel_cont_50$contMap$tree$initial_nodes_ID) 


# ----- Example 2: Categorical trait ----- #

 # (May take several minutes to run)
## Load categorical trait data mapped on a phylogeny
data(eel_cat_3lvl_data, package = "deepSTRAPP")

# Explore data
str(eel_cat_3lvl_data, 1)
eel_cat_3lvl_data$densityMaps # Three density maps: one per state

# Set focal time to 10 Mya
focal_time <- 10

## Extract trait data and update densityMaps for the given focal_time

# Extract from the densityMaps
eel_cat_3lvl_data_10My <- extract_most_likely_trait_values_for_focal_time(
   densityMaps = eel_cat_3lvl_data$densityMaps,
   trait_data_type = "categorical",
   focal_time = focal_time,
   update_Map = TRUE)

## Print trait data
str(eel_cat_3lvl_data_10My, 1)
eel_cat_3lvl_data_10My$trait_data

## Plot density maps as overlay of all state posterior probabilities

# Plot initial density maps with ACE pies
plot_densityMaps_overlay(densityMaps = eel_cat_3lvl_data$densityMaps)
abline(v = max(phytools::nodeHeights(eel_cat_3lvl_data$densityMaps[[1]]$tree)[,2]) - focal_time,
       col = "red", lty = 2, lwd = 2)

# Plot updated densityMaps with ACE pies
plot_densityMaps_overlay(eel_cat_3lvl_data_10My$densityMaps) 


# ----- Example 3: Biogeographic ranges ----- #

 # (May take several minutes to run)
## Load biogeographic range data mapped on a phylogeny
data(eel_biogeo_data, package = "deepSTRAPP")

# Explore data
str(eel_biogeo_data, 1)
eel_biogeo_data$densityMaps # Two density maps: one per unique area: A, B.
eel_biogeo_data$densityMaps_all_ranges # Three density maps: one per range: A, B, and AB.

# Set focal time to 10 Mya
focal_time <- 10

## Extract trait data and update densityMaps for the given focal_time

# Extract from the densityMaps
eel_biogeo_data_10My <- extract_most_likely_trait_values_for_focal_time(
   densityMaps = eel_biogeo_data$densityMaps,
   # ace = eel_biogeo_data$ace,
   trait_data_type = "biogeographic",
   focal_time = focal_time,
   update_Map = TRUE)

## Print trait data
str(eel_biogeo_data_10My, 1)
eel_biogeo_data_10My$trait_data

## Plot density maps as overlay of all range posterior probabilities

# Plot initial density maps with ACE pies
plot_densityMaps_overlay(densityMaps = eel_biogeo_data$densityMaps)
abline(v = max(phytools::nodeHeights(eel_biogeo_data$densityMaps[[1]]$tree)[,2]) - focal_time,
       col = "red", lty = 2, lwd = 2)

# Plot updated densityMaps with ACE pies
plot_densityMaps_overlay(eel_biogeo_data_10My$densityMaps)

---

Extract trait data from a trait_data object

Description

Extract trait data from a trait_data_list object as produced within the deepSTRAPP workflow for a specific time in the past (i.e. the focal_time) and produce a melted dataset recording trait values across (stochastic maps x) edges for the associated focal_time.

Usage

extract_trait_data_melted_df_for_focal_time(trait_data_list)

Arguments

trait_data_list

Object summarizing trait data extracted for a given focal_time, including a ⁠$trait_data⁠ element storing trait data recorded across (stochastic maps x) edges for the associated ⁠$focal_time⁠. Use extract_all_trait_values_for_focal_time() to obtain trait data across multiple stochastic maps. Use extract_most_likely_trait_values_for_focal_time() to obtain only the most likely trait data values.

Value

Returns a data.frame with five columns.

• ⁠$focal_time⁠ Integer. The time, in terms of time distance from the present, at which the trait data were extracted.

• ⁠$Map_ID⁠ Character string. ID of the stochastic map from which the trait data are extracted. If using 'rate_only' strategy to account for uncertainty in ancestral estimate, this is fixed to "Map_ML". If using 'paired' or 'full' strategies, this records either true stochastic maps "Map_X", or dummy maps "Dummy_map_X", depending on whether you provided respectively stochastic maps (contMaps/simmaps) or simply densityMaps as inputs.

• ⁠$tip_ID⁠ Character string. Tip labels of the branches cut-off at focal_time.

  • If keep_tip_labels = TRUE was used in update_rates_and_regimes_for_focal_time(), cut-off branches with a single descendant tip retain their initial tip.label.

  • If keep_tip_labels = FALSE was used in update_rates_and_regimes_for_focal_time(), all cut-off branches are labeled using their tipward node ID.

• ⁠$trait_data_type⁠ Character string. Records the type of traits: "continuous", "categorical", or "biogeographic".

• ⁠$trait_value⁠ Numerical or Character string.

  • For "continuous" traits: numerical values.

  • For "categorical" traits: states.

  • For "biogeographic" data: ranges.

Author(s)

Maël Doré

Examples

# ----- Example 1: Extract ML estimates data ----- #

## Load categorical trait data mapped on a phylogeny
data(eel_cat_3lvl_data, package = "deepSTRAPP")

# Explore data
str(eel_cat_3lvl_data, 1)
eel_cat_3lvl_data$densityMaps # Three density maps: one per state

# Set focal time to 10 Mya
focal_time <- 10

## Extract ML estimates of ancestral states for the given focal_time

# Extract from the densityMaps
eel_cat_3lvl_data_10My <- extract_most_likely_trait_values_for_focal_time(
   densityMaps = eel_cat_3lvl_data$densityMaps,
   trait_data_type = "categorical",
   focal_time = focal_time)

## Format ML states data as a melted df

melted_df <- extract_trait_data_melted_df_for_focal_time(trait_data_list = eel_cat_3lvl_data_10My)
head(melted_df)


# ----- Example 2: Extract data across multiple stochastic maps ----- #

## Load biogeographic range data mapped on a phylogeny
data(eel_biogeo_data, package = "deepSTRAPP")

# Explore data
str(eel_biogeo_data, 1)
length(eel_biogeo_data$simmaps) # 100 stochastic maps: one per simulated biogeographic history

# Set focal time to 10 Mya
focal_time <- 10

## Extract range data across stochastic maps for the given focal_time

# Extract from the simmaps
eel_biogeo_data_10My <- extract_all_trait_values_for_focal_time(
   simmaps = eel_biogeo_data$simmaps,
   trait_data_type = "biogeographic",
   focal_time = focal_time)

## Format range data as a melted df

melted_df <- extract_trait_data_melted_df_for_focal_time(trait_data_list = eel_biogeo_data_10My)
head(melted_df)

---

Phylogeny and body mass data for extant and extinct mammal families/genera from Slater, 2013

Description

A list containing two elements:

• ⁠$mammal.mass⁠ A data.frame with mean and standard error of mammal body masses in ln(mass in g).

• ⁠$mammal.phy⁠ A time-calibrated phylogeny of extinct and extant mammal families/genera.

Source: Slater, G. J. (2013). Phylogenetic evidence for a shift in the mode of mammalian body size evolution at the Cretaceous-Palaeogene boundary. Methods in Ecology and Evolution, 4(8), 734-744. doi:10.1111/2041-210X.12084

Usage

data(mammals)

Format

A list with 2 elements.

• ⁠$mammal.mass⁠ A data.frame of 213 observations and two columns. Values are numerical.

• ⁠$mammal.phy⁠ A time-calibrated phylogeny of 211 tips.

Details

Initial dataset from Slater, G. J. (2013). The R object was imported from R package motmot to include in deepSTRAPP.

---

Plot diversification rates and regime shifts from BAMM on phylogeny

Description

Plot on a time-calibrated phylogeny the evolution of diversification rates and the location of regime shifts estimated from a BAMM (Bayesian Analysis of Macroevolutionary Mixtures). Each branch is colored according to the estimated rates of speciation, extinction, or net diversification stored in an object of class "bammdata". Rates can vary along time, thus colors evolve along individual branches.

This function is a wrapper of original functions from the R package {BAMMtools}:

• Step 1: Use BAMMtools::plot.bammdata() to map rates on the phylogeny.

• Step 2: Add the location of regime shifts with BAMMtools::addBAMMshifts() (if add_regime_shifts = TRUE).

Usage

plot_BAMM_rates(
  BAMM_object,
  rate_type = "net_diversification",
  method = "phylogram",
  configuration_type = "MAP",
  sample_index = 1,
  regimes_fill = "grey",
  regimes_size = 1,
  regimes_pch = 21,
  regimes_border_col = "black",
  regimes_border_width = 1,
  ...,
  add_regime_shifts = TRUE,
  adjust_size_to_prob = TRUE,
  display_plot = TRUE,
  PDF_file_path = NULL
)

Arguments

BAMM_object	Object of class "bammdata", typically generated with prepare_diversification_data(), that contains a phylogenetic tree and associated diversification rate mapping across selected posterior samples. It works also for BAMM_object updated for a specific focal_time using update_rates_and_regimes_for_focal_time(), or the deepSTRAPP workflow with run_deepSTRAPP_for_focal_time and run_deepSTRAPP_over_time().
rate_type	A character string specifying the type of diversification rates to plot. Must be one of 'speciation', 'extinction' or 'net_diversification' (default).
method	A character string indicating the method for plotting the phylogenetic tree. method = "phylogram" (default) plots the phylogenetic tree using rectangular coordinates. method = "polar" plots the phylogenetic tree using polar coordinates.
configuration_type	A character string specifying how to select the location of regime shifts across posterior samples. configuration_type = "MAP": Use the average locations recorded in posterior samples with the Maximum A Posteriori probability (MAP) configuration. This regime shift configuration is the most frequent configuration among the posterior samples (See BAMMtools::getBestShiftConfiguration()). This is the default option. configuration_type = "MSC": Use the average locations recorded in posterior samples with the Maximum Shift Credibility (MSC) configuration. This regime shift configuration has the highest product of marginal probabilities across branches (See BAMMtools::maximumShiftCredibility()). configuration_type = "index": Use the configuration of a unique posterior sample whose index is provided in sample_index.
sample_index	Integer. Index of the posterior samples to use to plot the location of regime shifts. Used only if configuration_type = index. Default = 1.
regimes_fill	Character string. Set the color of the background of the symbols showing the location of regime shifts. Equivalent to the bg argument in BAMMtools::addBAMMshifts(). Default is "grey".
regimes_size	Numeric. Set the size of the symbols showing the location of regime shifts. Equivalent to the cex argument in BAMMtools::addBAMMshifts(). Default is 1.
regimes_pch	Integer. Set the shape of the symbols showing the location of regime shifts. Equivalent to the pch argument in BAMMtools::addBAMMshifts(). Default is 21.
regimes_border_col	Character string. Set the color of the border of the symbols showing the location of regime shifts. Equivalent to the col argument in BAMMtools::addBAMMshifts(). Default is "black".
regimes_border_width	Numeric. Set the width of the border of the symbols showing the location of regime shifts. Equivalent to the lwd argument in BAMMtools::addBAMMshifts(). Default is 1.
...	Additional graphical arguments to pass down to BAMMtools::plot.bammdata(), BAMMtools::addBAMMshifts(), and par(). Among them, par.reset is ignored, with a message: it would make BAMMtools::plot.bammdata() restore a full par() snapshot, discarding the coordinate system of the phylogeny and rewinding the panel counter of a multi-panel layout. Graphical parameters are instead restored by plot_BAMM_rates() when it exits: every parameter the call actually changed is put back, except those describing the plot itself, which are kept so that the phylogeny can be annotated afterwards with for instance graphics::abline() or graphics::title(), as in the examples.
add_regime_shifts	Logical. Whether to add the location of regime shifts on the phylogeny (Step 2). Default is TRUE. Provide the full argument name to avoid ambiguity with the add graphical argument.
adjust_size_to_prob	Logical. Whether to scale the size of the symbols showing the location of regime shifts according to the marginal shift probability of the shift happening on each location/branch. This will only work if there is an ⁠$MSP_tree⁠ element summarizing the marginal shift probabilities across branches in the BAMM_object. Default is TRUE. Provide the full argument name to avoid ambiguity with the adj graphical argument.
display_plot	Logical. Whether to display the plot generated in the R console. Default is TRUE.
PDF_file_path	Character string. If provided, the plot will be saved in a PDF file following the path provided here. The path must end with ".pdf".

Details

The main input BAMM_object is the typical output of prepare_diversification_data(). It provides information on rates and regime shifts across the posterior samples of a BAMM.

⁠$MAP_BAMM_object⁠ and ⁠$MSC_BAMM_object⁠ elements are required to plot regime shift locations following the "MAP" or "MSC" configuration_type respectively. A ⁠$MSP_tree⁠ element is required to scale the size of the symbols showing the location of regime shifts according to marginal shift probabilities. (If adjust_size_to_prob = TRUE).

The default option to display regime shift is to use the average locations from the posterior samples with the Maximum A Posteriori probability (MAP) configuration. However, sometimes, multiple configurations have similarly high frequency in the posterior samples (See BAMMtools::credibleShiftSet()). An alternative is to use the average locations from posterior samples with the Maximum Shift Credibility (MSC) configuration instead. This regime shift configuration has the highest product of marginal probabilities across branches where a shift is estimated. It may differ from the MAP configuration. (See BAMMtools::maximumShiftCredibility()).

Value

The function returns (invisibly) a list with three elements similarly to BAMMtools::plot.bammdata().

• ⁠$coords⁠: A matrix of plot coordinates. Rows correspond to branches. Columns 1-2 are starting (x,y) coordinates of each branch and columns 3-4 are ending (x,y) coordinates of each branch. If method = "polar" a fifth column gives the angle (in radians) of each branch.

• ⁠$colorbreaks⁠: A vector of percentiles used to group macroevolutionary rates into color bins.

• ⁠$colordens⁠: A matrix of the kernel density estimates (column 2) of evolutionary rates (column 1) and the color (column 3) corresponding to each rate value.

Author(s)

Maël Doré

Original functions by Mike Grundler & Pascal Title in R package {BAMMtools}.

See Also

Initial functions in BAMMtools: BAMMtools::plot.bammdata() BAMMtools::addBAMMshifts()

Associated functions in deepSTRAPP: prepare_diversification_data() update_rates_and_regimes_for_focal_time() run_deepSTRAPP_for_focal_time() run_deepSTRAPP_over_time()

Examples

# Load BAMM output
data(whale_BAMM_object, package = "deepSTRAPP")

mfrow_ini <- par()$mfrow
par(mfrow = c(1,2))

## Plot overall mean rates with MAP configuration for regime shifts
# (rates are averaged across all posterior samples)
plot_BAMM_rates(whale_BAMM_object, # Use the full object to plot mean overall rates
                add_regime_shifts = TRUE,
                configuration_type = "MAP",
                regimes_size = 3, bg = "black")
title("Overall mean - MAP shift config", col.main = "white")

## Plot overall mean rates with MSC configuration for regime shifts
# (rates are averaged across all posterior samples)
plot_BAMM_rates(whale_BAMM_object, add_regime_shifts = TRUE,
                configuration_type = "MSC",
                regimes_size = 3, bg = "black")
title("Overall mean - MSC shift config", col.main = "white")

par(mfrow_ini)


mfrow_ini <- par()$mfrow
par(mfrow = c(1,2))

## Plot mean MAP rates with regime shifts
# (rates are averaged only across MAP samples)
plot_BAMM_rates(whale_BAMM_object$MAP_BAMM_object, # Use the MAP object to plot mean MAP rates
                add_regime_shifts = TRUE,
                configuration_type = "index",
                # Set to index to use the regime shift location from MAP configuration
                regimes_size = 3)
title("MAP mean - MAP shift config")

## Plot mean MSC rates with regime shifts
# (rates averaged only across MSC samples)
plot_BAMM_rates(whale_BAMM_object$MSC_BAMM_object, # Use the MSC object to plot mean MSC rates
                add_regime_shifts = TRUE,
                configuration_type = "index",
                # Set to index to use the regime shift data from MSC configuration
                regimes_size = 3)
title("MSC mean - MSC shift config")

par(mfrow_ini)

---

Plot continuous trait evolution on the tree

Description

Plot on a time-calibrated phylogeny the evolution of a continuous trait as summarized in a contMap object typically generated with prepare_trait_data().

This function is a wrapper of original functions from the R package {phytools}:

• Step 1: Use phytools::setMap() to update the color scale if requested.

• Step 2: Use phytools::plot.contMap() to plot the mapped phylogeny.

Usage

plot_contMap(
  contMap,
  ...,
  color_scale = NULL,
  display_plot = TRUE,
  PDF_file_path = NULL
)

Arguments

contMap	List of class "contMap", typically generated with prepare_trait_data(), that contains a phylogenetic tree and associated ancestral trait estimates mapped along branches in contMap$tree$maps.
...	Additional arguments to pass down to phytools::plot.contMap() to control plotting.
color_scale	Character vector. List of colors to use to build the color scale with grDevices::colorRampPalette() showing the evolution of a continuous trait. From lowest values to highest values. Default (color_scale = NULL) is using the color palette recorded in the contMap$cols item. If none was provided, the rainbow() palette is used. Provide the full argument name to avoid ambiguity with the col graphical argument.
display_plot	Logical. Whether to display the plot generated in the R console. Default is TRUE.
PDF_file_path	Character string. If provided, the plot will be saved in a PDF file following the path provided here. The path must end with ".pdf".

Details

This function is a wrapper of original functions phytools::setMap() and phytools::plot.contMap(). Additions are listed below:

• The color scale can be controlled directly with the argument color_scale.

• The plot can be exported in PDF using PDF_file_path to define the output file.

Value

If display_plot = TRUE, the function plots the time-calibrated phylogeny displaying the evolution of a continuous trait. If PDF_file_path is provided, the function exports the plot into a PDF file.

An object of class "contMap" with an (optionally) updated color scale (⁠$cols⁠) is returned invisibly.

Author(s)

Maël Doré

Original functions by Liam Revell in R package {phytools}. Contact: [email protected]

See Also

phytools::plot.contMap() plot_densityMaps_overlay()

Examples

# Load phylogeny
data(Ponerinae_trait_tip_data, package = "deepSTRAPP")
# Load trait df
data(Ponerinae_tree, package = "deepSTRAPP")

## Prepare trait data

 # (May take several minutes to run)
# Extract continuous trait data as a named vector
Ponerinae_cont_tip_data <- setNames(object = Ponerinae_trait_tip_data$fake_cont_tip_data,
                                    nm = Ponerinae_trait_tip_data$Taxa)

# Get Ancestral Character Estimates based on a Brownian Motion model
# To obtain values at internal nodes
Ponerinae_ACE <- phytools::fastAnc(tree = Ponerinae_tree, x = Ponerinae_cont_tip_data)

# Infer Ancestral Character Estimates based on a Brownian Motion model
# and run a Stochastic Mapping to interpolate values along branches and obtain a "contMap" object
Ponerinae_contMap <- phytools::contMap(Ponerinae_tree, x = Ponerinae_cont_tip_data,
                                       res = 100, # Number of time steps
                                       plot = FALSE)

# Plot contMap with the phytools method
plot(x = Ponerinae_contMap, fsize = c(0.5, 1))

# Plot contMap with an updated color scale
plot_contMap(contMap = Ponerinae_contMap, fsize = c(0.5, 1),
             color_scale = c("darkgreen", "limegreen", "orange", "red"))
             # PDF_file_path = "Ponerinae_contMap.pdf")

---

Plot posterior probabilities of states/ranges on phylogeny from densityMaps

Description

Plot on a time-calibrated phylogeny the evolution of a categorical trait/biogeographic ranges summarized from densityMaps typically generated with prepare_trait_data(). Each branch is colored according to the posterior probability of being in a given state/range. Colors for each state/range are overlaid using transparency to produce a single plot for all states/ranges.

Usage

plot_densityMaps_overlay(
  densityMaps,
  ace = NULL,
  ...,
  colors_per_levels = NULL,
  add_ACE_pies = TRUE,
  cex_pies = 0.5,
  display_plot = TRUE,
  PDF_file_path = NULL
)

Arguments

densityMaps	List of objects of class "densityMap", typically generated with prepare_trait_data(), that contains a phylogenetic tree and associated posterior probability of being in a given state/range along branches. Each object (i.e., densityMap) corresponds to a state/range. If no color is provided for multi-area ranges, they will be interpolated.
ace	Numerical matrix. To provide the posterior probabilities of ancestral states/ranges (characters) estimates (ACE) at internal nodes used to plot the ACE pies. Rows are internal nodes. Columns are states/ranges. Values are posterior probabilities of each state per node. Typically generated with prepare_trait_data() in the ⁠$ace⁠ slot. If NULL (default), the ACE are extracted from the densityMaps with a possible slight discrepancy with the actual tip states and estimated posterior probabilities of ancestral states.
...	Additional arguments to pass down to phytools::plotSimmap() to control plotting.
colors_per_levels	Named character string. To set the colors to use to map each state/range posterior probabilities. Names = states/ranges; values = colors. If NULL (default), the color scale provided in the densityMaps will be used. Provide the full argument name to avoid ambiguity with the col graphical argument.
add_ACE_pies	Logical. Whether to add pies of posterior probabilities of states/ranges at internal nodes on the mapped phylogeny. Default = TRUE. Provide the full argument name to avoid ambiguity with the add graphical argument.
cex_pies	Numeric. To adjust the size of the ACE pies. Default = 0.5. Provide the full argument name to avoid ambiguity with the cex graphical argument.
display_plot	Logical. Whether to display the plot generated in the R console. Default is TRUE.
PDF_file_path	Character string. If provided, the plot will be saved in a PDF file following the path provided here. The path must end with ".pdf".

Value

If display_plot = TRUE, the function plots a time-calibrated phylogeny displaying the evolution of a categorical trait/biogeographic ranges. If PDF_file_path is provided, the function exports the plot into a PDF file.

Author(s)

Maël Doré

Original functions by Liam Revell in R package {phytools}. Contact: [email protected]

See Also

phytools::plot.densityMap() phytools::plotSimmap()

Examples

# Load phylogeny and tip data
library(phytools)
data(eel.tree)
data(eel.data)

# Transform feeding mode data into a 3-level factor
eel_data <- stats::setNames(eel.data$feed_mode, rownames(eel.data))
eel_data <- as.character(eel_data)
eel_data[c(1, 5, 6, 7, 10, 11, 15, 16, 17, 24, 25, 28, 30, 51, 52, 53, 55, 58, 60)] <- "kiss"
eel_data <- stats::setNames(eel_data, rownames(eel.data))
table(eel_data)

# Manually define a Q_matrix for rate classes of state transition to use in the 'matrix' model
# Does not allow transitions from state 1 ("bite") to state 2 ("kiss") or state 3 ("suction")
# Does not allow transitions from state 3 ("suction") to state 1 ("bite")
# Set symmetrical rates between state 2 ("kiss") and state 3 ("suction")
Q_matrix = rbind(c(NA, 0, 0), c(1, NA, 2), c(0, 2, NA))

# Set colors per state
colors_per_levels <- c("limegreen", "orange", "dodgerblue")
names(colors_per_levels) <- c("bite", "kiss", "suction")

 # (May take several minutes to run)
# Run evolutionary models to prepare trait data
eel_cat_3lvl_data <- prepare_trait_data(tip_data = eel_data, phylo = eel.tree,
    trait_data_type = "categorical",
    colors_per_levels = colors_per_levels,
    evolutionary_models = c("ER", "SYM", "ARD", "meristic", "matrix"),
    Q_matrix = Q_matrix,
    nb_simulations = 1000,
    plot_map = TRUE,
    plot_overlay = TRUE,
    return_best_model_fit = TRUE,
    return_model_selection_df = TRUE) 

# Load directly output
data(eel_cat_3lvl_data, package = "deepSTRAPP")

# Plot densityMaps one by one
plot(eel_cat_3lvl_data$densityMaps[[1]]) # densityMap for state n°1 ("bite")
plot(eel_cat_3lvl_data$densityMaps[[2]]) # densityMap for state n°1 ("kiss")
plot(eel_cat_3lvl_data$densityMaps[[3]]) # densityMap for state n°1 ("suction")

# Plot overlay of all densityMaps
plot_densityMaps_overlay(densityMaps = eel_cat_3lvl_data$densityMaps)

---

Plot histogram of STRAPP test statistics to assess results

Description

Plot a histogram of the distribution of the test statistics obtained from a STRAPP test carried out for a unique focal_time. (See compute_STRAPP_test_for_focal_time() and run_deepSTRAPP_for_focal_time()).

Returns a single histogram for overall tests. If plot_posthoc_tests = TRUE, it will return a faceted plot with a histogram per post hoc test.

If a PDF file path is provided in PDF_file_path, the plot will be saved directly in a PDF file.

Usage

plot_histogram_STRAPP_test_for_focal_time(
  deepSTRAPP_outputs,
  focal_time = NULL,
  display_plot = TRUE,
  plot_posthoc_tests = FALSE,
  PDF_file_path = NULL
)

Arguments

deepSTRAPP_outputs	List of elements generated with run_deepSTRAPP_for_focal_time(), that summarize the results of a STRAPP test for a specific time in the past (i.e. the focal_time). deepSTRAPP_outputs can also be extracted from the output of run_deepSTRAPP_over_time() that run the whole deepSTRAPP workflow and store results inside ⁠$STRAPP_results_over_time⁠.
focal_time	Numeric. (Optional) If deepSTRAPP_outputs comprises results over multiple time-steps (i.e., output of run_deepSTRAPP_over_time(), this is the time of the STRAPP test targeted for plotting.
display_plot	Logical. Whether to display the histogram(s) generated in the R console. Default is TRUE.
plot_posthoc_tests	Logical. For multinomial data only. Whether to plot the histogram for the overall Kruskal-Wallis test across all states (plot_posthoc_tests = FALSE), or plot the histograms for all the pairwise post hoc Dunn's tests across pairs of states (plot_posthoc_tests = TRUE). Default is FALSE.
PDF_file_path	Character string. If provided, the plot will be saved in a PDF file following the path provided here. The path must end with ".pdf".

Details

The main input deepSTRAPP_outputs is the typical output of run_deepSTRAPP_for_focal_time(). It provides information on results of a STRAPP test performed at a given focal_time.

Histograms are built based on the distribution of the test statistics. Such distributions are recorded in the outputs of a deepSTRAPP run carried out with run_deepSTRAPP_for_focal_time() when return_perm_data = TRUE so that the distributions of test stats computed across posterior samples are returned among the outputs under ⁠$STRAPP_results$perm_data_df⁠.

For multinomial data (categorical or biogeographic data with more than 2 states), it is possible to plot the histograms of post hoc pairwise tests. Set plot_posthoc_tests = TRUE to generate histograms for all the pairwise post hoc Dunn's tests across pairs of states. To achieve this, the deepSTRAPP_outputs input object must contain a ⁠$STRAPP_results$posthoc_pairwise_tests$perm_data_array⁠ element that summarizes test statistics computed across posterior samples for all pairwise post hoc tests. This is obtained from run_deepSTRAPP_for_focal_time() when setting both posthoc_pairwise_tests = TRUE to carry out post hoc tests, and return_perm_data = TRUE to record distributions of test statistics.

Alternatively, the main input deepSTRAPP_outputs can be the output of run_deepSTRAPP_over_time(), providing results of STRAPP tests over multiple time-steps. In this case, you must provide a focal_time to select the unique time-step used for plotting.

• return_perm_data must be set to TRUE so that the permutation data used to compute the tests are returned among the outputs under ⁠$STRAPP_results_over_time[[i]]$perm_data_df⁠.

• posthoc_pairwise_tests must be set to TRUE so that the permutation data used to perform the post hoc tests are also returned among the outputs under ⁠$STRAPP_results_over_time[[i]]$posthoc_pairwise_tests$perm_data_array⁠.

For plotting all time-steps at once, see plot_histograms_STRAPP_tests_over_time().

Value

By default, the function returns a list of classes gg and ggplot. This object is a ggplot that can be displayed on the console with print(output). It corresponds to the histogram being displayed on the console when the function is run, if display_plot = TRUE, and can be further modify for aesthetics using the ggplot2 grammar.

If using multinomial data and setting plot_posthoc_tests = TRUE, the function will return a list of objects. Each object is the ggplot associated with a pairwise post hoc test. To plot each histogram i individually, use print(output_list[[i]]). To plot all histograms at once in a multifaceted plot, as displayed on the console if display_plot = TRUE, use cowplot::plot_grid(plotlist = output_list).

Each plot also displays summary statistics for the STRAPP test associated with the data displayed.

• QX% = The quantile of null statistic distribution at the significance threshold used to define test significance. This is the value found on the red dashed line. The test will be considered significant (i.e., the null hypothesis is rejected) if this value is higher than zero (the black dashed line).

• P-value = The p-value of the STRAPP test which corresponds to the proportion of cases in which the statistic was lower than expected under the null hypothesis (i.e., the proportion of the histogram found below / on the left-side of the black dashed line).

• N = The number of test statistics in the null distribution.

If a PDF_file_path is provided, the function will also generate a PDF file of the plot. For post hoc tests, this will save the multifaceted plot.

Author(s)

Maël Doré

See Also

Associated functions in deepSTRAPP: run_deepSTRAPP_for_focal_time() plot_histograms_STRAPP_tests_over_time()

Examples

if (deepSTRAPP::is_dev_version())
{
 # ----- Example 1: Continuous trait ----- #

 # Load fake trait df
 data(Ponerinae_trait_tip_data, package = "deepSTRAPP")
 # Load phylogeny with old calibration
 data(Ponerinae_tree_old_calib, package = "deepSTRAPP")

 # Load the BAMM_object summarizing 1000 posterior samples of BAMM
 data(Ponerinae_BAMM_object_old_calib, package = "deepSTRAPP")
 ## This dataset is only available in development versions installed from GitHub.
 # It is not available in CRAN versions.
 # Use remotes::install_github(repo = "MaelDore/deepSTRAPP") to get the latest development version.

 ## Prepare trait data

 # Extract continuous trait data as a named vector
 Ponerinae_cont_tip_data <- setNames(object = Ponerinae_trait_tip_data$fake_cont_tip_data,
                                     nm = Ponerinae_trait_tip_data$Taxa)

 # Select a color scheme from lowest to highest values
 color_scale = c("darkgreen", "limegreen", "orange", "red")

 # Get Ancestral Character Estimates based on a Brownian Motion model
 # To obtain values at internal nodes
 Ponerinae_ACE <- phytools::fastAnc(tree = Ponerinae_tree_old_calib, x = Ponerinae_cont_tip_data)

  # (May take several minutes to run)
 # Run a Stochastic Mapping based on a Brownian Motion model
 # to interpolate values along branches and obtain a "contMap" object
 Ponerinae_contMap <- phytools::contMap(tree = Ponerinae_tree_old_calib,
                                        x = Ponerinae_cont_tip_data,
                                        res = 100, # Number of time steps
                                        plot = FALSE)
 # Plot contMap = stochastic mapping of continuous trait
 plot_contMap(contMap = Ponerinae_contMap,
              color_scale = color_scale)

 ## Set focal time to 10 Mya
 focal_time <- 10

 ## Run deepSTRAPP on net diversification rates for focal time = 10 Mya.

 Ponerinae_deepSTRAPP_cont_old_calib_10My <- run_deepSTRAPP_for_focal_time(
   contMap = Ponerinae_contMap,
   ace = Ponerinae_ACE,
   tip_data = Ponerinae_cont_tip_data,
   trait_data_type = "continuous",
   BAMM_object = Ponerinae_BAMM_object_old_calib,
   focal_time = focal_time,
   rate_type = "net_diversification",
   uncertainty_strategy = "rates_only",
   return_perm_data = TRUE)

 ## Explore output
 str(Ponerinae_deepSTRAPP_cont_old_calib_10My, max.level = 1)

 # ----- Plot histogram of STRAPP overall test results from run_deepSTRAPP_for_focal_time() ----- #

 histogram_ggplot <- plot_histogram_STRAPP_test_for_focal_time(
    deepSTRAPP_outputs = Ponerinae_deepSTRAPP_cont_old_calib_10My,
    display_plot = TRUE,
    # PDF_file_path = "./plot_STRAPP_histogram_10My.pdf",
    plot_posthoc_tests = FALSE)
 # Adjust aesthetics a posteriori
 histogram_ggplot_adj <- histogram_ggplot +
    ggplot2::theme(plot.title = ggplot2::element_text(color = "red", size = 15))
 print(histogram_ggplot_adj) 

 # ----- Plot histogram of STRAPP overall test results from run_deepSTRAPP_over_time() ----- #

 ## Load directly outputs from run_deepSTRAPP_over_time()
 Ponerinae_deepSTRAPP_cont_old_calib_0_40 <- readRDS(system.file("extdata",
    "Ponerinae_deepSTRAPP_cont_old_calib_0_40.rds", package = "deepSTRAPP"))
 ## This dataset is only available in development versions installed from GitHub.
 # It is not available in CRAN versions.
 # Use remotes::install_github(repo = "MaelDore/deepSTRAPP") to get the latest development version.

 # Select focal_time = 10My
 focal_time <- 10

 ## Plot histogram for overall test
 histogram_ggplot <- plot_histogram_STRAPP_test_for_focal_time(
    deepSTRAPP_outputs = Ponerinae_deepSTRAPP_cont_old_calib_0_40,
    focal_time = focal_time,
    display_plot = TRUE,
    # PDF_file_path = "./plot_STRAPP_histogram_10My.pdf",
    plot_posthoc_tests = FALSE)
 # Adjust aesthetics a posteriori
 histogram_ggplot_adj <- histogram_ggplot +
   ggplot2::theme(plot.title = ggplot2::element_text(color = "red", size = 15))
 print(histogram_ggplot_adj)

 # ----- Example 2: Categorical trait ----- #

 ## Load data

 # Load phylogeny
 data(Ponerinae_tree_old_calib, package = "deepSTRAPP")
 # Load trait df
 data(Ponerinae_trait_tip_data, package = "deepSTRAPP")

 # Load the BAMM_object summarizing 1000 posterior samples of BAMM
 data(Ponerinae_BAMM_object_old_calib, package = "deepSTRAPP")
 ## This dataset is only available in development versions installed from GitHub.
 # It is not available in CRAN versions.
 # Use remotes::install_github(repo = "MaelDore/deepSTRAPP") to get the latest development version.

 ## Prepare trait data

 # Extract categorical data with 3-levels
 Ponerinae_cat_3lvl_tip_data <- setNames(object = Ponerinae_trait_tip_data$fake_cat_3lvl_tip_data,
                                         nm = Ponerinae_trait_tip_data$Taxa)
 table(Ponerinae_cat_3lvl_tip_data)

 # Select color scheme for states
 colors_per_states <- c("forestgreen", "sienna", "goldenrod")
 names(colors_per_states) <- c("arboreal", "subterranean", "terricolous")

  # (May take several minutes to run)
 ## Produce densityMaps using stochastic character mapping based on an ARD Mk model
 Ponerinae_cat_3lvl_data_old_calib <- prepare_trait_data(
    tip_data = Ponerinae_cat_3lvl_tip_data,
    phylo = Ponerinae_tree_old_calib,
    trait_data_type = "categorical",
    colors_per_states = colors_per_states,
    evolutionary_models = "ARD", # Use default ARD model
    nb_simulations = 100, # Reduce number of simulations to save time
    seed = 1234, # Seet seed for reproducibility
    return_best_model_fit = TRUE,
    return_model_selection_df = TRUE,
    plot_map = FALSE) 

 # Load directly output
 data(Ponerinae_cat_3lvl_data_old_calib, package = "deepSTRAPP")

 ## Set focal time to 10 Mya
 focal_time <- 10

  # (May take several minutes to run)
 ## Run deepSTRAPP on net diversification rates for focal time = 10 Mya.

 Ponerinae_deepSTRAPP_cat_3lvl_old_calib_10My <- run_deepSTRAPP_for_focal_time(
    densityMaps = Ponerinae_cat_3lvl_data_old_calib$densityMaps,
    nb_simulations = 100,
    ace = Ponerinae_cat_3lvl_data_old_calib$ace,
    tip_data = Ponerinae_cat_3lvl_tip_data,
    trait_data_type = "categorical",
    BAMM_object = Ponerinae_BAMM_object_old_calib,
    focal_time = focal_time,
    rate_type = "net_diversification",
    uncertainty_strategy = "paired",
    posthoc_pairwise_tests = TRUE,
    return_perm_data = TRUE,
    return_updated_Maps = TRUE,
    return_updated_BAMM_object = TRUE)

 ## Explore output
 str(Ponerinae_deepSTRAPP_cat_3lvl_old_calib_10My, max.level = 1)

 # ------ Plot histogram of STRAPP overall test results ------ #

 histogram_ggplot <- plot_histogram_STRAPP_test_for_focal_time(
    deepSTRAPP_outputs = Ponerinae_deepSTRAPP_cat_3lvl_old_calib_10My,
    display_plot = TRUE,
    # PDF_file_path = "./plot_STRAPP_histogram_overall_test.pdf",
    plot_posthoc_tests = FALSE)
 # Adjust aesthetics a posteriori
 histogram_ggplot_adj <- histogram_ggplot +
    ggplot2::theme(plot.title = ggplot2::element_text(color = "red", size = 15))
 print(histogram_ggplot_adj)

 # ------ Plot histograms of STRAPP post hoc test results ------ #

 histograms_ggplot_list <- plot_histogram_STRAPP_test_for_focal_time(
    deepSTRAPP_outputs = Ponerinae_deepSTRAPP_cat_3lvl_old_calib_10My,
    display_plot = TRUE,
    # PDF_file_path = "./plot_STRAPP_histograms_posthoc_tests.pdf",
    plot_posthoc_tests = TRUE)
 # Plot all histograms one by one
 print(histograms_ggplot_list)
 # Plot all histograms on one faceted plot
 cowplot::plot_grid(plotlist = histograms_ggplot_list) 
}

---

Plot multiple histograms of STRAPP test statistics over time-steps

Description

Plot a histogram of the distribution of the test statistics obtained from a deepSTRAPP workflow carried out for each focal time in ⁠$time_steps⁠. Main input = output of a deepSTRAPP run over time using run_deepSTRAPP_over_time()).

Returns one histogram for overall tests for each focal time in ⁠$time_steps⁠. If plot_posthoc_tests = TRUE, it will return one faceted plot with a histogram per post hoc test for each focal time in ⁠$time_steps⁠.

If a PDF file path is provided in PDF_file_path, the plots will be saved directly in a PDF file, with one page per focal time in ⁠$time_steps⁠.

Usage

plot_histograms_STRAPP_tests_over_time(
  deepSTRAPP_outputs,
  display_plots = TRUE,
  plot_posthoc_tests = FALSE,
  PDF_file_path = NULL
)

Arguments

deepSTRAPP_outputs	List of elements generated with run_deepSTRAPP_over_time(), that summarize the results of multiple deepSTRAPP runs across ⁠$time_steps⁠. It needs to include the ⁠$STRAPP_results_over_time⁠ element with ⁠$perm_data_df⁠ obtained when setting both return_STRAPP_results = TRUE and return_perm_data = TRUE.
display_plots	Logical. Whether to display the histograms generated in the R console. Default is TRUE.
plot_posthoc_tests	Logical. For multinomial data only. Whether to plot the histograms for the overall Kruskal-Wallis test across all states (plot_posthoc_tests = FALSE), or plot the histograms for all the pairwise post hoc Dunn's tests across pairs of states (plot_posthoc_tests = TRUE). Time-steps at which the data does not yield more than two states/ranges will show a warning and generate no plot. Default is FALSE.
PDF_file_path	Character string. If provided, the plots will be saved in a unique PDF file following the path provided here. The path must end with ".pdf". Each page of the PDF corresponds to a focal time in ⁠$time_steps⁠.

Details

The main input deepSTRAPP_outputs is the typical output of run_deepSTRAPP_over_time(). It provides information on results of STRAPP tests performed over multiple time-steps.

Histograms are built based on the distribution of the test statistics. Such distributions are recorded in the outputs of a deepSTRAPP run carried out with run_deepSTRAPP_over_time() when return_STRAPP_results = TRUE AND return_perm_data = TRUE. The ⁠$STRAPP_results_over_time⁠ objects provided within the input are lists that must contain a ⁠$perm_data_df⁠ element that summarizes test statistics computed across posterior samples.

For multinomial data (categorical or biogeographic data with more than 2 states), it is possible to plot the histograms of post hoc pairwise tests. Set plot_posthoc_tests = TRUE to generate histograms for all the pairwise post hoc Dunn's tests across pairs of states. To achieve this, the ⁠$STRAPP_results_over_time⁠ objects must contain a ⁠$posthoc_pairwise_tests$perm_data_array⁠ element that summarizes test statistics computed across posterior samples for all pairwise post hoc tests. This is obtained from run_deepSTRAPP_over_time() when setting return_STRAPP_results = TRUE to return the STRAPP results, posthoc_pairwise_tests = TRUE to carry out post hoc tests, and return_perm_data = TRUE to record distributions of test statistics. Time-steps for which the data do not yield more than two states/ranges will show a warning and generate no plot.

Value

By default, the function returns a list of sub-lists of classes gg and ggplot ordered as in ⁠$time_steps⁠. Each sub-list corresponds to a ggplot for a given focal_time i that can be displayed on the console with print(output[[i]]). If display_plots = TRUE, the histograms are being displayed on the console one by one while generated.

If using multinomial data and set plot_posthoc_tests = TRUE, the function will return a list of sub-lists of objects ordered as in ⁠$time_steps⁠. Each sub-list is a list of the ggplots associated with pairwise post hoc tests carried out for a given focal_time. For a given focal_time i, to plot each histogram j individually, use print(output_list[[i]][[j]]). To plot all histograms of a given focal_time i at once in a multifaceted plot, as displayed sequentially on the console if display_plots = TRUE, use cowplot::plot_grid(plotlist = output_list[[i]]).

Each plot also displays summary statistics for the STRAPP test associated with the data displayed.

• QX% = The quantile of null statistic distribution at the significance threshold used to define test significance. This is the value found on the red dashed line. The test will be considered significant (i.e., the null hypothesis is rejected) if this value is higher than zero (the black dashed line).

• P-value = The p-value of the STRAPP test which corresponds to the proportion of cases in which the statistic was lower than expected under the null hypothesis (i.e., the proportion of the histogram found below / on the left-side of the black dashed line).

• N = The number of test statistics in the null distribution.

If a PDF_file_path is provided, the function will also generate a PDF file of the plots with one page per ⁠$time_steps⁠. For post hoc tests, this will save the multifaceted plots.

Author(s)

Maël Doré

See Also

Associated functions in deepSTRAPP: run_deepSTRAPP_over_time() plot_histogram_STRAPP_test_for_focal_time()

Examples

if (deepSTRAPP::is_dev_version())
{
 # ----- Example 1: Continuous trait ----- #

 # Load fake trait df
 data(Ponerinae_trait_tip_data, package = "deepSTRAPP")
 # Load phylogeny with old calibration
 data(Ponerinae_tree_old_calib, package = "deepSTRAPP")

 # Load the BAMM_object summarizing 1000 posterior samples of BAMM
 data(Ponerinae_BAMM_object_old_calib, package = "deepSTRAPP")
 ## This dataset is only available in development versions installed from GitHub.
 # It is not available in CRAN versions.
 # Use remotes::install_github(repo = "MaelDore/deepSTRAPP") to get the latest development version.

 ## Prepare trait data

 # Extract continuous trait data as a named vector
 Ponerinae_cont_tip_data <- setNames(object = Ponerinae_trait_tip_data$fake_cont_tip_data,
                                     nm = Ponerinae_trait_tip_data$Taxa)

 # Select a color scheme from lowest to highest values
 color_scale = c("darkgreen", "limegreen", "orange", "red")

 # Get Ancestral Character Estimates based on a Brownian Motion model
 # To obtain values at internal nodes
 Ponerinae_ACE <- phytools::fastAnc(tree = Ponerinae_tree_old_calib, x = Ponerinae_cont_tip_data)

  # (May take several minutes to run)
 # Run a Stochastic Mapping based on a Brownian Motion model
 # to interpolate values along branches and obtain a "contMap" object
 Ponerinae_contMap <- phytools::contMap(tree = Ponerinae_tree_old_calib,
                                        x = Ponerinae_cont_tip_data,
                                        res = 100, # Number of time steps
                                        plot = FALSE)
 # Plot contMap = stochastic mapping of continuous trait
 plot_contMap(contMap = Ponerinae_contMap,
              color_scale = color_scale)

 ## Set for time steps of 5 My. Will generate deepSTRAPP workflows for 0 to 40 Mya.
 # nb_time_steps <- 5
 time_step_duration <- 5
 time_range <- c(0, 40)

 ## Run deepSTRAPP on net diversification rates
 Ponerinae_deepSTRAPP_cont_old_calib_0_40 <- run_deepSTRAPP_over_time(
    contMap = Ponerinae_contMap,
    ace = Ponerinae_ACE,
    tip_data = Ponerinae_cont_tip_data,
    trait_data_type = "continuous",
    BAMM_object = Ponerinae_BAMM_object_old_calib,
    # nb_time_steps = nb_time_steps,
    time_range = time_range,
    time_step_duration = time_step_duration,
    uncertainty_strategy = "rates_only",
    return_perm_data = TRUE,
    extract_trait_data_melted_df = TRUE,
    extract_diversification_data_melted_df = TRUE,
    return_STRAPP_results = TRUE,
    return_updated_Maps = TRUE,
    return_updated_BAMM_objects = TRUE,
    verbose = TRUE,
    verbose_extended = TRUE) 

 ## Load directly trait data output
 Ponerinae_deepSTRAPP_cont_old_calib_0_40 <- readRDS(system.file("extdata",
    "Ponerinae_deepSTRAPP_cont_old_calib_0_40.rds", package = "deepSTRAPP"))
 ## This dataset is only available in development versions installed from GitHub.
 # It is not available in CRAN versions.
 # Use remotes::install_github(repo = "MaelDore/deepSTRAPP") to get the latest development version.

 ## Plot histograms of STRAPP overall test results
 # Tests are Spearman's rank correlation tests

 # Plot all histograms
 histogram_ggplots <- plot_histograms_STRAPP_tests_over_time(
    deepSTRAPP_outputs = Ponerinae_deepSTRAPP_cont_old_calib_0_40,
    display_plots = TRUE,
    # PDF_file_path = "./plot_STRAPP_histogram_overall_test.pdf",
    plot_posthoc_tests = FALSE)

 # Print histogram for time step 1 = 0 My
 print(histogram_ggplots[[1]])
 # Adjust aesthetics of plot for time step 1 a posteriori
 histogram_ggplot_adj <- histogram_ggplots[[1]] +
    ggplot2::theme(plot.title = ggplot2::element_text(color = "red", size = 15))
 print(histogram_ggplot_adj)

 # ----- Example 2: Categorical data ----- #

 ## Load data

 # Load trait df
 data(Ponerinae_trait_tip_data, package = "deepSTRAPP")
 # Load phylogeny
 data(Ponerinae_tree_old_calib, package = "deepSTRAPP")

 # Load the BAMM_object summarizing 1000 posterior samples of BAMM
 data(Ponerinae_BAMM_object_old_calib, package = "deepSTRAPP")
 ## This dataset is only available in development versions installed from GitHub.
 # It is not available in CRAN versions.
 # Use remotes::install_github(repo = "MaelDore/deepSTRAPP") to get the latest development version.

 ## Prepare trait data

 # Extract categorical data with 3-levels
 Ponerinae_cat_3lvl_tip_data <- setNames(object = Ponerinae_trait_tip_data$fake_cat_3lvl_tip_data,
                                         nm = Ponerinae_trait_tip_data$Taxa)
 table(Ponerinae_cat_3lvl_tip_data)

 # Select color scheme for states
 colors_per_states <- c("forestgreen", "sienna", "goldenrod")
 names(colors_per_states) <- c("arboreal", "subterranean", "terricolous")

  # (May take several minutes to run)
 ## Produce densityMaps using stochastic character mapping based on an ARD Mk model
 Ponerinae_cat_3lvl_data_old_calib <- prepare_trait_data(
    tip_data = Ponerinae_cat_3lvl_tip_data,
    phylo = Ponerinae_tree_old_calib,
    trait_data_type = "categorical",
    colors_per_levels = colors_per_states,
    evolutionary_models = "ARD",
    nb_simulations = 100,
    return_best_model_fit = TRUE,
    return_model_selection_df = TRUE,
    plot_map = FALSE) 

 # Load directly trait data output
 data(Ponerinae_cat_3lvl_data_old_calib, package = "deepSTRAPP")

 ## Set for time steps of 5 My. Will generate deepSTRAPP workflows for 0 to 40 Mya.
 # nb_time_steps <- 5
 time_step_duration <- 5
 time_range <- c(0, 40)

  # (May take several minutes to run)
 ## Run deepSTRAPP on net diversification rates across time-steps.
 Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40 <- run_deepSTRAPP_over_time(
    densityMaps = Ponerinae_cat_3lvl_data_old_calib$densityMaps,
    ace = Ponerinae_cat_3lvl_data_old_calib$ace,
    tip_data = Ponerinae_cat_3lvl_tip_data,
    nb_simulations = 100,
    trait_data_type = "categorical",
    BAMM_object = Ponerinae_BAMM_object_old_calib,
    # nb_time_steps = nb_time_steps,
    time_range = time_range,
    time_step_duration = time_step_duration,
    rate_type = "net_diversification",
    uncertainty_strategy = "paired",
    seed = 1234, # Set for reproducibility
    alpha = 0.10, # Select a generous level of significance for the sake of the example
    posthoc_pairwise_tests = TRUE,
    return_perm_data = TRUE,
    extract_trait_data_melted_df = TRUE,
    extract_diversification_data_melted_df = TRUE,
    return_STRAPP_results = TRUE,
    return_updated_Maps = TRUE,
    return_updated_BAMM_objects = TRUE,
    verbose = TRUE,
    verbose_extended = TRUE) 

 ## Load directly deepSTRAPP output
 Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40 <- readRDS(system.file("extdata",
    "Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40.rds", package = "deepSTRAPP"))
 ## This dataset is only available in development versions installed from GitHub.
 # It is not available in CRAN versions.
 # Use remotes::install_github(repo = "MaelDore/deepSTRAPP") to get the latest development version.

 ## Explore output
 str(Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40, max.level = 1)

 ## Plot histograms of STRAPP overall test results #
 # Tests are Kruskall-Wallis H tests when more than two states/ranges are present.
 # Tests are Mann-Whitney-Wilcoxon rank-sum tests when only two states/ranges are present.

 histogram_ggplots <- plot_histograms_STRAPP_tests_over_time(
    deepSTRAPP_outputs = Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40,
    display_plots = TRUE,
    # PDF_file_path = "./plot_STRAPP_histograms_overall_tests.pdf",
    plot_posthoc_tests = FALSE)

 # Print histogram for time step 1 = 0 My
 print(histogram_ggplots[[1]])
 # Adjust aesthetics of plot for time step 1 a posteriori
 histogram_ggplot_adj <- histogram_ggplots[[1]] +
     ggplot2::theme(plot.title = ggplot2::element_text(color = "red", size = 15))
 print(histogram_ggplot_adj)

 ## Plot histograms of STRAPP post hoc test results ------ #
 # Tests are Dunn's multiple comparison pairwise post hoc tests possible
 # only when more than two states/ranges are present.

 histograms_ggplots_list <- plot_histograms_STRAPP_tests_over_time(
     deepSTRAPP_outputs = Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40,
     display_plots = TRUE,
     # PDF_file_path = "./plot_STRAPP_histograms_posthoc_tests.pdf",
     plot_posthoc_tests = TRUE)

 # Print all histograms for time step 1 (= 0 My) one by one
 print(histograms_ggplots_list[[1]])
 # Plot all histograms for time step 1 (= 0 My) on one faceted plot
 cowplot::plot_grid(plotlist = histograms_ggplots_list[[1]])
}

---