## ----set_options, include = FALSE---------------------------------------------
knitr::opts_chunk$set(
  eval = FALSE, # Chunks of codes will not be evaluated by default
  collapse = TRUE,
  comment = "#>",
  fig.width = 7, fig.height = 5, # Set device size at rendering time (when plots are generated)
  fig.align = "center"
)

## ----setup, eval = TRUE, include = FALSE--------------------------------------
library(deepSTRAPP)

is_dev_version <- function (pkg = "deepSTRAPP")
{
  # # Check if ran on CRAN
  # not_cran <- identical(Sys.getenv("NOT_CRAN"), "true") # || interactive()

  # Version number check
  version <- tryCatch(as.character(utils::packageVersion(pkg)), error = function(e) "")
  dev_version <- grepl("\\.9000", version)

  # not_cran || dev_version
  
  return(dev_version)
}


## ----adjust_dpi_CRAN, include = FALSE, eval = !is_dev_version()---------------
knitr::opts_chunk$set(
  dpi = 50   # Lower DPI to save space
)

## ----adjust_dpi_dev, include = FALSE, eval = is_dev_version()-----------------
# knitr::opts_chunk$set(
#   dpi = 72   # Default DPI for the dev version
# )

## ----load_data_cat_3lvl-------------------------------------------------------
# # ------ 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).
# 

## ----load_data_cat_3lvl_eval, eval = TRUE, echo = FALSE-----------------------
## Select color scheme for states
colors_per_states <- c("forestgreen", "sienna", "goldenrod")
names(colors_per_states) <- c("arboreal", "subterranean", "terricolous")

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

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

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

## ----plot_pvalues_cat_3lvl----------------------------------------------------
# # ------ 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.
# 

## ----plot_pvalues_cat_3lvl_eval_dev, eval = is_dev_version(), echo = FALSE----
# 
# # Load the deepSTRAPP output summarizing results for 0 to 40 My
# data(Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40, package = "deepSTRAPP")
# 
# # Produce the results of overall Kruskal-Wallis tests over time
# ggplot_STRAPP_pvalues <- deepSTRAPP::plot_STRAPP_pvalues_over_time(
#    deepSTRAPP_outputs = Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40,
#    alpha = 0.10, display_plot = FALSE)
# # Adjust main title size
# ggplot_STRAPP_pvalues <- ggplot_STRAPP_pvalues +
#   ggplot2::theme(plot.title = ggplot2::element_text(size = 18))
# # Print plot
# print(ggplot_STRAPP_pvalues)
# 
# # Produce the results of posthoc pairwise Dunn's tests over time
# ggplot_STRAPP_pvalues_posthoc <- 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
#    display_plot = FALSE)
# # Adjust main title size + legend
# ggplot_STRAPP_pvalues_posthoc <- ggplot_STRAPP_pvalues_posthoc +
#   ggplot2::theme(
#     plot.title = ggplot2::element_text(size = 18),
#     legend.title = ggplot2::element_text(size  = 12),
#     legend.position.inside = c(0.30, 0.40),
#     legend.text = ggplot2::element_text(size = 9))
# # Print plot
# print(ggplot_STRAPP_pvalues_posthoc)
# 

## ----plot_pvalues_cat_3lvl_eval_CRAN, eval = !is_dev_version(), echo = FALSE, out.width = "100%"----

# Plot pre-rendered graph
knitr::include_graphics("figures/1.2_deepSTRAPP_categorical_3lvl_data_4.1_plot_pvalues_1.PNG")
knitr::include_graphics("figures/1.2_deepSTRAPP_categorical_3lvl_data_4.1_plot_pvalues_2.PNG")


## ----plot_histogram_STRAPP_tests_overall_cat_3lvl-----------------------------
# ### 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
# 

## ----plot_histogram_STRAPP_tests_cat_3lvl_eval_dev, fig.width = 8.5, fig.height = 6, out.width = "100%", eval = is_dev_version(), echo = FALSE----
# # Plot the histogram of test 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)
# 
# # Plot the histogram of test stats for time-step n°3 = 10 My
# ggplot_STRAPP_pvalues_posthoc <- 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)

## ----plot_histogram_STRAPP_tests_cat_3lvl_eval_CRAN, eval = !is_dev_version(), echo = FALSE, out.width = "100%"----

# Plot pre-rendered graph
knitr::include_graphics("figures/1.2_deepSTRAPP_categorical_3lvl_data_4.2_plot_histograms_1.PNG")
knitr::include_graphics("figures/1.2_deepSTRAPP_categorical_3lvl_data_4.2_plot_histograms_2.PNG")


