Performance analysis of DAISIEprep::extract_island_species()
In this article we test and analyse the performance, in terms of time
consumption/complexity of the main function in the
DAISIEprep R package:
extract_island_species(). This function takes in a
phylogenetic tree and species endemicity data in the form of a
phylo4d object (an S4 class from the R package
phylobase).
This is not a thorough examination of all possible uses of the
extract_island_species() function, but rather gives an
indication of the usability of the function for data sets of different
sizes, as well as detecting any features of the data that may slow the
process of extracting and formatting the island community data.
All of the work for the performance analysis is carried out by the
benchmark() function from the DAISIEprep
package. But before calling this function we explain how the
benchmarking is set up.
The first argument of benchmark() is
phylod. If this is NULL then the function will
simulate the phylogenies and the endemicity data given the:
tree_size_range, num_points,
prob_on_island, and prob_endemic arguments.
These specify the lower and upper range of the tree size to simulate
(the sequence breaks can be in linear or log space depending on the
argument log_scale), the number of breaks in the sequence
between lower and upper tree sizes, the probability that a species will
be on the island, and if a species is on the island whether it is
endemic (1 - prob_endemic is the probability a species on
the island is non-endemic). If a phylo4d object is supplied
to the phylod argument in benchmark() then
this tree will be used to perform the benchmarking.
In the case of simulating data the parameter space of the performance
analysis is then the combination of these variables (using the
expand.grid() function). The phylogeny is simulated to a
fixed size using the rphylo() function from the
ape package.
As we are stochastically simulating the endemicity statuses of the species in the tree it may be that there is not an outgroup species not on the island in order to correctly extract the colonisation times from the stem age of the species or clade on the island. Therefore, we add an outgroup which we ensure is not present on the island.
When we want to test the performance of the
extract_island_species() function using the
asr algorithm we need to first run an ancestral state
reconstruction. This can easily be achieved if the phylod
data with the phylogeny and endemicity statuses at each tip of the tree
are ready (see tutorial vignette for more information on the
asr algorithm).
To quantify the time elapsed while the function runs there are
several methods in both base R and in various R packages
(e.g. microbenchmark or rbenchmark). However
here we use system.time() from base R (see https://radfordneal.wordpress.com/2014/02/02/inaccurate-results-from-microbenchmark/)
for reasoning).
Side note: by default DAISIEprep::benchmark() conditions
on each simulated data set having a non-empty island, and thus the
function is not tested for the trivial case that no species need to be
extracted.
Due to the computational time not being deterministic and to ensure
the results are not spurious we replicate timing three times and the
mean “real time elapsed” is calculated, as well as replicating the
simulation (given by argument replicates) to account for
stochasticity in the simulation of the data.
The range of tree sizes we use encompasses those of common empirical phylogenies (10 tips to 10,000 tips). We then generate a random sample of tip states for each species, with the possible states being: endemic to the island, non-endemic and not present on the island. We tested a low probability of each species in the tree being on the island (\(P(0.2)\)) and a high probability of species being on the island (\(P(0.5)\)). For each of these scenarios we tested a low (\(P_E(0.2)\), \(P_{NE}(0.8)\)) and high probability (\(P_E(0.8)\), \(P_{NE}(0.2)\)) of island species being endemic to test whether this affects performance. For each scenario we ran ten replicates. For isolated islands where cladogenesis is high, island species will likely not be spread across uniformly across the phylogeny as we assumed in our simulations, but instead will be clustered. To check for scenarios where the number of species per colonisation events is high (i.e. radiations on the island) we used the mammals of Madagascar data set. For this we used the (Upham et al. 2019) complete and DNA-only phylogenies, and the island endemicity data of Madagascar (Michielsen et al. 2023).
The results presented in this vignette are not computed each time the
vignette is rendered due to the large computation time required.
Instead, the analyses are run on a cluster computer and saved in the
package. The analysis script run to produce the results can be found in
the DAISIEprepExtra
package here.
The performance analysis uses the benchmark() function in
the DAISIEprep package.
The output produces results for the DNA-only phylogeny and the complete phylogeny. The raw data of parameter estimates for the different parameter settings is tidied into a tibble containing the data we need for both the DNA and complete phylogeny.
# Internal function - Not to be called for regular 'DAISIEprep' operation,
# but merely to load in data
performance_data <- DAISIEprep:::read_performance()We can plot the time consumption of each simulated and empirical
(DNA-only and complete phylogeny) data set and group the results by the
probability of species being on the island
(prob_on_island). Now the results can be plotted with
plot_performance implemented in DAISIEprep.
This function follows the tidyverse convention of giving variable names
as variables (as opposed to strings) and uses tidy evaluation to group
by the variable given, either prob_on_island or
prob_endemic. First for the min algorithm:
plot_performance(
performance_data = performance_data$performance_data_min,
group_by = prob_on_island
)
#> Warning: `position_dodge()` requires non-overlapping x intervals.
#> `position_dodge()` requires non-overlapping x intervals.
And secondly for the asr algorithm:
plot_performance(
performance_data = performance_data$performance_data_asr,
group_by = prob_on_island
)
#> Warning: `position_dodge()` requires non-overlapping x intervals.
#> `position_dodge()` requires non-overlapping x intervals.
Alternatively, the data can be grouped by the probability that
species on the island are endemic (prob_endemic) for the
min (first) and asr (second) algorithms:
plot_performance(
performance_data = performance_data$performance_data_min,
group_by = prob_endemic
)
#> Warning: `position_dodge()` requires non-overlapping x intervals.
#> `position_dodge()` requires non-overlapping x intervals.
plot_performance(
performance_data = performance_data$performance_data_min,
group_by = prob_endemic
)
#> Warning: `position_dodge()` requires non-overlapping x intervals.
#> `position_dodge()` requires non-overlapping x intervals.
We find that even for large trees (10,000 tips) the scale of extraction less than seconds, whereas running the DAISIE inference model is on the scale of minutes to days, and thus pre-processing does not present a bottleneck to the pipeline. The empirical trees we ran (the Madagascar mammal example) were quicker to process than trees simulated with uniformly random island presence (figure 1). This suggests that when the ratio of species to colonisation events on the island is higher it should be faster to extract. However, even for large trees with many colonisations extraction should pose a computational problem. The speed of extraction facilitates extraction of island data across many phylogenies, for example extracting island community data across the posterior distribution of inferred phylogenetic trees to account for uncertainty in branching times and tree topology.
Given the approximately linear relationship between the size of the phylogeny and time required for extraction in the above log-log plots we fit a power law to the mean run time (mean across the replicates).
grouped_performance_data <- dplyr::group_by(
performance_data$performance_data_min,
tree_size,
"prob_on_island"
)
mean_performance_data <- dplyr::summarise(
grouped_performance_data,
mean = mean(median_time),
.groups = "drop"
)
fit_min <- lm(log(mean_performance_data$mean) ~ log(mean_performance_data$tree_size))
fit_min$coefficients
#> (Intercept) log(mean_performance_data$tree_size)
#> -28.297771 1.399895
grouped_performance_data <- dplyr::group_by(
performance_data$performance_data_asr,
tree_size,
"prob_on_island"
)
mean_performance_data <- dplyr::summarise(
grouped_performance_data,
mean = mean(median_time),
.groups = "drop"
)
fit_asr <- lm(log(mean_performance_data$mean) ~ log(mean_performance_data$tree_size))
fit_asr$coefficients
#> (Intercept) log(mean_performance_data$tree_size)
#> -28.0177 1.3008The growth in time for both the min and asr
algorithms follows a power function (\(y=ax^k\)) with exponents \(k =\) 1.3998949 and \(k =\) 1.3008005, respectively, when fitted
on mean run time against tree size. Therefore the time complexity scales
relatively well with tree size, and thus unless needing to be applied to
extremely large phylogenies (>10,000 tips) the
DAISIEprep extraction functionality should be
applicable.