# deepSTRAPP: Categorical trait data

#### Maël Doré

#### 2026-01-13

This is a simple example that shows how deepSTRAPP can be used to test for differences in diversification rates between multiple (three) categorical states along evolutionary times. It presents the main functions in a typical deepSTRAPP workflow.

For an example with binary data (2 levels), please see the example in the Main tutorial: `vignette("main_tutorial")`.

For an example with continuous data, see this vignette: `vignette("deepSTRAPP_continuous_data")`

Please note that the trait data and phylogeny calibration used in this example are NOT valid biological data. They were modified in order to provide results illustrating the usefulness of deepSTRAPP.

    # ------ Step 0: Load data ------ #

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

    dim(Ponerinae_trait_tip_data)
    View(Ponerinae_trait_tip_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)
    # Here, data represent three types of habitats
    table(Ponerinae_cat_3lvl_tip_data)

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

    ## Load phylogeny with old time-calibration
    data(Ponerinae_tree_old_calib, package = "deepSTRAPP")

    plot(Ponerinae_tree_old_calib)
    ape::Ntip(Ponerinae_tree_old_calib) == length(Ponerinae_cat_3lvl_tip_data)

    ## Check that trait data and phylogeny are named and ordered similarly
    all(names(Ponerinae_cat_3lvl_tip_data) == Ponerinae_tree_old_calib$tip.label)

    ## Reorder trait data as in phylogeny
    Ponerinae_cat_3lvl_tip_data <- Ponerinae_cat_3lvl_tip_data[match(x = Ponerinae_tree_old_calib$tip.label,
                                                                     table = names(Ponerinae_cat_3lvl_tip_data))]

    ## Plot data on tips for visualization
    pdf(file = "./Ponerinae_cat_3lvl_data_old_calib_on_phylo.pdf", width = 20, height = 150)

    # Set plotting parameters
    old_par <- par(no.readonly = TRUE)
    par(mar = c(0.1,0.1,0.1,0.1), oma = c(0,0,0,0)) # bltr

    # Graph presence/absence using plotTree.datamatrix
    range_map <- phytools::plotTree.datamatrix(
      tree = Ponerinae_tree_old_calib,
      X = as.data.frame(Ponerinae_cat_3lvl_tip_data),
      fsize = 0.7, yexp = 1.1,
      header = TRUE, xexp = 1.25,
      colors = colors_per_states)

    # Get plot info in "last_plot.phylo"
    plot_info <- get("last_plot.phylo", envir=.PlotPhyloEnv)

    # Add time line

    # Extract root age
    root_age <- max(phytools::nodeHeights(Ponerinae_tree_old_calib))

    # Define ticks
    # ticks_labels <- seq(from = 0, to = 100, by = 20)
    ticks_labels <- seq(from = 0, to = 120, by = 20)
    axis(side = 1, pos = 0, at = (-1 * ticks_labels) + root_age, labels = ticks_labels, cex.axis = 1.5)
    legend(x = root_age/2,
           y = 0 - 5, adj = 0,
           bty = "n", legend = "", title = "Time  [My]", title.cex = 1.5)

    # Add a legend
    legend(x = plot_info$x.lim[2] - 10,
           y = mean(plot_info$y.lim),
           # adj = c(0,0),
           # x = "topleft",
           legend = c("Absence", "Presence"),
           pch = 22, pt.bg = c("white","gray30"), pt.cex =  1.8,
           cex = 1.5, bty = "n")

    dev.off()

    # Reset plotting parameters
    par(old_par)

    ## Inputs needed for Step 1 are the tip_data (Ponerinae_cat_3lvl_tip_data) and the phylogeny
    ## (Ponerinae_tree_old_calib), and optionally, a color scheme (colors_per_states).

    # ------ Step 1: Prepare trait data ------ #

    ## Goal: Map trait evolution on the time-calibrated phylogeny

    # 1.1/ Fit evolutionary models to trait data using Maximum Likelihood.
    # 1.2/ Select the best fitting model comparing AICc.
    # 1.3/ Infer ancestral characters estimates (ACE) at nodes.
    # 1.4/ Run stochastic mapping simulations to generate evolutionary histories
    #      compatible with the best model and inferred ACE.
    # 1.5/ Infer ancestral states along branches.
    #  - Compute posterior frequencies of each state to produce a `densityMap` for each state.

    library(deepSTRAPP)

    # All these actions are performed by a single function: deepSTRAPP::prepare_trait_data()
    ?deepSTRAPP::prepare_trait_data()

    # Run prepare_trait_data with default options
    # For categorical trait, an ARD model is assumed by default.
    Ponerinae_trait_object <- 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,
       nb_simulations = 100, # Reduce number of simulations to save time
       seed = 1234) # Set seed for reproducibility

    # Explore output
    str(Ponerinae_trait_object, 1)

    # Extract the densityMaps representing the posterior probabilities of states on the phylogeny
    Ponerinae_densityMaps <- Ponerinae_trait_object$densityMaps

    # Plot ancestral states as a single continuously mapped phylogeny overlaying
    # all state posterior probabilities
    plot_densityMaps_overlay(Ponerinae_densityMaps,
                             colors_per_levels = colors_per_states)

    # Plot posterior probabilities of each state on an independent densityMap
    # Plot densityMap for state = "arboreal"
    plot(Ponerinae_densityMaps[[1]])
    # Plot densityMap for state = "subterranean"
    plot(Ponerinae_densityMaps[[2]])
    # Plot densityMap for state = "terricolous"
    plot(Ponerinae_densityMaps[[3]])

    # Extract the Ancestral Character Estimates (ACE) = trait values at nodes
    Ponerinae_ACE <- Ponerinae_trait_object$ace
    head(Ponerinae_ACE)

    ## Inputs needed for Step 2 are the densityMaps, and optionally, the tip_data
    ## (Ponerinae_cat_3lvl_tip_data), and the ACE (Ponerinae_ACE)

    # ------ Step 2: Prepare diversification data ------ #

    ## Goal: Map evolution of diversification rates and regime shifts on the time-calibrated phylogeny

    # Run a BAMM (Bayesian Analysis of Macroevolutionary Mixtures)

    # You need the BAMM C++ program installed in your machine to run this step.
    # See the BAMM website: http://bamm-project.org/ and the companion R package [BAMMtools].

    # 2.1/ Set BAMM - Record BAMM settings and generate all input files needed for BAMM.
    # 2.2/ Run BAMM - Run BAMM and move output files in dedicated directory.
    # 2.3/ Evaluate BAMM - Produce evaluation plots and ESS data.
    # 2.4/ Import BAMM outputs - Load `BAMM_object` in R and subset posterior samples.
    # 2.5/ Clean BAMM files - Remove files generated during the BAMM run.

    # All these actions are performed by a single function: deepSTRAPP::prepare_diversification_data()
    ?deepSTRAPP::prepare_diversification_data()

    # Run BAMM workflow with deepSTRAPP
    ## This step is time-consuming. You can skip it and load directly the result if needed
    Ponerinae_BAMM_object_old_calib <- prepare_diversification_data(
       BAMM_install_directory_path = "./software/bamm-2.5.0/", # To adjust to your own path to BAMM
       phylo = Ponerinae_tree_old_calib,
       prefix_for_files = "Ponerinae_old_calib",
       seed = 1234, # Set seed for reproducibility
       numberOfGenerations = 10^7, # Set high for optimal run, but will take a long time
       BAMM_output_directory_path =  "./BAMM_outputs/")

    # Load directly the result
    data(Ponerinae_BAMM_object_old_calib)
    # 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_BAMM_object_old_calib, 1)
    # Record the regime shift events and macroevolutionary regimes parameters across posterior samples
    str(Ponerinae_BAMM_object_old_calib$eventData, 1)
    # Mean speciation rates at tips aggregated across all posterior samples
    head(Ponerinae_BAMM_object_old_calib$meanTipLambda)
    # Mean extinction rates at tips aggregated across all posterior samples
    head(Ponerinae_BAMM_object_old_calib$meanTipMu)

    # Plot mean net diversification rates and regime shifts on the phylogeny
    plot_BAMM_rates(Ponerinae_BAMM_object_old_calib,
                    labels = FALSE, legend = TRUE)

    ## Input needed for Step 3 is the BAMM_object (Ponerinae_BAMM_object)

    # ------ Step 3: Run a deepSTRAPP workflow ------ #

    ## Goal: Extract traits, diversification rates and regimes at a given time in the past
    ## to test for differences with a STRAPP test

    # 3.1/ Extract trait data at a given time in the past ('focal_time')
    # 3.2/ Extract diversification rates and regimes at a given time in the past ('focal_time')
    # 3.3/ Compute STRAPP test
    # 3.4/ Repeat previous actions for many timesteps along evolutionary time

    # Because we have three levels as trait data, two types of tests can be performed:
    #  - Overall Kruskal-Wallis tests that test for rate differences across all states at once.
    #  - post hoc pairwise Dunn's tests that test for rate differences between pairs of states.
    # Here, we select 'posthoc_pairwise_tests = TRUE' to conduct post hoc pairwise tests
    # in addition to overall Kruskal-Wallis tests.

    # All these actions are performed by a single function:
    #  For a single 'focal_time': deepSTRAPP::run_deepSTRAPP_for_focal_time()
    #  For multiple 'time_steps': deepSTRAPP::run_deepSTRAPP_over_time()
    ?deepSTRAPP::run_deepSTRAPP_for_focal_time()
    ?deepSTRAPP::run_deepSTRAPP_over_time()

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

    # Run deepSTRAPP on net diversification rates
    ## This step is time-consuming. You can skip it and load directly the result if needed
    Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40 <- run_deepSTRAPP_over_time(
        densityMaps = Ponerinae_densityMaps,
        ace = Ponerinae_ACE,
        tip_data = Ponerinae_cat_3lvl_data,
        trait_data_type = "categorical",
        BAMM_object = Ponerinae_BAMM_object_old_calib,
        time_range = time_range,
        time_step_duration = time_step_duration,
        seed = 1234, # Set seed for reproducibility
        alpha = 0.10, # Set significance threshold to use for tests
        posthoc_pairwise_tests = TRUE, # To run pairwise posthoc tests between pairs of states
        # Needed to obtain STRAPP stats and plot evaluation histograms (See 4.2)
        return_perm_data = TRUE,
        # Needed to get trait data and plot rates through time (See 4.3)
        extract_trait_data_melted_df = TRUE,
        # Needed to get diversification data and plot rates through time (See 4.3)
        extract_diversification_data_melted_df = TRUE,
        # Needed to obtain STRAPP stats and plot evaluation histograms (See 4.2)
        return_STRAPP_results = TRUE,
        # Needed to plot updated densityMaps (See 4.4)
        return_updated_trait_data_with_Map = TRUE,
        # Needed to map diversification rates on updated phylogenies (See 4.5)
        return_updated_BAMM_object = TRUE,
        verbose = TRUE,
        verbose_extended = TRUE)

    # Load the deepSTRAPP output summarizing results for 0 to 40 My
    data(Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40, 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)

    # See next step for how to generate plots from those outputs

    # Display test summaries
    # Can be passed down to [deepSTRAPP::plot_STRAPP_pvalues_over_time()] to generate a plot
    # showing the evolution of the test results across time

    # For overall Kruskal-Wallis tests over time-steps
    Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40$pvalues_summary_df
    # For posthoc pairwise Dunn's tests over time-steps
    Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40$pvalues_summary_df_for_posthoc_pairwise_tests

    # Access STRAPP test results
    # Can be passed down to [deepSTRAPP::plot_histograms_STRAPP_tests_over_time()] to generate plot
    # showing the null distribution of the test statistics

    # For overall Kruskal-Wallis tests over time-steps
    str(Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40$STRAPP_results, max.level = 2)
    # For posthoc pairwise Dunn's tests over time-steps
    # Results are found in the '$posthoc_pairwise_tests' element of each 'STRAPP_result'.
    str(lapply(X = Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40$STRAPP_results,
               FUN = function (x) { x$posthoc_pairwise_tests } ), max.level = 3)

    # Access trait data in a melted data.frame
    head(Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40$trait_data_df_over_time)
    # Access the diversification data in a melted data.frame
    head(Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40$diversification_data_df_over_time)
    # Both can be passed down to [deepSTRAPP::plot_rates_through_time()] to generate a plot
    # showing the evolution of diversification rates though time in relation to trait values

    # Access updated densityMaps for each focal time
    # Can be used to plot densityMaps with branch cut-off at focal time
    # with [deepSTRAPP::plot_densityMaps_overlay()]
    str(Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40$updated_trait_data_with_Map_over_time, max.level = 2)

    # Access updated BAMM_object for each focal time
    # Can be used to map rates and regime shifts on phylogeny with branch cut-off
    # at focal time with [deepSTRAPP::plot_BAMM_rates()]
    str(Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40$updated_BAMM_objects_over_time, max.level = 2)

    ## Input needed for Step 4 is the deepSTRAPP object (Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40)

    # ------ Step 4: Plot results ------ #

    ## Goal: Summarize the outputs in meaningful plots

    # 4.1/ Plot evolution of STRAPP tests p-values through time
    # 4.2/ Plot histogram of STRAPP test stats
    # 4.3/ Plot evolution of rates through time in relation to trait values
    # 4.4/ Plot rates vs. states across branches for a given 'focal_time'
    # 4.5/ Plot updated densityMaps mapping trait evolution for a given 'focal_time'
    # 4.6/ Plot updated diversification rates and regimes for a given 'focal_time'
    # 4.7/ Combine 4.5 and 4.6 to plot both mapped phylogenies with trait evolution (4.5)
    #      and diversification rates and regimes (4.6).

    # Each plot is achieved through a dedicated function

    # Load the deepSTRAPP output summarizing results for 0 to 40 My
    data(Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40, 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.

    ### 4.1/ Plot evolution of STRAPP tests p-values through time ####

    # ?deepSTRAPP::plot_STRAPP_pvalues_over_time()

    ## 4.1.1/ Plot results of overall Kruskal-Wallis tests over time

    deepSTRAPP::plot_STRAPP_pvalues_over_time(
       deepSTRAPP_outputs = Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40,
       alpha = 0.10)

    # This is the main output of deepSTRAPP. They show the evolution of the significance of
    # the STRAPP tests over time.
    # Here, overall Kruskal-Wallis tests for rate difference across all states (i.e., habitats) are shown.
    # This example highlights the importance of deepSTRAPP as the significance of
    # STRAPP tests change over time.
    # Differences in net diversification rates are not significant in the present
    # (assuming a significant threshold of alpha = 0.10).
    # Meanwhile, rates are significantly different in the past between 5 My to 15 My (the green area).
    # This result supports the idea that differences in biodiversity across habitats
    # (i.e., "arboreal" vs. , "subterranean" vs. "terricolous" ants) can be explained
    # by differences of diversification rates that was detected in the past. Without use of deepSTRAPP,
    # this conclusion would not have been supported by current diversification rates alone.

    # Note: This is NOT true ecological data. It is not a valid scientific result,
    # but an illustration of the use of deepSTRAPP.

    # A next step is to look in details into rate differences across pairs of states (i.e., habitats).
    # For this, we can plot the results of the post hoc pairwise tests.

    ## 4.1.2/ Plot results of posthoc pairwise Dunn's tests over time

    deepSTRAPP::plot_STRAPP_pvalues_over_time(
       deepSTRAPP_outputs = Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40,
       plot_posthoc_tests = TRUE) # To plot results of post hoc pairwise tests instead

    # Here, post hoc pairwise Dunn's tests for rate difference between pairs of states are shown.
    # These results show that differences in rates were only detected between "arboreal"
    # and "terricolous" ants between 2 My to 15 My (the green area), providing more detailed insights on
    # how type of habitats may affect diversification rates.
    # Note: This is NOT true ecological data. It is not a valid scientific result,
    # but an illustration of the use of deepSTRAPP.
    # This highlights the critical use of deepSTRAPP in revealing differences in diversification rates
    # occurring in the past, that may drive current biodiversity patterns.

    ### 4.2/ Plot histogram of STRAPP test stats ####

    # Plot an histogram of the distribution of the test statistics used to assess
    # the significance of STRAPP tests
      #  For a single 'focal_time': deepSTRAPP::plot_histogram_STRAPP_test_for_focal_time()
      #  For multiple 'time_steps': deepSTRAPP::plot_histograms_STRAPP_tests_over_time()

    # ?deepSTRAPP::plot_histogram_STRAPP_test_for_focal_time
    # ?deepSTRAPP::plot_histograms_STRAPP_tests_over_time

    ## These functions are used to provide visual illustration of the results of each STRAPP test.
    # They can be used to complement the provision of the statistical results summarized in Step 3.

    # Display the time-steps
    Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40$time_steps

    ## 4.2.1/ Plot results from overall Kruskal-Wallis tests across all states ####

    # Plot the histogram of overall Kruskal-Wallis stats for time-step n°3 = 10 My
    plot_histogram_STRAPP_test_for_focal_time(
       deepSTRAPP_outputs = Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40,
       focal_time = 10)

    # The black line represents the expected value under the null hypothesis H0
    #   => Δ Kruskal-Wallis H-stat = 0.
    # The histogram shows the distribution of the test statistics as observed across
    # the 1000 posterior samples from BAMM.
    # The red line represents the significance threshold for which 90% of the observed data
    # exhibited a higher value than expected.
    # Since this red line is below the null expectation (quantile 10% = 6.942),
    # the test is significant for a value of alpha = 0.10.
    # However, this significance must be discussed in regards to the relatively generous
    # significance threshold chosen here (alpha = 0.10).

    # Plot the histograms of overall Kruskal-Wallis stats for all time-steps
    plot_histograms_STRAPP_tests_over_time(
       deepSTRAPP_outputs = Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40)

    ## 4.2.2/ Plot results from posthoc pairwise Dunn's tests between pairs of states ####

    # Plot the histogram of posthoc pairwise Dunn's stats for time-step n°3 = 10 My
    plot_histogram_STRAPP_test_for_focal_time(
       deepSTRAPP_outputs = Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40,
       plot_posthoc_tests = TRUE, # To plot results of post hoc pairwise tests instead
       focal_time = 10)

    # Each facet represent a pairwise post hoc test conducted across a given pair of states.
    # In each facet, the black line represents the expected value under the null hypothesis H0
    #  => Δ Dunn's Z-stat = 0.
    # The red line represents the significance threshold for which 90% of the observed data
    # exhibited a higher value than expected.
    # This red line is below the null expectation for the "arboreal != subterranean" and
    # "subterranean != terricolous" pairs. This means the test is not significant for these pairs of habitats.
    # The red line is above the null expectation for the "arboreal != terricolous" pair
    # (Q10% = 1.695, p = 0.025). This means the test is significant for this pair of habitat.
    # This is the pair that is driving the significance detected in the previous plot
    # when looking at differences across all habitats.
    # This significance must still be discussed in regards to the relatively generous
    # significance threshold chosen here (alpha = 0.10).

    # Plot the histograms of posthoc pairwise Dunn's stats for all time-steps
    plot_histograms_STRAPP_tests_over_time(
       deepSTRAPP_outputs = Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40,
       plot_posthoc_tests = TRUE) # To plot results of post hoc pairwise tests instead

    ### 4.3/ Plot evolution of rates through time ~ trait data ####

    # ?deepSTRAPP::plot_rates_through_time()

    # Generate ggplot
    plot_rates_through_time(deepSTRAPP_outputs = Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40,
                            colors_per_levels = colors_per_states,
                            plot_CI = TRUE)

    # This plot helps to visualize how differences in rates evolved over time.

    # You can see that both type of ants "arboreal" and "terricolous" had fairly different rates over time,
    # with differences detected as significant between 2 to 15 My.
    # Meanwhile, "subterranean" ants exhibited intermediate diversification levels.
    # This plot, alongside results of deepSTRAPP, supports the Diversification Rate Hypothesis in showing
    # how "terricolous" ant lineages may have accumulated faster, especially between 2 to 15 My.
    # It hints that "terricolous" ant lineages are fairly recent as no lineage in this state/habitat
    # is inferred to have existed before 25 Mya.
    # The larger uncertainty across estimates of diversification rates for "terricolous" ant lineages
    # also hints at their relatively lower number due to their recent emergence.

    # Note: This is NOT true ecological data. It is not a valid scientific result,
    # but an illustration of the use of deepSTRAPP.

    ### 4.4/ Plot rates vs. states across branches for a given focal time ####

    # ?deepSTRAPP::plot_rates_vs_trait_data_for_focal_time()
    # ?deepSTRAPP::plot_rates_vs_trait_data_over_time()

    # This plot help to visualize differences in rates vs. states across all branches
    # found at specific time-steps (i.e., 'focal_time').

    # Generate ggplot for time = 10 My
    plot_rates_vs_trait_data_for_focal_time(
       deepSTRAPP_outputs = Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40,
       focal_time = 10,
       colors_per_levels = colors_per_states)

    # Here we focus on T = 10 My to highlight the differences detected in the previous steps.
    # You can see that "terricolous" ants tend to have higher rates than "subterranean" ants,
    # who tends to have higher rates than "arboreal" ants, at this time-step.
    # This plot, alongside other results of deepSTRAPP, supports the Diversification Rate Hypothesis in showing
    # how "terricolous" ant lineages may have accumulated faster, especially between 5 to 15 My.
    # Additionally, the plot displays summary statistics for the STRAPP test associated with the data shown:
    #   * An observed statistic computed across the mean rates and trait states (i.e., habitats) shown on the plot.
    #     Here, H-stat obs = 374.82. Please note that this is not the statistic of the STRAPP test itself,
    #     which is conducted across all BAMM posterior samples.
    #   * The quantile of null statistic distribution at the significant threshold used to define test significance.
    #     The test will be considered significant (i.e., the null hypothesis is rejected)
    #     if this value is higher than zero, as shown on the histogram in Section 4.2.
    #     Here, Q10% = 6.942, so the test is significant (according to this significance threshold).
    #   * The p-value of the associated STRAPP test. Here, p = 0.071.

    # Plot rates vs. trait data for all time-steps
    plot_rates_vs_trait_data_over_time(
       deepSTRAPP_outputs = Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40,
       colors_per_levels = colors_per_states)

    ### 4.5/ Plot updated densityMaps mapping trait evolution for a given 'focal_time' ####

    # ?deepSTRAPP::plot_densityMaps_overlay()

    ## These plots help to visualize the evolution of states across the phylogeny,
    ## and to focus on tip values at specific time-steps.

    # Display the time-steps
    Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40$time_steps

    ## The next plot shows the evolution of states across the whole phylogeny (100-0 My).

    # Plot initial densityMaps (t = 0)
    densityMaps_0My <- Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40$updated_trait_data_with_Map_over_time[[1]]
    plot_densityMaps_overlay(densityMaps_0My$densityMaps,
                             colors_per_levels = colors_per_states,
                             fsize = 0.1) # Reduce tip label size
    title(main = "Trait evolution for 100-0 My")

    # It highlights the relatively recent emergence of "terricolous" ants (in this fake illustrative dataset),
    # where no lineages exhibit this state in deep times.

    ## The next plot shows the evolution of states from root to 10 Mya (100-10 My).

    # Plot updated densityMaps for time-step n°3 = 10 My
    densityMaps_10My <- Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40$updated_trait_data_with_Map_over_time[[3]]
    plot_densityMaps_overlay(densityMaps_10My$densityMaps,
                             colors_per_levels = colors_per_states,
                             fsize = 0.1) # Reduce tip label size
    title(main = "Trait evolution for 100-10 My")

    ## The next plot shows the evolution of states from root to 40 Mya (100-40 My).

    # Plot updated densityMaps for time-step n°9 = 40 My
    densityMaps_40My <- Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40$updated_trait_data_with_Map_over_time[[9]]
    plot_densityMaps_overlay(densityMaps = densityMaps_40My$densityMaps,
                             colors_per_levels = colors_per_states,
                             fsize = 0.2) # Reduce tip label size
    title(main = "Trait evolution for 100-40 My")

    # In this simulated illustrative dataset, no ant lineages are inferred in "terricolous" habitats 40 Mya.

    ### 4.6/ Plot updated diversification rates and regimes for a given 'focal_time' ####

    # ?deepSTRAPP::plot_BAMM_rates()

    ## These plots help to visualize the evolution of diversification rates across the phylogeny,
    ## and to focus on tip values at specific time-steps.

    # Display the time-steps
    Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40$time_steps

    # Extract root age
    root_age <- max(phytools::nodeHeights(Ponerinae_tree_old_calib)[,2])

    ## The next plot shows the evolution of diversification rates across the whole phylogeny (100-0 My).

    # Plot diversification rates on initial phylogeny (t = 0)
    BAMM_map_0My <- Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40$updated_BAMM_objects_over_time[[1]]
    plot_BAMM_rates(BAMM_map_0My, labels = FALSE, par.reset = FALSE)
    abline(v = root_age - 10, col = "red", lty = 2) # Show where the phylogeny will be cut at 10 Mya
    abline(v = root_age - 40, col = "red", lty = 2) # Show where the phylogeny will be cut at 40 Mya
    title(main = "BAMM rates for 100-0 My")

    ## The next plot shows the evolution of diversification rates from root to 10 Mya (100-10 My).

    # Plot diversification rates on updated phylogeny for time-step n°3 = 10 My
    BAMM_map_10My <- Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40$updated_BAMM_objects_over_time[[3]]
    plot_BAMM_rates(BAMM_map_10My, labels = FALSE,
                    colorbreaks = BAMM_map_10My$initial_colorbreaks$net_diversification)
    title(main = "BAMM rates for 100-10 My")

    ## The next plot shows the evolution of diversification rates from root to 40 Mya (100-40 My).

    # Plot diversification rates on updated phylogeny for time-step n°9 = 40 My
    BAMM_map_40My <- Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40$updated_BAMM_objects_over_time[[9]]
    plot_BAMM_rates(BAMM_map_40My, labels = FALSE,
                    colorbreaks = BAMM_map_40My$initial_colorbreaks$net_diversification)
    title(main = "BAMM rates for 100-40 My")

    ### 4.7/ Plot both trait evolution and diversification rates and regimes updated for a given 'focal_time' ####

    # ?deepSTRAPP::plot_traits_vs_rates_on_phylogeny_for_focal_time()

    ## These plots help to visualize simultaneously the evolution of trait and diversification rates
    ## across the phylogeny, and to focus on tip values at specific time-steps.

    # Display the time-steps
    Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40$time_steps

    ## The next plot shows the evolution of states and rates across the whole phylogeny (100-0 My).

    # Plot both mapped phylogenies in the present (t = 0)
    plot_traits_vs_rates_on_phylogeny_for_focal_time(
      deepSTRAPP_outputs = Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40,
      focal_time = 0,
      ftype = "off", lwd = 0.7,
      colors_per_levels = colors_per_states,
      labels = FALSE, legend = FALSE,
      par.reset = FALSE)

    ## The next plot shows the evolution of states and rates from root to 10 Mya (100-10 My).

    # Plot both mapped phylogenies for time-step n°3 = 10 My
    plot_traits_vs_rates_on_phylogeny_for_focal_time(
      deepSTRAPP_outputs = Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40,
      focal_time = 10,
      ftype = "off", lwd = 1.2,
      colors_per_levels = colors_per_states,
      labels = FALSE, legend = FALSE,
      par.reset = FALSE)

    ## The next plot shows the evolution of states and rates from root to 40 Mya (100-40 My).

    # Plot both mapped phylogenies for time-step n°9 = 40 My
    plot_traits_vs_rates_on_phylogeny_for_focal_time(
      deepSTRAPP_outputs = Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40,
      focal_time = 40,
      ftype = "off", lwd = 1.2,
      colors_per_levels = colors_per_states,
      labels = FALSE, legend = FALSE,
      par.reset = FALSE)