## ----plot_rates_through_time_cat_3lvl-----------------------------------------
# ### 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.
# 

## ----plot_rates_through_time_cat_3lvl_eval_dev, fig.width = 8.5, out.width = "100%", eval = is_dev_version(), echo = FALSE----
# # Produce RTT plot
# ggplot_RTT_list <- plot_rates_through_time(deepSTRAPP_outputs = Ponerinae_deepSTRAPP_cat_3lvl_old_calib_0_40,
#                         colors_per_levels = colors_per_states,
#                         plot_CI = TRUE, display_plot = FALSE)
# # Adjust title size
# ggplot_RTT <- ggplot_RTT_list$rates_TT_ggplot +
#   ggplot2::theme(plot.title = ggplot2::element_text(size = 18),
#                  axis.title = ggplot2::element_text(size = 16))
# # Print plot
# print(ggplot_RTT)
# 

## ----plot_rates_through_time_cat_3lvl_eval_CRAN, eval = !is_dev_version(), echo = FALSE, out.width = "100%"----

# Plot pre-rendered graph
knitr::include_graphics("figures/1.2_deepSTRAPP_categorical_3lvl_data_4.3_plot_rates_through_time.PNG")


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

## ----plot_rates_vs_traits_cat_3lvl_eval_dev, fig.height = 7, fig.width = 8.5, out.width = "100%", eval = is_dev_version(), echo = FALSE----
# # 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)

## ----plot_rates_vs_traits_cat_3lvl_eval_CRAN, eval = !is_dev_version(), echo = FALSE, out.width = "100%"----

# Plot pre-rendered graph
knitr::include_graphics("figures/1.2_deepSTRAPP_categorical_3lvl_data_4.4_plot_rates_vs_traits.PNG")


## ----plot_updated_densityMaps_cat_3lvl----------------------------------------
# ### 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.
# 

## ----plot_updated_densityMaps_cat_3lvl_eval_dev, fig.height = 7, eval = is_dev_version(), echo = FALSE----
# # 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,
#                          cex_pies = 0.3,
#                          fsize = 0.1) # Reduce tip label size
# title(main = "Trait evolution for 100-0 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_40My$densityMaps,
#                          colors_per_levels = colors_per_states,
#                          cex_pies = 0.4,
#                          fsize = 0.2) # Reduce tip label size
# title(main = "Trait evolution for 100-40 My")

## ----plot_updated_densityMaps_cat_3lvl_eval_CRAN, eval = !is_dev_version(), echo = FALSE, out.width = "100%"----

# Plot pre-rendered graph
knitr::include_graphics("figures/1.2_deepSTRAPP_categorical_3lvl_data_4.5_plot_updated_densityMaps_1.PNG")
knitr::include_graphics("figures/1.2_deepSTRAPP_categorical_3lvl_data_4.5_plot_updated_densityMaps_2.PNG")


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

## ----plot_BAMM_rates_cat_3lvl_eval_dev, eval = is_dev_version(), echo = FALSE----
# old_par <- par(no.readonly = TRUE)
# par(mfrow = c(1, 2))
# 
# # 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, legend = TRUE, par.reset = FALSE)
# abline(v = max(phytools::nodeHeights(Ponerinae_tree_old_calib)[,2]) - 10, col = "red", lty = 2) # Show where the phylogeny will be cut at 10 Mya
# title(main = "BAMM rates for 100-0 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, legend = TRUE,
#                 colorbreaks = BAMM_map_10My$initial_colorbreaks$net_diversification)
# title(main = "BAMM rates for 100-10 My")
# 
# par(old_par)

## ----plot_BAMM_rates_cat_3lvl_eval_CRAN, eval = !is_dev_version(), echo = FALSE, out.width = "100%"----

# Plot pre-rendered graph
knitr::include_graphics("figures/1.2_deepSTRAPP_categorical_3lvl_data_4.6_plot_BAMM_rates.PNG")


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

## ----plot_traits_vs_rates_on_phylogeny_cat_3lvl_eval_dev, fig.height = 7, eval = is_dev_version(), echo = FALSE----
# 
# # 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)
# 
# # 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)
# 

## ----plot_traits_vs_rates_on_phylogeny_cat_3lvl_eval_CRAN, eval = !is_dev_version(), echo = FALSE, out.width = "100%"----

# Plot pre-rendered graph
knitr::include_graphics("figures/1.2_deepSTRAPP_categorical_3lvl_data_4.7_plot_traits_vs_rate_maps_1.PNG")
knitr::include_graphics("figures/1.2_deepSTRAPP_categorical_3lvl_data_4.7_plot_traits_vs_rate_maps_2.PNG")