In TGuillerme/dispRity: Measuring Disparity

Details of specific functions

The following section contains information specific to some functions. If any of your questions are not covered in these sections, please refer to the function help files in R, send me an email (guillert@tcd.ie), or raise an issue on GitHub. The several tutorials below describe specific functionalities of certain functions; please always refer to the function help files for the full function documentation!

Before each section, make sure you loaded the @beckancient2014 data (see example data for more details).

## Loading the data
data(BeckLee_mat50)
data(BeckLee_mat99)
data(BeckLee_tree)
data(BeckLee_ages)

Time slicing {#chrono-subsets}

The function chrono.subsets allows users to divide the matrix into different time subsets or slices given a dated phylogeny that contains all the elements (i.e. taxa) from the matrix. Each subset generated by this function will then contain all the elements present at a specific point in time or during a specific period in time.

Two types of time subsets can be performed by using the method option:

Discrete time subsets (or time-binning) using method = discrete
Continuous time subsets (or time-slicing) using method = continuous

For the time-slicing method details see @time-slice. For both methods, the function takes the time argument which can be a vector of numeric values for:

Defining the boundaries of the time bins (when method = discrete)
Defining the time slices (when method = continuous)

Otherwise, the time argument can be set as a single numeric value for automatically generating a given number of equidistant time-bins/slices. Additionally, it is also possible to input a dataframe containing the first and last occurrence data (FAD/LAD) for taxa that span over a longer time than the given tips/nodes age, so taxa can appear in more than one time bin/slice.

Time-binning

Here is an example for the time binning method (method = discrete):

## Generating three time bins containing the taxa present every 40 Ma
chrono.subsets(data = BeckLee_mat50, tree = BeckLee_tree,
               method = "discrete",
               time = c(120, 80, 40, 0))

Note that we can also generate equivalent results by just telling the function that we want three time-bins as follow:

## Automatically generate three equal length bins:
chrono.subsets(data = BeckLee_mat50, tree = BeckLee_tree,
               method = "discrete",
               time = 3)

In this example, the taxa were split inside each time-bin according to their age. However, the taxa here are considered as single points in time. It is totally possible that some taxa could have had longer longevity and that they exist in multiple time bins. In this case, it is possible to include them in more than one bin by providing a table of first and last occurrence dates (FAD/LAD). This table should have the taxa names as row names and two columns for respectively the first and last occurrence age:

## Displaying the table of first and last occurrence dates
## for each taxa
head(BeckLee_ages)

## Generating time bins including taxa that might span between them
chrono.subsets(data = BeckLee_mat50, tree = BeckLee_tree,
               method = "discrete",
               time = c(120, 80, 40, 0), FADLAD = BeckLee_ages)

When using this method, the oldest boundary of the first bin (or the first slice, see below) is automatically generated as the root age plus 1\% of the tree length, as long as at least three elements/taxa are present at that point in time. The algorithm adds an extra 1\% tree length until reaching the required minimum of three elements. It is also possible to include nodes in each bin by using inc.nodes = TRUE and providing a matrix that contains the ordinated distance among tips and nodes.

If you want to generate time subsets based on stratigraphy, the package proposes a useful functions to do it for you: get.bin.ages (check out the function's manual in R)!

Time-slicing {#time-slicing}

For the time-slicing method (method = continuous), the idea is fairly similar. This option, however, requires a matrix that contains the ordinated distance among taxa and nodes and an extra argument describing the assumed evolutionary model (via the model argument). This model argument is used when the time slice occurs along a branch of the tree rather than on a tip or a node, meaning that a decision must be made about what the value for the branch should be. The model can be one of the following:

Punctuated models
acctran where the data chosen along the branch is always the one of the descendant
deltran where the data chosen along the branch is always the one of the ancestor
random where the data chosen along the branch is randomly chosen between the descendant or the ancestor
proximity where the data chosen along the branch is either the descendant or the ancestor depending on branch length
Gradual models
equal.split where the data chosen along the branch is both the descendant and the ancestor with an even probability
gradual.split where the data chosen along the branch is both the descendant and the ancestor with a probability depending on branch length

Note that the four first models are a proxy for punctuated evolution: the selected data is always either the one of the descendant or the ancestor. In other words, changes along the branches always occur at either ends of it. The two last models are a proxy for gradual evolution: the data from both the descendant and the ancestor is used with an associate probability. These later models perform better when bootstrapped, effectively approximating the "intermediate" state between and the ancestor and the descendants.

More details about the differences between these methods can be found in @time-slice.

## Generating four time slices every 40 million years
## under a model of proximity evolution
chrono.subsets(data = BeckLee_mat99, tree = BeckLee_tree, 
               method = "continuous", model = "proximity",
               time = c(120, 80, 40, 0),
               FADLAD = BeckLee_ages)

## Generating four time slices automatically
chrono.subsets(data = BeckLee_mat99, tree = BeckLee_tree,
               method = "continuous", model = "proximity",
               time = 4, FADLAD = BeckLee_ages)

Customised subsets {#custom-subsets}

Another way of separating elements into different categories is to use customised subsets as briefly explained above. This function simply takes the list of elements to put in each group (whether they are the actual element names or their position in the matrix).

## Creating the two groups (crown and stems)
mammal_groups <- crown.stem(BeckLee_tree, inc.nodes = FALSE)

## Separating the dataset into two different groups
custom.subsets(BeckLee_mat50, group = mammal_groups)

Like in this example, you can use the utility function crown.stem that allows to automatically separate the crown and stems taxa given a phylogenetic tree. Also, elements can easily be assigned to different groups if necessary!

## Creating the three groups as a list
weird_groups <- list("even" = seq(from = 1, to = 49, by = 2),
                      "odd" = seq(from = 2, to = 50, by = 2),
                      "all" = c(1:50))

The custom.subsets function can also take a phylogeny (as a phylo object) as an argument to create groups as clades:

## Creating groups as clades
custom.subsets(BeckLee_mat50, group = BeckLee_tree)

This automatically creates 49 (the number of nodes) groups containing between two and 50 (the number of tips) elements.

Bootstraps and rarefactions

One important step in analysing ordinated matrices is to pseudo-replicate the data to see how robust the results are, and how sensitive they are to outliers in the dataset. This can be achieved using the function boot.matrix to bootstrap and/or rarefy the data. The default options will bootstrap the matrix 100 times without rarefaction using the "full" bootstrap method (see below):

## Default bootstrapping
boot.matrix(data = BeckLee_mat50)

The number of bootstrap replicates can be defined using the bootstraps option. The method can be modified by controlling which bootstrap algorithm to use through the boot.type argument. Currently two algorithms are implemented:

"full" where the bootstrapping is entirely stochastic (n elements are replaced by any m elements drawn from the data)
"single" where only one random element is replaced by one other random element for each pseudo-replicate
"null" where every element is resampled across the whole matrix (not just the subsets). I.e. for each subset of n elements, this algorithm resamples n elements across ALL subsets (not just the current one). If only one subset (or none) is used, this does the same as the "full" algorithm.

## Bootstrapping with the single bootstrap method
boot.matrix(BeckLee_mat50, boot.type = "single")

This function also allows users to rarefy the data using the rarefaction argument. Rarefaction allows users to limit the number of elements to be drawn at each bootstrap replication. This is useful if, for example, one is interested in looking at the effect of reducing the number of elements on the results of an analysis.

This can be achieved by using the rarefaction option that draws only n-x at each bootstrap replicate (where x is the number of elements not sampled). The default argument is FALSE but it can be set to TRUE to fully rarefy the data (i.e. remove x elements for the number of pseudo-replicates, where x varies from the maximum number of elements present in each subset to a minimum of three elements). It can also be set to one or more numeric values to only rarefy to the corresponding number of elements.

## Bootstrapping with the full rarefaction
boot.matrix(BeckLee_mat50, bootstraps = 20,
            rarefaction = TRUE)

## Or with a set number of rarefaction levels
boot.matrix(BeckLee_mat50, bootstraps = 20,
            rarefaction = c(6:8, 3))

Note that using the rarefaction argument also bootstraps the data. In these examples, the function bootstraps the data (without rarefaction) AND also bootstraps the data with the different rarefaction levels.

## Creating subsets of crown and stem mammals
crown_stem <- custom.subsets(BeckLee_mat50,
                             group = crown.stem(BeckLee_tree,
                                                inc.nodes = FALSE))
## Bootstrapping and rarefying these groups
boot.matrix(crown_stem, bootstraps = 200, rarefaction = TRUE)

## Creating time slice subsets
time_slices <- chrono.subsets(data = BeckLee_mat99,
                              tree = BeckLee_tree, 
                              method = "continuous",
                              model = "proximity", 
                              time = c(120, 80, 40, 0),
                              FADLAD = BeckLee_ages)

## Bootstrapping the time slice subsets
boot.matrix(time_slices, bootstraps = 100)

Bootstrapping with probabilities

It is also possible to specify the sampling probability in the bootstrap for each elements. This can be useful for weighting analysis for example (i.e. giving more importance to specific elements). These probabilities can be passed to the prob argument individually with a vector with the elements names or with a matrix with the rownames as elements names. The elements with no specified probability will be assigned a probability of 1 (or 1/maximum weight if the argument is weights rather than probabilities).

## Attributing a weight of 0 to Cimolestes and 10 to Maelestes
boot.matrix(BeckLee_mat50,
            prob = c("Cimolestes" = 0, "Maelestes" = 10))

Bootstrapping dimensions

In some cases, you might also be interested in bootstrapping dimensions rather than observations. I.e. bootstrapping the columns of a matrix rather than the rows.

It's pretty easy! By default, boot.matrix uses the option boot.by = "rows" which you can toggle to boot.by = "columns"

## Bootstrapping the observations (default)
set.seed(1)
boot_obs <- boot.matrix(data = crown_stem, boot.by = "rows")

## Bootstrapping the columns rather than the rows
set.seed(1)
boot_dim <- boot.matrix(data = crown_stem, boot.by = "columns")

In these two examples, the first one boot_obs bootstraps the rows as showed before (default behaviour). But the second one, boot_dim bootstraps the dimensions. That means that for each bootstrap sample, the value calculated is actually obtained by reshuffling the dimensions (columns) rather than the observations (rows).

## Measuring disparity and summarising
summary(dispRity(boot_obs, metric = sum))
summary(dispRity(boot_dim, metric = sum))

Note here how the observed sum is the same (no bootstrapping) but the bootstrapping distributions are quiet different even though the same seed was used.

Disparity metrics {#disparity-metrics}

There are many ways of measuring disparity! In brief, disparity is a summary metric that will represent an aspect of an ordinated space (e.g. a MDS, PCA, PCO, PCoA). For example, one can look at ellipsoid hyper-volume of the ordinated space (Donohue et al. 2013), the sum and the product of the ranges and variances (Wills et al. 1994) or the median position of the elements relative to their centroid (Wills et al. 1994). Of course, there are many more examples of metrics one can use for describing some aspect of the ordinated space, with some performing better than other ones at particular descriptive tasks, and some being more generalist. Check out this paper on selecting the best metric for your specific question in Ecology and Evolution. You can also use the moms shiny app to test which metric captures which aspect of traitspace occupancy regarding your specific space and your specific question.

Regardless, and because of this great diversity of metrics, the package dispRity does not have one way to measure disparity but rather proposes to facilitate users in defining their own disparity metric that will best suit their particular analysis. In fact, the core function of the package, dispRity, allows the user to define any metric with the metric argument. However the metric argument has to follow certain rules:

It must be composed from one to three function objects;
The function(s) must take as a first argument a matrix or a vector;
The function(s) must be of one of the three dimension-levels described below;
At least one of the functions must be of dimension-level 1 or 2 (see below).

The function dimension-levels

The metric function dimension-levels determine the "dimensionality of decomposition" of the input matrix. In other words, each dimension-level designates the dimensions of the output, i.e. either three (a matrix); two (a vector); or one (a single numeric value) dimension.

Illustration of the different dimension-levels of functions with an input matrix

Dimension-level 1 functions

A dimension-level 1 function will decompose a matrix or a vector into a single value:

## Creating a dummy matrix
dummy_matrix <- matrix(rnorm(12), 4, 3)

## Example of dimension-level 1 functions
mean(dummy_matrix)
median(dummy_matrix)

Any summary metric such as mean or median are good examples of dimension-level 1 functions as they reduce the matrix to a single dimension (i.e. one value).

Dimension-level 2 functions

A dimension-level 2 function will decompose a matrix into a vector.

## Defining the function as the product of rows
prod.rows <- function(matrix) apply(matrix, 1, prod)

## A dimension-level 2 metric
prod.rows(dummy_matrix)

Several dimension-level 2 functions are implemented in dispRity (see ?dispRity.metric) such as the variances or ranges functions that calculate the variance or the range of each dimension of the ordinated matrix respectively.

Dimension-level 3 functions

Finally a dimension-level 3 function will transform the matrix into another matrix. Note that the dimension of the output matrix doesn't need to match the the input matrix:

## A dimension-level 3 metric
var(dummy_matrix)

## A dimension-level 3 metric with a forced matrix output
as.matrix(dist(dummy_matrix))

Between groups metrics {#betweengroupmetricsexplain}

One specific category of metrics in the dispRity package is the between groups metrics. As the name suggest, these metrics can be used to calculate the disparity between groups rather than within the groups. These metrics follow the same classifications as the "normal" (within group) metrics with dimension-level 1, 2 and 3 between groups metrics. However, at the difference of the "normal" metrics, their input arguments must be matrix and matrix2 (and of course any other additional arguments). For example, this metric measures the difference in mean between two matrices:

## A simple example
mean.difference <- function(matrix, matrix2) {
  mean(matrix) - mean(matrix2)
}

You can find the list of implemented between groups metric here or design them yourself for your specific needs (potentially using make.metric for help).

The function works by simply using the two available matrices, with no restriction in terms of dimensions (although you'd probably want both matrices to have the same number of dimensions)

## A second matrix
dummy_matrix2 <- matrix(runif(12), 4, 3)

## The difference between groups
mean.difference(dummy_matrix, dummy_matrix2)

Beyond this super simple example, it might probably be interesting to use this metric on dispRity objects, especially the ones from custom.subsets and chrono.subsets. In fact, the dispRity function allows to apply the between groups metric directly to the dispRity objects using the between.groups = TRUE option. For example:

## Combining both matrices
big_matrix <- rbind(dummy_matrix, dummy_matrix2)
rownames(big_matrix) <- 1:8

## Making a dispRity object with both groups
grouped_matrix <- custom.subsets(big_matrix,
                                 group = c(list(1:4), list(1:4)))

## Calculating the mean difference between groups
(mean_differences <- dispRity(grouped_matrix,
                              metric = mean.difference,
                              between.groups = TRUE))

## Summarising the object
summary(mean_differences)
## Note how the summary table now indicates
## the number of elements for each group

For dispRity objects generated by custom.subsets, the dispRity function will by default apply the metric on the groups in a pairwise fashion. For example, if the object contains multiple groups, all groups will be compared to each other:

## A dispRity object with multiple groups
grouped_matrix <- custom.subsets(big_matrix,
                                 group = c("A" = list(1:4),
                                           "B" = list(1:4),
                                           "C" = list(2:6), 
                                           "D" = list(1:8)))

## Measuring disparity between all groups
summary(dispRity(grouped_matrix, metric = mean.difference,
                 between.groups = TRUE))

For dispRity objects generated by chrono.subsets (not shown here), the dispRity function will by default apply the metric on the groups in a serial way (group 1 vs. group 2, group 2 vs. group 3, group 3 vs. group 4, etc...). However, in both cases (for objects from custom.subsets or chrono.subsets) it is possible to manually specific the list of pairs of comparisons through their ID numbers:

## Measuring disparity between specific groups
summary(dispRity(grouped_matrix, metric = mean.difference,
                 between.groups = list(c(1,3), c(3,1), c(4,1))))

Note that in any case, the order of the comparison can matter. In our example, it is obvious that mean(matrix) - mean(matrix2) is not the same as mean(matrix2) - mean(matrix).

`make.metric` {#makemetric}

Of course, functions can be more complex and involve multiple operations such as the centroids function (see ?dispRity.metric) that calculates the Euclidean distance between each element and the centroid of the ordinated space. The make.metric function implemented in dispRity is designed to help test and find the dimension-level of the functions. This function tests:

If your function can deal with a matrix or a vector as an input;
Your function's dimension-level according to its output (dimension-level 1, 2 or 3, see above);
Whether the function can be implemented in the dispRity function (the function is fed into a lapply loop).

For example, let's see if the functions described above are the right dimension-levels:

## Which dimension-level is the mean function?
## And can it be used in dispRity?
make.metric(mean)

## Which dimension-level is the prod.rows function?
## And can it be used in dispRity?
make.metric(prod.rows)

## Which dimension-level is the var function?
## And can it be used in dispRity?
make.metric(var)

A non verbose version of the function is also available. This can be done using the option silent = TRUE and will simply output the dimension-level of the metric.

## Testing whether mean is dimension-level 1
if(make.metric(mean, silent = TRUE)$type != "level1") {
    message("The metric is not dimension-level 1.")
}
## Testing whether var is dimension-level 1
if(make.metric(var, silent = TRUE)$type != "level1") {
    message("The metric is not dimension-level 1.")
}

Metrics in the `dispRity` function

Using this metric structure, we can easily use any disparity metric in the dispRity function as follows:

## Measuring disparity as the standard deviation
## of all the values of the
## ordinated matrix (dimension-level 1 function).
summary(dispRity(BeckLee_mat50, metric = sd))

## Measuring disparity as the standard deviation
## of the variance of each axis of
## the ordinated matrix (dimension-level 1 and 2 functions).
summary(dispRity(BeckLee_mat50, metric = c(sd, variances)))

## Measuring disparity as the standard deviation
## of the variance of each axis of
## the variance covariance matrix (dimension-level 1, 2 and 3 functions).
summary(dispRity(BeckLee_mat50, metric = c(sd, variances, var)), round = 10)

Note that the order of each function in the metric argument does not matter, the dispRity function will automatically detect the function dimension-levels (using make.metric) and apply them to the data in decreasing order (dimension-level 3 > 2 > 1).

## Disparity as the standard deviation of the variance of each axis of the
## variance covariance matrix:
disparity1 <- summary(dispRity(BeckLee_mat50,
                               metric = c(sd, variances, var)),
                      round = 10)

## Same as above but using a different function order for the metric argument
disparity2 <- summary(dispRity(BeckLee_mat50,
                               metric = c(variances, sd, var)),
                      round = 10)

## Both ways output the same disparity values:
disparity1 == disparity2

In these examples, we considered disparity to be a single value. For example, in the previous example, we defined disparity as the standard deviation of the variances of each column of the variance/covariance matrix (metric = c(variances, sd, var)). It is, however, possible to calculate disparity as a distribution.

Metrics implemented in `dispRity`

Several disparity metrics are implemented in the dispRity package. The detailed list can be found in ?dispRity.metric along with some description of each metric.

Level | Name | Description | Source | ------|------|---------------------------------------------------|--------| 2 | ancestral.dist | The distance between an element and its ancestor | dispRity | 2 | angles | The angle of main variation of each dimensions | dispRity | 2 | centroids¹ | The distance between each element and the centroid of the ordinated space | dispRity | 1 | convhull.surface | The surface of the convex hull formed by all the elements | geometry::convhulln$area | 1 | convhull.volume | The volume of the convex hull formed by all the elements | geometry::convhulln$vol | 2 | count.neighbours | The number of neigbhours to each element in a specified radius | dispRity | 2 | deviations | The minimal distance between each element and a hyperplane | dispRity | 1 | diagonal | The longest distance in the ordinated space (like the diagonal in two dimensions) | dispRity | 1 | disalignment | The rejection of the centroid of a matrix from the major axis of another (typically an "as.covar" metric) | dispRity | 2 | displacements | The ratio between the distance from a reference and the distance from the centroid | dispRity | 1 | edge.length.tree | The edge lengths of the elements on a tree | ape | 1 | ellipsoid.volume¹ | The volume of the ellipsoid of the space | Donohue et al. (2013) | 1 | func.div | The functional divergence (the ratio of deviation from the centroid) | dispRity (similar to FD::dbFD$FDiv but without abundance)| 1 | func.eve | The functional evenness (the minimal spanning tree distances evenness) | dispRity (similar to FD::dbFD$FEve but without abundance)| 1 | group.dist | The distance between two groups | dispRity | 1 | mode.val | The modal value | dispRity | 1 | n.ball.volume | The hyper-spherical (n-ball) volume | dispRity | 2 | neighbours | The distance to specific neighbours (e.g. the nearest neighbours - by default) | dispRity | 2 | pairwise.dist | The pairwise distances between elements | vegan::vegist | 2 | point.dist | The distance between one group and the point of another group | dispRity | 2 | projections | The distance on (projection) or from (rejection) an arbitrary vector | dispRity | 1 | projections.between | projections metric applied between groups | dispRity | 2 | projections.tree | The projections metric but where the vector can be based on a tree | dispRity | 2 | quantiles | The nth quantile range per axis | dispRity | 2 | radius | The radius of each dimensions | dispRity | 2 | ranges | The range of each dimension | dispRity | 1 | roundness | The integral of the ranked scaled eigenvalues of a variance-covariance matrix | dispRity 2 | span.tree.length | The minimal spanning tree length | vegan::spantree | 2 | variances | The variance of each dimension | dispRity |

1: Note that by default, the centroid is the centroid of the elements. It can, however, be fixed to a different value by using the centroid argument centroids(space, centroid = rep(0, ncol(space))), for example the origin of the ordinated space.

2: This function uses an estimation of the eigenvalue that only works for MDS or PCoA ordinations (not PCA).

You can find more informations on the vast variety of metrics that you can use in your analysis in this paper.

Equations and implementations

Some of the functions described below are implemented in the dispRity package and do not require any other packages to calculate (see implementation here).

\begin{equation} ancestral.dist = \sqrt{\sum_{i=1}^{n}{({d}{n}-Ancestor{n})^2}} \end{equation}

\begin{equation} centroids = \sqrt{\sum_{i=1}^{n}{({d}{n}-Centroid{d})^2}} \end{equation}

\begin{equation} diagonal = \sqrt{\sum_{i=1}^{d}|max(d_i) - min(k_i)|} \end{equation}

\begin{equation} deviations = \frac{|Ax + By + ... + Nm + Intercept|}{\sqrt{A^2 + B^2 + ... + N^2}} \end{equation}

\begin{equation} displacements = \frac{\sqrt{\sum_{i=1}^{n}{({d}{n}-Reference{d})^2}}}{\sqrt{\sum_{i=1}^{n}{({d}{n}-Centroid{k})^2}}} \end{equation}

\begin{equation} ellipsoid.volume = \frac{\pi^{d/2}}{\Gamma(\frac{d}{2}+1)}\displaystyle\prod_{i=1}^{d} (\lambda_{i}^{0.5}) \end{equation}

\begin{equation} n.ball.volume = \frac{\pi^{d/2}}{\Gamma(\frac{d}{2}+1)}\displaystyle\prod_{i=1}^{d} R \end{equation}

\begin{equation} projection_{on} = \| \overrightarrow{i} \cdot \overrightarrow{b} \| \end{equation} \begin{equation} projection_{from} = \| \overrightarrow{i} - \overrightarrow{i} \cdot \overrightarrow{b} \| \end{equation}

\begin{equation} radius = |\frac{\sum_{i=1}^{n}d_i}{n} - f(\mathbf{v}d)| \end{equation}

\begin{equation} ranges = |max(d_i) - min(d_i)| \end{equation}

\begin{equation} roundness = \int_{i = 1}^{n}{\frac{\lambda_{i}}{\text{max}(\lambda)}} \end{equation}

\begin{equation} variances = \sigma^{2}{d_i} \end{equation}

\begin{equation} span.tree.length = \mathrm{branch\ length} \end{equation}

Where d is the number of dimensions, n the number of elements, $\Gamma$ is the Gamma distribution, $\lambda_i$ is the eigenvalue of each dimensions, $\sigma^{2}$ is their variance and $Centroid_{k}$ is their mean, $Ancestor_{n}$ is the coordinates of the ancestor of element $n$, $f(\mathbf{v}k)$ is function to select one value from the vector $\mathbf{v}$ of the dimension $k$ (e.g. it's maximum, minimum, mean, etc.), R is the radius of the sphere or the product of the radii of each dimensions ($\displaystyle\prod_{i=1}^{k}R_{i}$ - for a hyper-ellipsoid), $Reference_{k}$ is an arbitrary point's coordinates (usually 0), $\overrightarrow{b}$ is the vector defined by ((point1, point2)), and $\overrightarrow{i}$ is the vector defined by ((point1, i) where i is any row of the matrix).

Using the different disparity metrics

Here is a brief demonstration of the main metrics implemented in dispRity. First, we will create a dummy/simulated ordinated space using the space.maker utility function (more about that here:

## Creating a 10*5 normal space
set.seed(1)
dummy_space <- space.maker(10, 5, rnorm)
rownames(dummy_space) <- 1:10

We will use this simulated space to demonstrate the different metrics.

Volumes and surface metrics

The functions ellipsoid.volume, convhull.surface, convhull.volume and n.ball.volume all measure the surface or the volume of the ordinated space occupied:

Because there is only one subset (i.e. one matrix) in the dispRity object, the operations below are the equivalent of metric(dummy_space) (with rounding).

## Calculating the ellipsoid volume
summary(dispRity(dummy_space, metric = ellipsoid.volume))

WARNING: in such dummy space, this gives the estimation of the ellipsoid volume, not the real ellipsoid volume! See the cautionary note in ?ellipsoid.volume.

## Calculating the convex hull surface
summary(dispRity(dummy_space, metric = convhull.surface))

## Calculating the convex hull volume
summary(dispRity(dummy_space, metric = convhull.volume))

## Calculating the convex hull volume
summary(dispRity(dummy_space, metric = n.ball.volume))

The convex hull based functions are a call to the geometry::convhulln function with the "FA" option (computes total area and volume). Also note that they are really sensitive to the size of the dataset.

Cautionary note: measuring volumes in a high number of dimensions can be strongly affected by the curse of dimensionality that often results in near 0 disparity values. I strongly recommend reading this really intuitive explanation from Toph Tucker.

Ranges, variances, quantiles, radius, pairwise distance, neighbours (and counting them), modal value and diagonal

The functions ranges, variances radius, pairwise.dist, mode.val and diagonal all measure properties of the ordinated space based on its dimensional properties (they are also less affected by the "curse of dimensionality"):

ranges, variances quantiles and radius work on the same principle and measure the range/variance/radius of each dimension:

## Calculating the ranges of each dimension in the ordinated space
ranges(dummy_space)

## Calculating disparity as the distribution of these ranges
summary(dispRity(dummy_space, metric = ranges))

## Calculating disparity as the sum and the product of these ranges
summary(dispRity(dummy_space, metric = c(sum, ranges)))
summary(dispRity(dummy_space, metric = c(prod, ranges)))

## Calculating the variances of each dimension in the
## ordinated space
variances(dummy_space)

## Calculating disparity as the distribution of these variances
summary(dispRity(dummy_space, metric = variances))

## Calculating disparity as the sum and
## the product of these variances
summary(dispRity(dummy_space, metric = c(sum, variances)))
summary(dispRity(dummy_space, metric = c(prod, variances)))

## Calculating the quantiles of each dimension
## in the ordinated space
quantiles(dummy_space)

## Calculating disparity as the distribution of these variances
summary(dispRity(dummy_space, metric = quantiles))

## By default, the quantile calculated is the 95%
## (i.e. 95% of the data on each axis)
## this can be changed using the option quantile:
summary(dispRity(dummy_space, metric = quantiles, quantile = 50))

## Calculating the radius of each dimension in the ordinated space
radius(dummy_space)

## By default the radius is the maximum distance from the centre of
## the dimension. It can however be changed to any function:
radius(dummy_space, type = min)
radius(dummy_space, type = mean)

## Calculating disparity as the mean average radius
summary(dispRity(dummy_space,
                 metric = c(mean, radius),
                 type = mean))

The pairwise distances and the neighbours distances uses the function vegan::vegdist and can take the normal vegdist options:

## The average pairwise euclidean distance
summary(dispRity(dummy_space, metric = c(mean, pairwise.dist)))

## The distribution of the Manhattan distances
summary(dispRity(dummy_space, metric = pairwise.dist,
                 method = "manhattan"))

## The average nearest neighbour distances
summary(dispRity(dummy_space, metric = neighbours))

## The average furthest neighbour manhattan distances
summary(dispRity(dummy_space, metric = neighbours,
                 which = max, method = "manhattan"))

## The overall number of neighbours per point
summary(dispRity(dummy_space, metric = count.neighbours,
                 relative = FALSE))

## The relative number of neigbhours
## two standard deviations of each element
summary(dispRity(dummy_space, metric = count.neighbours,
                 radius = function(x)(sd(x)*2),
                 relative = TRUE))

Note that this function is a direct call to vegan::vegdist(matrix, method = method, diag = FALSE, upper = FALSE, ...).

The diagonal function measures the multidimensional diagonal of the whole space (i.e. in our case the longest Euclidean distance in our five dimensional space). The mode.val function measures the modal value of the matrix:

## Calculating the ordinated space's diagonal
summary(dispRity(dummy_space, metric = diagonal))

## Calculating the modal value of the matrix
summary(dispRity(dummy_space, metric = mode.val))

This metric is only a Euclidean diagonal (mathematically valid) if the dimensions within the space are all orthogonal!

Centroids, displacements and ancestral distances metrics {#centroids}

The centroids metric allows users to measure the position of the different elements compared to a fixed point in the ordinated space. By default, this function measures the distance between each element and their centroid (centre point):

## The distribution of the distances between each element and their centroid
summary(dispRity(dummy_space, metric = centroids))

## Disparity as the median value of these distances
summary(dispRity(dummy_space, metric = c(median, centroids)))

It is however possible to fix the coordinates of the centroid to a specific point in the ordinated space, as long as it has the correct number of dimensions:

## The distance between each element and the origin
## of the ordinated space
summary(dispRity(dummy_space, metric = centroids, centroid = 0))

## Disparity as the distance between each element
## and a specific point in space
summary(dispRity(dummy_space, metric = centroids,
                 centroid = c(0,1,2,3,4)))

If you have subsets in your dispRity object, you can also use the matrix.dispRity (see utilities) and colMeans to get the centre of a specific subgroup. For example

## Create a custom subsets object
dummy_groups <- custom.subsets(dummy_space,
                               group = list("group1" = 1:5,
                                            "group2" = 6:10))
summary(dispRity(dummy_groups, metric = centroids,
    centroid = colMeans(get.matrix(dummy_groups, "group1"))))

The displacements distance is the ratio between the centroids distance and the centroids distance with centroid = 0. Note that it is possible to measure a ratio from another point than 0 using the reference argument. It gives indication of the relative displacement of elements in the multidimensional space: a score >1 signifies a displacement away from the reference. A score of >1 signifies a displacement towards the reference.

## The relative displacement of the group in space to the centre
summary(dispRity(dummy_space, metric = displacements))

## The relative displacement of the group to an arbitrary point
summary(dispRity(dummy_space, metric = displacements,
                 reference = c(0,1,2,3,4)))

The ancestral.dist metric works on a similar principle as the centroids function but changes the centroid to be the coordinates of each element's ancestor (if to.root = FALSE; default) or to the root of the tree (to.root = TRUE). Therefore this function needs a matrix that contains tips and nodes and a tree as additional argument.

## A generating a random tree with node labels
my_tree <- makeNodeLabel(rtree(5), prefix = "n")
## Adding the tip and node names to the matrix
dummy_space2 <- dummy_space[-1,]
rownames(dummy_space2) <- c(my_tree$tip.label,
                            my_tree$node.label)

## Calculating the distances from the ancestral nodes
ancestral_dist <- dispRity(dummy_space2, metric = ancestral.dist,
                           tree = my_tree)

## The ancestral distances distributions
summary(ancestral_dist)

## Calculating disparity as the sum of the distances from all the ancestral nodes
summary(dispRity(ancestral_dist, metric = sum))

Minimal spanning tree length

The span.tree.length uses the vegan::spantree function to heuristically calculate the minimum spanning tree (the shortest multidimensional tree connecting each elements) and calculates its length as the sum of every branch lengths.

## The length of the minimal spanning tree
summary(dispRity(dummy_space, metric = c(sum, span.tree.length)))

Note that because the solution is heuristic, this metric can take a long time to compute for big matrices.

Functional divergence and evenness

The func.div and func.eve functions are based on the FD::dpFD package. They are the equivalent to FD::dpFD(matrix)$FDiv and FD::dpFD(matrix)$FEve but a bit faster (since they don't deal with abundance data). They are pretty straightforward to use:

## The ratio of deviation from the centroid 
summary(dispRity(dummy_space, metric = func.div))

## The minimal spanning tree distances evenness
summary(dispRity(dummy_space, metric = func.eve))

## The minimal spanning tree manhanttan distances evenness
summary(dispRity(dummy_space, metric = func.eve,
                 method = "manhattan"))

Orientation: angles and deviations

The angles performs a least square regression (via the lm function) and returns slope of the main axis of variation for each dimension. This slope can be converted into different units, "slope", "degree" (the default) and "radian". This can be changed through the unit argument. By default, the angle is measured from the slope 0 (the horizontal line in a 2D plot) but this can be changed through the base argument (using the defined unit):

## The distribution of each angles in degrees for each
## main axis in the matrix
summary(dispRity(dummy_space, metric = angles))

## The distribution of slopes deviating from the 1:1 slope:
summary(dispRity(dummy_space, metric = angles, unit = "slope",
                 base = 1))

The deviations function is based on a similar algorithm as above but measures the deviation from the main axis (or hyperplane) of variation. In other words, it finds the least square line (for a 2D dataset), plane (for a 3D dataset) or hyperplane (for a >3D dataset) and measures the shortest distances between every points and the line/plane/hyperplane. By default, the hyperplane is fitted using the least square algorithm from stats::glm:

## The distribution of the deviation of each point
## from the least square hyperplane
summary(dispRity(dummy_space, metric = deviations))

It is also possible to specify the hyperplane equation through the hyperplane equation. The equation must contain the intercept first and then all the slopes and is interpreted as $intercept + Ax + By + ... + Nd = 0$. For example, a 2 line defined as beta + intercept (e.g. $y = 2x + 1$) should be defined as hyperplane = c(1, 2, 1) ($2x - y + 1 = 0$).

## The distribution of the deviation of each point
## from a slope (with only the two first dimensions)
summary(dispRity(dummy_space[, c(1:2)], metric = deviations,
                 hyperplane = c(1, 2, -1)))

Since both the functions angles and deviations effectively run a lm or glm to estimate slopes or hyperplanes, it is possible to use the option significant = TRUE to only consider slopes or intercepts that have a slope significantly different than zero using an aov with a significant threshold of $p = 0.05$. Note that depending on your dataset, using and aov could be completely inappropriate! In doubt, it's probably better to enter your base (for angles) or your hyperplane (for deviations) manually so you're sure you know what the function is measuring.

Projections and phylo projections: elaboration and exploration {#projections}

The projections metric calculates the geometric projection and corresponding rejection of all the rows in a matrix on an arbitrary vector (respectively the distance on and the distance from that vector). The function is based on @aguilera2004's n-dimensional rotation algorithm to use linear algebra in mutidimensional spaces. The projection or rejection can be seen as respectively the elaboration and exploration scores on a trajectory (sensu @endler2005).

By default, the vector (e.g. a trajectory, an axis), on which the data is projected is the one going from the centre of the space (coordinates 0,0, ...) and the centroid of the matrix. However, we advice you do define this axis to something more meaningful using the point1 and point2 options, to create the vector (the vector's norm will be dist(point1, point2) and its direction will be from point1 towards point2).

## The elaboration on the axis defined by the first and
## second row in the dummy_space
summary(dispRity(dummy_space, metric = projections,
                              point1 = dummy_space[1,],
                              point2 = dummy_space[2,]))

## The exploration on the same axis
summary(dispRity(dummy_space, metric = projections,
                              point1 = dummy_space[1,],
                              point2 = dummy_space[2,],
                              measure = "distance"))

By default, the vector (point1, point2) is used as unit vector of the projections (i.e. the Euclidean distance between (point1, point2) is set to 1) meaning that a projection value ("distance" or "position") of X means X times the distance between point1 and point2. If you want use the unit vector of the input matrix or are using a space where Euclidean distances are non-sensical, you can remove this option using scale = FALSE:

## The elaboration on the same axis using the dummy_space's
## unit vector
summary(dispRity(dummy_space, metric = projections,
                              point1 = dummy_space[1,],
                              point2 = dummy_space[2,],
                              scale = FALSE))

The projections.tree is the same as the projections metric but allows to determine the vector ((point1, point2)) using a tree rather than manually entering these points. The function intakes the exact same options as the projections function described above at the exception of point1 and point2. Instead it takes a the argument type that designates the type of vector to draw from the data based on a phylogenetic tree phy. The argument type can be a pair of any of the following inputs:

"root": to automatically use the coordinates of the root of the tree (the first element in phy$node.label);
"ancestor": to automatically use the coordinates of the elements' (i.e. any row in the matrix) most recent ancestor;
"tips": to automatically use the coordinates from the centroid of all tips;
"nodes": to automatically use the coordinates from the centroid of all nodes;
"livings": to automatically use the coordinates from the centroid of all "living" tips (i.e. the tips that are the furthest away from the root);
"fossils": to automatically use the coordinates from the centroid of all "fossil" tips and nodes (i.e. not the "living" ones);
any numeric values that can be interpreted as point1 and point2 in projections (e.g. 0, c(0, 1.2, 3/4), etc.);
or a user defined function that with the inputs matrix and phy and row (the element's ID, i.e. the row number in matrix).

For example, if you want to measure the projection of each element in the matrix (tips and nodes) on the axis from the root of the tree to each element's most recent ancestor, you can define the vector as type = c("root", "ancestor").

## Adding a extra row to dummy matrix (to match dummy_tree)
tree_space <- rbind(dummy_space, root = rnorm(5))
## Creating a random dummy tree (with labels matching the ones from tree_space)
dummy_tree <- rtree(6)
dummy_tree$tip.label <- rownames(tree_space)[1:6]
dummy_tree$node.label <- rownames(tree_space)[rev(7:11)]

## Measuring the disparity as the projection of each element
## on its root-ancestor vector
summary(dispRity(tree_space, metric = projections.tree,
                             tree   = dummy_tree,
                             type   = c("root", "ancestor")))

Of course you can also use any other options from the projections function:

## A user defined function that's returns the centroid of
## the first three nodes
fun.root <- function(matrix, tree, row = NULL) {
  return(colMeans(matrix[tree$node.label[1:3], ]))
}
## Measuring the unscaled rejection from the vector from the
## centroid of the three first nodes
## to the coordinates of the first tip
summary(dispRity(tree_space, metric  = projections.tree,
                             tree    = dummy_tree,
                             measure = "distance",
                             type    = list(fun.root,
                                            tree_space[1, ])))

Roundness

The roundness coefficient (or metric) ranges between 0 and 1 and expresses the distribution of and ellipse' major axis ranging from 1, a totally round ellipse (i.e. a circle) to 0 a totally flat ellipse (i.e. a line). A value of $0.5$ represents a regular ellipse where each major axis is half the size of the previous major axis. A value $> 0.5$ describes a pancake where the major axis distribution is convex (values close to 1 can be pictured in 3D as a cr`{e}pes with the first two axis being rather big - a circle - and the third axis being particularly thin; values closer to $0.5$ can be pictured as flying saucers). Conversely, a value $< 0.5$ describes a cigar where the major axis distribution is concave (values close to 0 can be pictured in 3D as a spaghetti with the first axis rather big and the two next ones being small; values closer to $0.5$ can be pictured in 3D as a fat cigar).

This is what it looks for example for three simulated variance-covariance matrices in 3D:

## Note that the plotting functions requires the ellipse, rgl and VCVtools packages (but not need to load them)
require("ellipse")
require("rgl")
require("VCVtools")
## To install them
# install.packages(c("ellipse", "rgl"))
# devtools::install_github("TGuillerme/VCVtools")

## knitr RGL setup
# knitr::knit_hooks$set(webgl = hook_webgl)
options(rgl.useNULL = TRUE) # Suppress the separate window.

## A 3D sphere
sphere  <- VCVtools::make.VCV(shape = 1, dimensions = 3, min.thick = 0.01)
pancake <- VCVtools::make.VCV(shape = 0.75, dimensions = 3, min.thick = 0.01)
cigar   <- VCVtools::make.VCV(shape = 0.25, dimensions = 3, min.thick = 0.01)

# options(rgl.printRglwidget = TRUE)
VCVtools::plot.VCV(sphere, dimensions = c(1,2,3), col.major.axes = c("orange", "blue", "darkgreen"), main = "Sphere's roundness: 1", lwd.major.axes = 3)
rglwidget()

VCVtools::plot.VCV(pancake, dimensions = c(1,2,3), col.major.axes = c("orange", "blue", "darkgreen"), main = "Pancake's roundness: ~0.70", lwd.major.axes = 3)
rglwidget()

VCVtools::plot.VCV(cigar, dimensions = c(1,2,3), col.major.axes = c("orange", "blue", "darkgreen"), main = "Cigar's roundness: ~0.25", lwd.major.axes = 3)
rglwidget()

par(mfrow = c(1,3))
VCVtools::plot.roundness(sphere, dimensions = c(1,2,3), col.lines = c("orange", "blue", "darkgreen"), main = "Sphere's roundness: 1")
VCVtools::plot.roundness(pancake, dimensions = c(1,2,3), col.lines = c("orange", "blue", "darkgreen"), main = "Pancake's roundness: ~0.70")
VCVtools::plot.roundness(cigar, dimensions = c(1,2,3), col.lines = c("orange", "blue", "darkgreen"), main = "Cigar's roundness: ~0.25")

Between group metrics {#betweengroupmetricslist}

You can find detailed explanation on how between group metrics work here.

`group.dist`

The group.dist metric allows to measure the distance between two groups in the multidimensional space. This function needs to intake several groups and use the option between.groups = TRUE in the dispRity function. It calculates the vector normal distance (euclidean) between two groups and returns 0 if that distance is negative. Note that it is possible to set up which quantiles to consider for calculating the distances between groups. For example, one might be interested in only considering the 95% CI for each group. This can be done through the option probs = c(0.025, 0.975) that is passed to the quantile function. It is also possible to use this function to measure the distance between the groups centroids by calculating the 50% quantile (probs = c(0.5)).

## Creating a dispRity object with two groups
grouped_space <- custom.subsets(dummy_space,
                      group = list(c(1:5), c(6:10)))

## Measuring the minimum distance between both groups
summary(dispRity(grouped_space, metric = group.dist,
                 between.groups = TRUE))

## Measuring the centroid distance between both groups
summary(dispRity(grouped_space, metric = group.dist,
                 between.groups = TRUE, probs = 0.5))

## Measuring the distance between both group's 75% CI
summary(dispRity(grouped_space, metric = group.dist,
                 between.groups = TRUE, probs = c(0.25, 0.75)))

`point.dist`

The metric measures the distance between the elements in one group (matrix) and a point calculated from a second group (matrix2). By default this point is the centroid but can be any point defined by a function passed to the point argument. For example, the centroid of matrix2 is the mean of each column of that matrix so point = colMeans (default). This function also takes the method argument like previous one described above to measure either the "euclidean" (default) or the "manhattan" distances:

## Measuring the distance between the elements of the first group
## and the centroid of the second group
summary(dispRity(grouped_space, metric = point.dist,
                 between.groups = TRUE))

## Measuring the distance between the elements of the second group
## and the centroid of the first group
summary(dispRity(grouped_space, metric = point.dist,
                 between.groups = list(c(2,1))))

## Measuring the distance between the elements of the first group
## a point defined as the standard deviation of each column
## in the second group
sd.point <- function(matrix2) {apply(matrix2, 2, sd)}
summary(dispRity(grouped_space, metric = point.dist,
                 point = sd.point, method = "manhattan",
                 between.groups = TRUE))

`projections.between` and `disalignment`

These two metrics are typically based on variance-covariance matrices from a dispRity object that has a $covar component (see more about that here). Both are based on the projections metric and can take the same optional arguments (more info here). The examples and explanations below are based on the default arguments but it is possible (and easy!) to change them.

We are going to use the charadriiformes example for both metrics (see more about that here).

## Loading the charadriiformes data
data(charadriiformes)

## Creating the dispRity object (see the #covar section in the manual for more info)
my_covar <- MCMCglmm.subsets(n = 50,
                             data = charadriiformes$data,
                             posteriors = charadriiformes$posteriors,
                             group = MCMCglmm.levels(charadriiformes$posteriors)[1:4],
                             tree = charadriiformes$tree,
                             rename.groups = c(levels(charadriiformes$data$clade), "phylogeny"))

The first metric, projections.between projects the major axis of one group (matrix) onto the major axis of another one (matrix2). For example we might want to know how some groups compare in terms of angle (orientation) to a base group:

## Creating the list of groups to compare
comparisons_list <- list(c("gulls", "phylogeny"),
                         c("plovers", "phylogeny"),
                         c("sandpipers", "phylogeny"))

## Measuring the angles between each groups
## (note that we set the metric as.covar, more on that in the #covar section below)
groups_angles <- dispRity(data = my_covar,
                          metric = as.covar(projections.between),
                          between.groups = comparisons_list,
                          measure = "degree")
## And here are the angles in degrees:
summary(groups_angles)

The second metric, disalignment rejects the centroid of a group (matrix) onto the major axis of another one (matrix2). This allows to measure wether the center of a group is aligned with the major axis of another. A disalignement value of 0 means that the groups are aligned. A higher disalignment value means the groups are more and more disaligned. We can use the same set of comparisons as in the projections.between examples to measure which group is most aligned (less disaligned) with the phylogenetic major axis:

## Measuring the disalignement of each group
groups_alignement <- dispRity(data = my_covar,
                              metric = as.covar(disalignment),
                              between.groups = comparisons_list)
## And here are the groups alignment (0 = aligned)
summary(groups_alignement)

Which disparity metric to choose?

The disparity metric that gives the most consistent results is the following one:

best.metric <- function() return(42)

Joke aside, this is a legitimate question that has no simple answer: it depends on the dataset and question at hand. Thoughts on which metric to choose can be find in @moms and @Guillerme2020 but again, will ultimately depend on the question and dataset. The question should help figuring out which type of metric is desired: for example, in the question "does the extinction released niches for mammals to evolve", the metric in interest should probably pick up a change in size in the trait space (the release could result in some expansion of the mammalian morphospace); or if the question is "does group X compete with group Y", maybe the metric of interested should pick up changes in position (group X can be displaced by group Y).

In order to visualise what signal different disparity metrics are picking, you can use the moms that come with a detailed manual on how to use it.

Alternatively, you can use the test.metric function:

`test.metric` {#test-metric}

This function allows to test whether a metric picks different changes in disparity. It intakes the space on which to test the metric, the disparity metric and the type of changes to apply gradually to the space. Basically this is a type of biased data rarefaction (or non-biased for "random") to see how the metric reacts to specific changes in trait space.

## Creating a 2D uniform space
example_space <- space.maker(300, 2, runif)

## Testing the product of ranges metric on the example space
example_test <- test.metric(example_space, metric = c(prod, ranges),
                           shifts = c("random", "size"))

By default, the test runs three replicates of space reduction as described in @moms by gradually removing 10% of the data points following the different algorithms from @moms (here the "random" reduction and the "size") reduction, resulting in a dispRity object that can be summarised or plotted. The number of replicates can be changed using the replicates option. Still by default, the function then runs a linear model on the simulated data to measure some potential trend in the changes in disparity. The model can be changed using the model option. Finally, the function runs 10 reductions by default from keeping 10% of the data (removing 90%) and way up to keeping 100% of the data (removing 0%). This can be changed using the steps option. A good disparity metric for your dataset will typically have no trend in the "random" reduction (the metric is ideally not affected by sample size) but should have a trend for the reduction of interest.

## The results as a dispRity object
example_test

## Summarising these results
summary(example_test)

## Or visualising them
plot(example_test)

Summarising dispRity data (plots)

Because of its architecture, printing dispRity objects only summarises their content but does not print the disparity value measured or associated analysis (more about this here). To actually see what is in a dispRity object, one can either use the summary function for visualising the data in a table or plot to have a graphical representation of the results.

Summarising `dispRity` data

This function is an S3 function (summary.dispRity) allowing users to summarise the content of dispRity objects that contain disparity calculations.

## Example data from previous sections
crown_stem <- custom.subsets(BeckLee_mat50,
                             group = crown.stem(BeckLee_tree,
                             inc.nodes = FALSE))
## Bootstrapping and rarefying these groups
boot_crown_stem <- boot.matrix(crown_stem, bootstraps = 100,
                               rarefaction = TRUE)
## Calculate disparity
disparity_crown_stem <- dispRity(boot_crown_stem,
                                 metric = c(sum, variances))

## Creating time slice subsets
time_slices <- chrono.subsets(data = BeckLee_mat99,
                              tree = BeckLee_tree,
                              method = "continuous",
                              model = "proximity",
                              time = c(120, 80, 40, 0),
                              FADLAD = BeckLee_ages)
## Bootstrapping the time slice subsets
boot_time_slices <- boot.matrix(time_slices, bootstraps = 100)
## Calculate disparity
disparity_time_slices <- dispRity(boot_time_slices,
                                  metric = c(sum, variances))

## Creating time bin subsets
time_bins <- chrono.subsets(data = BeckLee_mat99,
                            tree = BeckLee_tree, 
                            method = "discrete",
                            time = c(120, 80, 40, 0),
                            FADLAD = BeckLee_ages,
                            inc.nodes = TRUE)
## Bootstrapping the time bin subsets
boot_time_bins <- boot.matrix(time_bins, bootstraps = 100)
## Calculate disparity
disparity_time_bins <- dispRity(boot_time_bins,
                                metric = c(sum, variances))

These objects are easy to summarise as follows:

## Default summary
summary(disparity_time_slices)

Information about the number of elements in each subset and the observed (i.e. non-bootstrapped) disparity are also calculated. This is specifically handy when rarefying the data for example:

head(summary(disparity_crown_stem))

The summary functions can also take various options such as:

quantiles values for the confidence interval levels (by default, the 50 and 95 quantiles are calculated)
cent.tend for the central tendency to use for summarising the results (default is median)
digits option corresponding to the number of decimal places to print (default is 2)
recall option for printing the call of the dispRity object as well (default is FALSE)

These options can easily be changed from the defaults as follows:

## Same as above but using the 88th quantile and the standard deviation as the summary 
summary(disparity_time_slices, quantiles = 88, cent.tend = sd)

## Printing the details of the object and digits the values to the 5th decimal place
summary(disparity_time_slices, recall = TRUE, digits = 5)

Note that the summary table is a data.frame, hence it is as easy to modify as any dataframe using dplyr. You can also export it in csv format using write.csv or write_csv or even directly export into LaTeX format using the following;

## Loading the xtable package
require(xtable)
## Converting the table in LaTeX
xtable(summary(disparity_time_slices))

Plotting `dispRity` data

An alternative (and more fun!) way to display the calculated disparity is to plot the results using the S3 method plot.dispRity. This function takes the same options as summary.dispRity along with various graphical options described in the function help files (see ?plot.dispRity).

The plots can be of five different types:

preview for a 2d preview of the trait-space.
continuous for displaying continuous disparity curves
box, lines, and polygons to display discrete disparity results in respectively a boxplot, confidence interval lines, and confidence interval polygons.

This argument can be left empty. In this case, the algorithm will automatically detect the type of subsets from the dispRity object and plot accordingly.

It is also possible to display the number of elements in each subset (as a horizontal dotted line) using the option elements = TRUE. Additionally, when the data is rarefied, one can indicate which level of rarefaction to display (i.e. only display the results for a certain number of elements) by using the rarefaction argument.

## Graphical parameters
op <- par(mfrow = c(2, 2), bty = "n")

## Plotting continuous disparity results
plot(disparity_time_slices, type = "continuous")

## Plotting discrete disparity results
plot(disparity_crown_stem, type = "box")

## As above but using lines for the rarefaction level of 20 elements only
plot(disparity_crown_stem, type = "line", rarefaction = 20)

## As above but using polygons while also displaying the number of elements
plot(disparity_crown_stem, type = "polygon", elements = TRUE)

## Resetting graphical parameters
par(op)

Since plot.dispRity uses the arguments from the generic plot method, it is of course possible to change pretty much everything using the regular plot arguments:

## Graphical options
op <- par(bty = "n")

## Plotting the results with some classic options from plot
plot(disparity_time_slices, col = c("blue", "orange", "green"),
    ylab = c("Some measurement"), xlab = "Some other measurement",
    main = "Many options...", ylim = c(10, 0), xlim = c(4, 0))

## Adding a legend
legend("topleft", legend = c("Central tendency",
                             "Confidence interval 1",
                             "Confidence interval 2"),
      col = c("blue", "orange", "green"), pch = 19)

## Resetting graphical parameters
par(op)

In addition to the classic plot arguments, the function can also take arguments that are specific to plot.dispRity like adding the number of elements or rarefaction level (as described above), and also changing the values of the quantiles to plot as well as the central tendency.

## Graphical options
op <- par(bty = "n")

## Plotting the results with some plot.dispRity arguments
plot(disparity_time_slices,
    quantiles = c(seq(from = 10, to = 100, by = 10)),
    cent.tend = sd, type = "c", elements = TRUE,
    col = c("black", rainbow(10)),
    ylab = c("Disparity", "Diversity"),
    xlab = "Time (in in units from past to present)",
    observed = TRUE,
    main = "Many more options...")

## Resetting graphical parameters
par(op)

Note that the argument observed = TRUE allows to plot the disparity values calculated from the non-bootstrapped data as crosses on the plot.

For comparing results, it is also possible to add a plot to the existent plot by using add = TRUE:

## Graphical options
op <- par(bty = "n")

## Plotting the continuous disparity with a fixed y axis
plot(disparity_time_slices, ylim = c(3, 9))
## Adding the discrete data
plot(disparity_time_bins, type = "line", ylim = c(3, 9),
     xlab = "", ylab = "", add = TRUE)

## Resetting graphical parameters
par(op)

Finally, if your data has been fully rarefied, it is also possible to easily look at rarefaction curves by using the rarefaction = TRUE argument:

## Graphical options
op <- par(bty = "n")

## Plotting the rarefaction curves
plot(disparity_crown_stem, rarefaction = TRUE)

## Resetting graphical parameters
par(op)

`type = preview`

Note that all the options above are plotting disparity objects for which a disparity metric has been calculated. This makes totally sense for dispRity objects but sometimes it might be interesting to look at what the trait-space looks like before measuring the disparity. This can be done by plotting dispRity objects with no calculated disparity!

For example, we might be interested in looking at how the distribution of elements change as a function of the distributions of different sub-settings. For example custom subsets vs. time subsets:

## Making the different subsets
cust_subsets <- custom.subsets(BeckLee_mat99,
                               crown.stem(BeckLee_tree,
                                          inc.nodes = TRUE))
time_subsets <- chrono.subsets(BeckLee_mat99,
                               tree = BeckLee_tree,
                               method = "discrete",
                               time = 5)

## Note that no disparity has been calculated here:
is.null(cust_subsets$disparity)
is.null(time_subsets$disparity)

## But we can still plot both spaces by using the default plot functions
par(mfrow = c(1,2))
## Default plotting
plot(cust_subsets)
## Plotting with more arguments
plot(time_subsets, specific.args = list(dimensions = c(1,2)),
     main = "Some \"low\" dimensions")

DISCLAIMER: This functionality can be handy for exploring the data (e.g. to visually check whether the subset attribution worked) but it might be misleading on how the data is actually distributed in the multidimensional space! Groups that don't overlap on two set dimensions can totally overlap in all other dimensions!

For dispRity objects that do contain disparity data, the default option is to plot your disparity data. However you can always force the preview option using the following:

par(mfrow = c(2,1))
## Default plotting
plot(disparity_time_slices, main = "Disparity through time")
## Plotting with more arguments
plot(disparity_time_slices, type = "preview",
     main = "Two first dimensions of the trait space")

Graphical options with `...`

As mentioned above all the plots using plot.dispRity you can use the ... options to add any type of graphical parameters recognised by plot. However, sometimes, plotting more advanced "dispRity" objects also calls other generic functions such as lines, points or legend. You can fine tune which specific function should be affected by ... by using the syntax <function>.<argument> where <function> is usually the function to plot a specific element in the plot (e.g. points) and the <argument> is the specific argument you want to change for that function. For example, in a plot containing several elements, including circles (plotted internally with points), you can decide to colour everything in blue using the normal col = "blue" option. But you can also decide to only colour the circles in blue using points.col = "blue"!

Here is an example with multiple elements (lines and points) taken from the disparity with trees section below:

## Loading some demo data:
## An ordinated matrix with node and tip labels
data(BeckLee_mat99)
## The corresponding tree with tip and node labels
data(BeckLee_tree)
## A list of tips ages for the fossil data
data(BeckLee_ages)

## Time slicing through the tree using the equal split algorithm
time_slices <- chrono.subsets(data   = BeckLee_mat99,
                              tree   = BeckLee_tree,
                              FADLAD = BeckLee_ages,
                              method = "continuous",
                              model  = "acctran",
                              time   = 15)

par(mfrow = c(2,2))
## The preview plot with the tree using only defaults
plot(time_slices, type = "preview", specific.args = list(tree = TRUE))
## The same plot but by applying general options
plot(time_slices, type = "preview", specific.args = list(tree = TRUE),
     col = "blue", main = "General options")
## The same plot but by applying the colour only to the lines
## and change of shape only to the points
plot(time_slices, type = "preview", specific.args = list(tree = TRUE),
     lines.col = "blue", points.pch = 15, main = "Specific options")
## And now without the legend
plot(time_slices, type = "preview", specific.args = list(tree = TRUE),
     lines.col = "blue", points.pch = 15, legend = FALSE)

Testing disparity hypotheses

The dispRity package allows users to apply statistical tests to the calculated disparity to test various hypotheses. The function test.dispRity works in a similar way to the dispRity function: it takes a dispRity object, a test and a comparisons argument.

The comparisons argument indicates the way the test should be applied to the data:

pairwise (default): to compare each subset in a pairwise manner
referential: to compare each subset to the first subset
sequential: to compare each subset to the following subset
all: to compare all the subsets together (like in analysis of variance)

It is also possible to input a list of pairs of numeric values or characters matching the subset names to create personalised tests. Some other tests implemented in dispRity such as the dispRity::null.test have a specific way they are applied to the data and therefore ignore the comparisons argument.

The test argument can be any statistical or non-statistical test to apply to the disparity object. It can be a common statistical test function (e.g. stats::t.test), a function implemented in dispRity (e.g. see ?null.test) or any function defined by the user.

This function also allows users to correct for Type I error inflation (false positives) when using multiple comparisons via the correction argument. This argument can be empty (no correction applied) or can contain one of the corrections from the stats::p.adjust function (see ?p.adjust).

Note that the test.dispRity algorithm deals with some classical test outputs (h.test, lm and numeric vector) and summarises the test output. It is, however, possible to get the full detailed output by using the options details = TRUE.

Here we are using the variables generated in the section above:

## T-test to test for a difference in disparity between crown and stem mammals
test.dispRity(disparity_crown_stem, test = t.test)

## Performing the same test but with the detailed t.test output
test.dispRity(disparity_crown_stem, test = t.test, details = TRUE)

## Wilcoxon test applied to time sliced disparity with sequential comparisons,
## with Bonferroni correction
test.dispRity(disparity_time_slices, test = wilcox.test,
              comparisons = "sequential", correction = "bonferroni")

## Measuring the overlap between distributions in the time bins (using the
## implemented Bhattacharyya Coefficient function - see ?bhatt.coeff)
test.dispRity(disparity_time_bins, test = bhatt.coeff)

Because of the modular design of the package, tests can always be made by the user (the same way disparity metrics can be user made). The only condition is that the test can be applied to at least two distributions. In practice, the test.dispRity function will pass the calculated disparity data (distributions) to the provided function in either pairs of distributions (if the comparisons argument is set to pairwise, referential or sequential) or a table containing all the distributions (comparisons = all; this should be in the same format as data passed to lm-type functions for example).

NPMANOVA in `dispRity` {#adonis}

One often useful test to apply to multidimensional data is the permutational multivariate analysis of variance based on distance matrices vegan::adonis2. This can be done on dispRity objects using the adonis.dispRity wrapper function. Basically, this function takes the exact same arguments as adonis and a dispRity object for data and performs a PERMANOVA based on the distance matrix of the multidimensional space (unless the multidimensional space was already defined as a distance matrix). The adonis.dispRity function uses the information from the dispRity object to generate default formulas:

If the object contains customised subsets, it applies the default formula matrix ~ group testing the effect of group as a predictor on matrix (called from the dispRity object as data$matrix see dispRity object details)
If the object contains time subsets, it applies the default formula matrix ~ time testing the effect of time as a predictor (were the different levels of time are the different time slices/bins)

set.seed(1)
## Generating a random character matrix
character_matrix <- sim.morpho(rtree(20), 50,
                               rates = c(rnorm, 1, 0))

## Calculating the distance matrix
distance_matrix <- as.matrix(dist(character_matrix))

## Creating two groups
random_groups <- list("group1" = 1:10, "group2" = 11:20)

## Generating a dispRity object
random_disparity <- custom.subsets(distance_matrix, random_groups)

## Running a default NPMANOVA
adonis.dispRity(random_disparity)

Of course, it is possible to pass customised formulas if the disparity object contains more more groups. In that case the predictors must correspond to the names of the groups explained data must be set as matrix:

## Creating two groups with two states each
groups <- as.data.frame(matrix(data = c(rep(1,10),
                                        rep(2,10),
                                        rep(c(1,2), 10)),
                        nrow = 20, ncol = 2,
                        dimnames = list(paste0("t", 1:20),
                                        c("g1", "g2"))))

## Creating the dispRity object
multi_groups <- custom.subsets(distance_matrix, groups)

## Running the NPMANOVA
adonis.dispRity(multi_groups, matrix ~ g1 + g2)

Finally, it is possible to use objects generated by chrono.subsets. In this case, adonis.dispRity will applied the matrix ~ time formula by default:

## Creating time series
time_subsets <- chrono.subsets(BeckLee_mat50, BeckLee_tree,
                               method = "discrete",
                               inc.nodes = FALSE,
                               time = c(100, 85, 65, 0),
                               FADLAD = BeckLee_ages)

## Running the NPMANOVA with time as a predictor
adonis.dispRity(time_subsets)

Note that the function warns you that the input data was transformed into a distance matrix. This is reflected in the Call part of the output (formula = dist(matrix) ~ time).

To use each time subset as a separate predictor, you can use the matrix ~ chrono.subsets formula; this is equivalent to matrix ~ first_time_subset + second_time_subset + ...:

## Running the NPMANOVA with each time bin as a predictor
adonis.dispRity(time_subsets, matrix ~ chrono.subsets)

`geiger::dtt` model fitting in `dispRity` {#dtt}

The dtt function from the geiger package is also often used to compare a trait's disparity observed in living taxa to the disparity of a simulated trait based on a given phylogeny. The dispRity package proposes a wrapper function for geiger::dtt, dtt.dispRity that allows the use of any disparity metric. Unfortunately, this implementation is slower that geiger::dtt (so if you're using the metrics implemented in geiger prefer the original version) and, as the original function, is limited to ultrametric trees (only living taxa!)...

require(geiger)
geiger_data <- get(data(geospiza))

## Calculate the disparity of the dataset using the sum of variance
dispRity_dtt <- dtt.dispRity(data = geiger_data$dat,
                             metric = c(sum, variances),
                             tree = geiger_data$phy,
                             nsim = 100)

## Plotting the results
plot(dispRity_dtt)

Note that, like in the original dtt function, it is possible to change the evolutionary model (see ?geiger::sim.char documentation).

null morphospace testing with `null.test` {#null-test}

This test is equivalent to the test performed in @diaz2016global. It compares the disparity measured in the observed space to the disparity measured in a set of simulated spaces. These simulated spaces can be built with based on the hypothesis assumptions: for example, we can test whether our space is normal.

set.seed(123)
## A "normal" multidimensional space with 50 dimensions and 10 elements
normal_space <- matrix(rnorm(1000), ncol = 50)

## Calculating the disparity as the average pairwise distances
obs_disparity <- dispRity(normal_space,
                          metric = c(mean, pairwise.dist))

## Testing against 100 randomly generated normal spaces
(results <- null.test(obs_disparity, replicates = 100,
                      null.distrib = rnorm))

Here the results show that disparity measured in our observed space is not significantly different than the one measured in a normal space. We can then propose that our observed space is normal!

These results have an attributed dispRity and randtest class and can be plotted as randtest objects using the dispRity S3 plot method:

## Plotting the results
plot(results, main = "Is this space normal?")

For more details on generating spaces see the space.maker function tutorial.

Fitting modes of evolution to disparity data {#model-fitting}

The code used for these models is based on those developed by Gene Hunt [@hunt2006fitting; @hunt2012measuring; @hunt2015simple]. So we acknowledge and thank Gene Hunt for developing these models and writing the original R code that served as inspiration for these models.

DISCLAIMER: this method of analysing disparity has not been published yet and has not been peer reviewed. Caution should be used in interpreting these results: it is unclear what "a disparity curve fitting a Brownian motion" actually means biologically.

As Malcolm said in Jurassic Park: "although the examples within this chapter all work and produce solid tested results (from an algorithm point of view), that doesn't mean you should use it" (or something along those lines).

Simple modes of disparity change through time

`model.test`

Changes in disparity-through-time can follow a range of models, such as random walks, stasis, constrained evolution, trends, or an early burst model of evolution. We will start with by fitting the simplest modes of evolution to our data. For example we may have a null expectation of time-invariant change in disparity in which values fluctuate with a variance around the mean - this would be best describe by a Stasis model:

## Loading premade disparity data
data(BeckLee_disparity)
disp_time <- model.test(data = BeckLee_disparity, model = "Stasis")

We can see the standard output from model.test. The first output message tells us it has tested for equal variances in each sample. The model uses Bartlett's test of equal variances to assess if variances are equal, so if p > 0.05 then variance is treated as the same for all samples, but if (p < 0.05) then each bin variance is unique. Here we have p < 0.05, so variance is not pooled between samples.

By default model.test will use Bartlett's test to assess for homogeneity of variances, and then use this to decide to pool variances or not. This is ignored if the argument pool.variance in model.test is changed from the default NULL to TRUE or FALSE. For example, to ignore Bartlett's test and pool variances manually we would do the following:

disp_time_pooled <- model.test(data = BeckLee_disparity,
                               model = "Stasis",
                               pool.variance = TRUE)

However, unless you have good reason to choose otherwise it is recommended to use the default of pool.variance = NULL:

disp_time <- model.test(data = BeckLee_disparity,
                        model = "Stasis",
                        pool.variance = NULL)
disp_time

The remaining output gives us the log-likelihood of the Stasis model of r summary(disp_time)[,"log.lik"] (you may notice this change when we pooled variances above). The output also gives us the small sample Akaike Information Criterion (AICc), the delta AICc (the distance from the best fitting model), and the AICc weights (~the relative support of this model compared to all models, scaled to one).

These are all metrics of relative fit, so when we test a single model they are not useful. By using the function summary in dispRity we can see the maximum likelihood estimates of the model parameters:

summary(disp_time)

So we again see the AICc, delta AICc, AICc weight, and the log-likelihood we saw previously. We now also see the number of parameters from the model (2: theta and omega), and their estimates so the variance (omega = r summary(disp_time)[,"omega"]) and the mean (theta.1 = r summary(disp_time)[,"theta.1"]).

The model.test function is designed to test relative model fit, so we need to test more than one model to make relative comparisons. So let's compare to the fit of the Stasis model to another model with two parameters: the Brownian motion. Brownian motion assumes a constant mean that is equal to the ancestral estimate of the sequence, and the variance around this mean increases linearly with time. The easier way to compare these models is to simply add "BM" to the models vector argument:

disp_time <- model.test(data = BeckLee_disparity,
                        model = c("Stasis", "BM"))
disp_time

Et voilà! Here we can see by the log-likelihood, AICc, delta AICc, and AICc weight Brownian motion has a much better relative fit to these data than the Stasis model. Brownian motion has a relative AICc fitr summary(disp_time)[,"delta_aicc"][[1]] units better than Stasis, and has a AICc weight of 1.

We can also all the information about the relative fit of models alongside the maximum likelihood estimates of model parameters using the summary function

summary(disp_time)

Not that because the parameters per models differ, the summary includes NA for inapplicable parameters per models (e.g. the theta and omega parameters from the Stasis models are inapplicable for a Brownian motion model).

We can plot the relative fit of our models using the plot function

plot(disp_time)

Here we see and overwhelming support for the Brownian motion model.

Alternatively, we could test all available models single modes: Stasis, Brownian motion, Ornstein-Uhlenbeck (evolution constrained to an optima), Trend (increasing or decreasing mean through time), and Early Burst (exponentially decreasing rate through time)

disp_time <- model.test(data = BeckLee_disparity,
                  model = c("Stasis", "BM", "OU", "Trend", "EB"))
summary(disp_time)

These models indicate support for a Trend model, and we can plot the relative support of all model AICc weights.

plot(disp_time)

Note that although AIC values are indicator of model best fit, it is also important to look at the parameters themselves. For example OU can be really well supported but with an alpha parameter really close to 0, making it effectively a BM model [@Cooper2016].

Is this a trend of increasing or decreasing disparity through time? One way to find out is to look at the summary function for the Trend model:

summary(disp_time)["Trend",]

This show a positive trend (r unname(summary(disp_time)["Trend",]["trend"])) of increasing disparity through time.

Plot and run simulation tests in a single step

`model.test.wrapper`

Patterns of evolution can be fit using model.test, but the model.test.wrapper fits the same models as model.test as well as running predictive tests and plots.

The predictive tests use the maximum likelihood estimates of model parameters to simulate a number of datasets (default = 1000), and analyse whether this is significantly different to the empirical input data using the Rank Envelope test [@murrell2018global]. Finally we can plot the empirical data, simulated data, and the Rank Envelope test p values. This can all be done using the function model.test.wrapper, and we will set the argument show.p = TRUE so p values from the Rank Envelope test are printed on the plot:

disp_time <- model.test.wrapper(data = BeckLee_disparity,
                    model = c("Stasis", "BM", "OU", "Trend", "EB"),
                                show.p = TRUE)
disp_time

From this plot we can see the empirical estimates of disparity through time (pink) compared to the predictive data based upon the simulations using the estimated parameters from each model. There is no significant differences between the empirical data and simulated data, except for the Early Burst model.

Trend is the best-fitting model but the plot suggests the OU model also follows a trend-like pattern. This is because the optima for the OU model (r disp_time["OU", "optima.1"]) is different to the ancestral state (r disp_time["OU", "ancestral state"]) and outside the observed value. This is potentially unrealistic, and one way to alleviate this issue is to set the optima of the OU model to equal the ancestral estimate - this is the normal practice for OU models in comparative phylogenetics. To set the optima to the ancestral value we change the argument fixed.optima = TRUE:

disp_time <- model.test.wrapper(data = BeckLee_disparity,
                    model = c("Stasis", "BM", "OU", "Trend", "EB"),
                                show.p = TRUE, fixed.optima = TRUE)
disp_time

The relative fit of the OU model is decreased by constraining the fit of the optima to equal the ancestral state value. In fact as the OU attraction parameter (alpha) is zero, the model is equal to a Brownian motion model but is penalised by having an extra parameter. Note that indeed, the plots of the BM model and the OU model look nearly identical.

Multiple modes of evolution (time shifts)

As well as fitting a single model to a sequence of disparity values we can also allow for the mode of evolution to shift at a single or multiple points in time. The timing of a shift in mode can be based on an a prior expectation, such as a mass extinction event, or the model can test multiple points to allow to find time shift point with the highest likelihood.

Models can be fit using model.test but it can be more convenient to use model.test.wrapper. Here we will compare the relative fit of Brownian motion, Trend, Ornstein-Uhlenbeck and a multi-mode Ornstein Uhlenbck model in which the optima changes at 66 million years ago, the Cretaceous-Palaeogene boundary.

For example, we could be testing the hypothesis that the extinction of non-avian dinosaurs allowed mammals to go from scurrying in the undergrowth (low optima/low disparity) to dominating all habitats (high optima/high disparity). We will constrain the optima of OU model in the first time begin (i.e, pre-66 Mya) to equal the ancestral value:

disp_time <- model.test.wrapper(data = BeckLee_disparity,
                        model = c("BM", "Trend", "OU", "multi.OU"),
                                time.split = 66,
                                pool.variance = NULL,
                                show.p = TRUE,
                                fixed.optima = TRUE)
disp_time

The multi-OU model shows an increase an optima at the Cretaceous-Palaeogene boundary, indicating a shift in disparity. However, this model does not fit as well as a model in which there is an increasing trend through time. We can also fit a model in which the we specify a heterogeneous model but we do not give a time.split. In this instance the model will test all splits that have at least 10 time slices on either side of the split. That's 102 potential time shifts in this example dataset so be warned, the following code will estimate 105 models!

## An example of a time split model in which all potential splits are tested
## WARNING: this will take between 20 minutes and half and hour to run!
disp_time <- model.test.wrapper(data = BeckLee_disparity,
                        model = c("BM", "Trend", "OU", "multi.OU"),
                                show.p = TRUE, fixed.optima = TRUE)

As well as specifying a multi-OU model we can run any combination of models. For example we could fit a model at the Cretaceous-Palaeogene boundary that goes from an OU to a BM model, a Trend to an OU model, a Stasis to a Trend model or any combination you want to use. The only model that can't be used in combination is a multi-OU model.

These can be introduced by changing the input for the models into a list, and supplying a vector with the two models. This is easier to see with an example:

## The models to test
my_models <- list(c("BM", "OU"),
                  c("Stasis", "OU"),
                  c("BM", "Stasis"),
                  c("OU", "Trend"),
                  c("Stasis", "BM"))

## Testing the models
disp_time <- model.test.wrapper(data = BeckLee_disparity,
                                model = my_models, time.split = 66,
                                show.p = TRUE, fixed.optima = TRUE)
disp_time

`model.test.sim`

Note that all the models above where run using the model.test.wrapper function that is a... wrapping function! In practice, this function runs two main functions from the dispRity package and then plots the results:

model.test and
model.test.sim

The model.test.sim allows to simulate disparity evolution given a dispRity object input (as in model.test.wrapper) or given a model and its specification. For example, it is possible to simulate a simple Brownian motion model (or any of the other models or models combination described above):

## A simple BM model
model_simulation <- model.test.sim(sim = 1000, model = "BM",
                                   time.span = 50, variance = 0.1,
                                   sample.size = 100,
                                   parameters = list(ancestral.state = 0))
model_simulation

This will simulate 1000 Brownian motions for 50 units of time with 100 sampled elements, a variance of 0.1 and an ancestral state of 0. We can also pass multiple models in the same way we did it for model.test This model can then be summarised and plotted as most dispRity objects:

## Displaying the 5 first rows of the summary
head(summary(model_simulation))

## Plotting the simulations
plot(model_simulation)

Note that these functions can take all the arguments that can be passed to plot, summary, plot.dispRity and summary.dispRity.

Simulating tested models

Maybe more interestingly though, it is possible to pass the output of model.test directly to model.test.sim to simulate the models that fits the data the best and calculate the Rank Envelope test p value. Let's see that using the simple example from the start:

## Fitting multiple models on the data set
disp_time <- model.test(data = BeckLee_disparity,
                    model = c("Stasis", "BM", "OU", "Trend", "EB"))
summary(disp_time)

As seen before, the Trend model fitted this dataset the best. To simulate what 1000 Trend models would look like using the same parameters as the ones estimated with model.test (here the ancestral state being r summary(disp_time)["Trend", "ancestral state"], the sigma squared being r summary(disp_time)["Trend", "sigma squared"] and the trend of r summary(disp_time)["Trend", "trend"]), we can simply pass this model to model.test.sim:

## Simulating 1000 Trend model with the observed parameters
sim_trend <- model.test.sim(sim = 1000, model = disp_time)
sim_trend

By default, the model simulated is the one with the lowest AICc (model.rank = 1) but it is possible to choose any ranked model, for example, the OU (second one):

## Simulating 1000 OU model with the observed parameters
sim_OU <- model.test.sim(sim = 1000, model = disp_time,
                         model.rank = 2)
sim_OU

And as the example above, the simulated data can be plotted or summarised:

head(summary(sim_trend))
head(summary(sim_OU))

## The trend model with some graphical options
plot(sim_trend, xlab = "Time (Mya)", ylab = "sum of variances",
    col = c("#F65205", "#F38336", "#F7B27E"))

## Adding the observed disparity through time
plot(BeckLee_disparity, add = TRUE, col = c("#3E9CBA", "#98D4CF90", "#BFE4E390"))

Disparity as a distribution {#disparity-distribution}

Disparity is often regarded as a summary value of the position of the all elements in the ordinated space. For example, the sum of variances, the product of ranges or the median distance between the elements and their centroid will summarise disparity as a single value. This value can be pseudo-replicated (bootstrapped) to obtain a distribution of the summary metric with estimated error. However, another way to perform disparity analysis is to use the whole distribution rather than just a summary metric (e.g. the variances or the ranges).

This is possible in the dispRity package by calculating disparity as a dimension-level 2 metric only! Let's have a look using our previous example of bootstrapped time slices but by measuring the distances between each taxon and their centroid as disparity.

## Measuring disparity as a whole distribution
disparity_centroids <- dispRity(boot_time_slices, 
                                metric = centroids)

The resulting disparity object is of dimension-level 2, so it can easily be transformed into a dimension-level 1 object by, for example, measuring the median distance of all these distributions:

## Measuring median disparity in each time slice
disparity_centroids_median <- dispRity(disparity_centroids,
                                       metric = median)

And we can now compare the differences between these methods:

## Summarising both disparity measurements:
## The distributions:
summary(disparity_centroids)
## The summary of the distributions (as median)
summary(disparity_centroids_median)

We can see that the summary message for the distribution is slightly different than before. Here summary also displays the observed central tendency (i.e. the central tendency of the measured distributions). Note that, as expected, this central tendency is the same in both metrics!

Another, maybe more intuitive way, to compare both approaches for measuring disparity is to plot the distributions:

## Graphical parameters
op <- par(bty = "n", mfrow = c(1, 2))

## Plotting both disparity measurements
plot(disparity_centroids,
     ylab = "Distribution of all the distances")
plot(disparity_centroids_median,
     ylab = "Distribution of the medians of all the distances")

par(op)

We can then test for differences in the resulting distributions using test.dispRity and the bhatt.coeff test as described above.

## Probability of overlap in the distribution of medians
test.dispRity(disparity_centroids_median, test = bhatt.coeff)

In this case, we are looking at the probability of overlap of the distribution of median distances from centroids among each pair of time slices. In other words, we are measuring whether the medians from each bootstrap pseudo-replicate for each time slice overlap. But of course, we might be interested in the actual distribution of the distances from the centroid rather than simply their central tendencies. This can be problematic depending on the research question asked since we are effectively comparing non-independent medians distributions (because of the pseudo-replication).

One solution, therefore, is to look at the full distribution:

## Probability of overlap for the full distributions
test.dispRity(disparity_centroids, test = bhatt.coeff)

These results show the actual overlap among all the measured distances from centroids concatenated across all the bootstraps. For example, when comparing the slices 120 and 80, we are effectively comparing the 5 $\times$ 100 distances (the distances of the five elements in slice 120 bootstrapped 100 times) to the 19 $\times$ 100 distances from slice 80. However, this can also be problematic for some specific tests since the n $\times$ 100 distances are also pseudo-replicates and thus are still not independent.

A second solution is to compare the distributions to each other for each replicate:

## Boostrapped probability of overlap for the full distributions
test.dispRity(disparity_centroids, test = bhatt.coeff,
              concatenate = FALSE)

These results show the median overlap among pairs of distributions in the first column (bhatt.coeff) and then the distribution of these overlaps among each pair of bootstraps. In other words, when two distributions are compared, they are now compared for each bootstrap pseudo-replicate, thus effectively creating a distribution of probabilities of overlap. For example, when comparing the slices 120 and 80, we have a mean probability of overlap of 0.28 and a probability between 0.18 and 0.43 in 50\% of the pseudo-replicates. Note that the quantiles and central tendencies can be modified via the conc.quantiles option.

Disparity from other matrices {#other-matrices}

In the example so far, disparity was measured from an ordinated multidimensional space (i.e. a PCO of the distances between taxa based on discrete morphological characters). This is a common approach in palaeobiology, morphometrics or ecology but ordinated matrices are not mandatory for the dispRity package! It is totally possible to perform the same analysis detailed above using other types of matrices as long as your elements are rows in your matrix.

For example, we can use the data set eurodist, an R inbuilt dataset that contains the distances (in km) between European cities. We can check for example, if Northern European cities are closer to each other than Southern ones:

## Making the eurodist data set into a matrix (rather than "dist" object)
eurodist <- as.matrix(eurodist)
eurodist[1:5, 1:5]

## The two groups of cities
Northern <- c("Brussels", "Calais", "Cherbourg", "Cologne", "Copenhagen",
              "Hamburg", "Hook of Holland", "Paris", "Stockholm")
Southern <- c("Athens", "Barcelona", "Geneva", "Gibraltar", "Lisbon", "Lyons",
              "Madrid", "Marseilles", "Milan", "Munich", "Rome", "Vienna")

## Creating the subset dispRity object
eurodist_subsets <- custom.subsets(eurodist, group = list("Northern" = Northern,
                                                        "Southern" = Southern))

## Bootstrapping and rarefying to 9 elements (the number of Northern cities)
eurodist_bs <- boot.matrix(eurodist_subsets, rarefaction = 9)

## Measuring disparity as the median distance from group's centroid
euro_disp <- dispRity(eurodist_bs, metric = c(median, centroids))

## Testing the differences using a simple wilcox.test
euro_diff <- test.dispRity(euro_disp, test = wilcox.test)
euro_diff_rar <- test.dispRity(euro_disp, test = wilcox.test, rarefaction = 9)

We can compare this approach to an ordination one:

## Ordinating the eurodist matrix (with 11 dimensions)
euro_ord <- cmdscale(eurodist, k = 11)

## Calculating disparity on the bootstrapped and rarefied subset data
euro_ord_disp <- dispRity(boot.matrix(custom.subsets(euro_ord, group =
        list("Northern" = Northern, "Southern" = Southern)), rarefaction = 9),
        metric = c(median, centroids))

## Testing the differences using a simple wilcox.test
euro_ord_diff <- test.dispRity(euro_ord_disp, test = wilcox.test)
euro_ord_diff_rar <- test.dispRity(euro_ord_disp, test = wilcox.test, rarefaction = 9)

And visualise the differences:

## Plotting the differences
par(mfrow = c(2,2), bty = "n")
## Plotting the normal disparity
plot(euro_disp, main = "Distance differences")
## Adding the p-value
text(1.5, 4000, paste0("p=",round(euro_diff[[2]][[1]], digit = 5)))
## Plotting the rarefied disparity
plot(euro_disp, rarefaction = 9, main = "Distance differences (rarefied)")
## Adding the p-value
text(1.5, 4000, paste0("p=",round(euro_diff_rar[[2]][[1]], digit = 5)))

## Plotting the ordinated disparity
plot(euro_ord_disp, main = "Ordinated differences")
## Adding the p-value
text(1.5, 1400, paste0("p=",round(euro_ord_diff[[2]][[1]], digit = 5) ))
## Plotting the rarefied disparity
plot(euro_ord_disp, rarefaction = 9, main = "Ordinated differences (rarefied)")
## Adding the p-value
text(1.5, 1400, paste0("p=",round(euro_ord_diff_rar[[2]][[1]], digit = 5) ))

As expected, the results are pretty similar in pattern but different in terms of scale. The median centroids distance is expressed in km in the "Distance differences" plots and in Euclidean units of variation in the "Ordinated differences" plots.

Disparity from multiple matrices (and multiple trees!) {#multi.input}

Since the version 1.4 of this package, it is possible to use multiple trees and multiple matrices in dispRity objects. To use multiple matrices, this is rather easy: just supply a list of matrices to any of the dispRity functions and, as long as they have the same size and the same rownames they will be handled as a distribution of matrices.

set.seed(1)
## Creating 3 matrices with 4 dimensions and 10 elements each (called t1, t2, t3, etc...)
matrix_list <- replicate(3, matrix(rnorm(40), 10, 4, dimnames = list(paste0("t", 1:10))),
                         simplify = FALSE)
class(matrix_list) # This is a list of matrices

## Measuring some disparity metric on one of the matrices
summary(dispRity(matrix_list[[1]], metric = c(sum, variances)))

## Measuring the same disparity metric on the three matrices
summary(dispRity(matrix_list, metric = c(sum, variances)))

As you can see, when measuring the sum of variances on multiple matrices, we now have a distribution of sum of variances rather than a single observed value.

Similarly as running disparity analysis using multiple matrices, you can run the chrono.subsets function using multiple trees. This can be useful if you want to use a tree posterior distribution rather than a single consensus tree. These trees can be passed to chrono.subsets as a "multiPhylo" object (with the same node and tip labels in each tree). First let's define a function to generate multiple trees with the same labels and root ages:

set.seed(1)
## Matches the trees and the matrices
## A bunch of trees
make.tree <- function(n, fun = rtree) {
    ## Make the tree
    tree <- fun(n)
    tree <- chronos(tree, quiet = TRUE,
                    calibration = makeChronosCalib(tree, age.min = 10, age.max = 10))
    class(tree) <- "phylo"
    ## Add the node labels
    tree$node.label <- paste0("n", 1:Nnode(tree))
    ## Add the root time
    tree$root.time <- max(tree.age(tree)$ages)
    return(tree)
}
trees <- replicate(3, make.tree(10), simplify = FALSE)
class(trees) <- "multiPhylo"
trees

We can now simulate some ancestral states for the matrices in the example above to have multiple matrices associated with the multiple trees.

## A function for running the ancestral states estimations
do.ace <- function(tree, matrix) {
    ## Run one ace
    fun.ace <- function(character, tree) {
        results <- ace(character, phy = tree)$ace
        names(results) <- paste0("n", 1:Nnode(tree))
        return(results)
    }
    ## Run all ace
    return(rbind(matrix, apply(matrix, 2, fun.ace, tree = tree)))
}

## All matrices
matrices <- mapply(do.ace, trees, matrix_list, SIMPLIFY = FALSE)

Let's first see an example of time-slicing with one matrix and multiple trees. This assumes that your tip values (observed) and node values (estimated) are fixed with no error on them. It also assumes that the nodes in the matrix always corresponds to the node in the trees (in other words, the tree topologies are fixed):

## Making three "proximity" time slices across one tree
one_tree <- chrono.subsets(matrices[[1]], trees[[1]],
                           method = "continuous",
                           model = "proximity", time = 3)
## Making three "proximity" time slices across the three trees
three_tree <- chrono.subsets(matrices[[1]], trees,
                             method = "continuous",
                             model = "proximity", time = 3)
## Measuring disparity as the sum of variances and summarising it
summary(dispRity(one_tree, metric = c(sum, variances)))
summary(dispRity(three_tree, metric = c(sum, variances)))

This results show the effect of considering a tree distribution: in the first case (one_tree) the time slice at 3.95 Mya has a sum of variances of 2.9 but this values goes down to 0.256 in the second case (three_tree) which is due to the differences in branch lengths distributions:

par(mfrow = c(3,1))
slices <- c(7.9, 3.95, 0)
fun.plot <- function(tree) {
  plot(tree)
  nodelabels(tree$node.label, cex = 0.8)
  axisPhylo()
  abline(v = tree$root.time - slices)
}
silent <- lapply(trees, fun.plot)

Note that in this example, the nodes are actually even different in each tree! The node n4 for example, is not direct descendent of t4 and t6 in all trees! To fix that, it is possible to input a list of trees and a list of matrices that correspond to each tree in chrono.subsets by using the bind.data = TRUE option. In this case, the matrices need to all have the same row names and the trees all need the same labels as before:

## Making three "proximity" time slices across three trees and three bound matrices
bound_data <- chrono.subsets(matrices, trees,
                             method = "continuous",
                             model = "proximity",
                             time = 3,
                             bind.data = TRUE)
## Making three "proximity" time slices across three trees and three matrices
unbound_data <- chrono.subsets(matrices, trees,
                               method = "continuous",
                               model = "proximity",
                               time = 3,
                               bind.data = FALSE)

## Measuring disparity as the sum of variances and summarising it
summary(dispRity(bound_data, metric = c(sum, variances)))
summary(dispRity(unbound_data, metric = c(sum, variances)))

Note here that the results are again rather different: with the bound data, the slices are done across the three trees and each of their corresponding matrix (resulting in three observation) which is more accurate than the previous results from three_trees above. With the unbound data, the slices are done across the three trees and applied to the three matrices (resulting in 9 observations). As we've seen before, this is incorrect in this case since the trees don't have the same topology (so the nodes selected by a slice through the second tree are not equivalent to the nodes in the first matrix) but it can be useful if the topology is fixed to integrate both uncertainty in branch length (slicing through different trees) and uncertainty from, say, ancestral states estimations (applying the slices on different matrices).

Note that since the version 1.8 the trees and the matrices don't have to match allowing to run disparity analyses with variable matrices and trees. This can be useful when running ancestral states estimations from a tree distribution where not all trees have the same topology.

Disparity with trees: dispRitree! {#dispRitree}

Since the package's version 1.5.10, trees can be directly attached to dispRity objects. This allows any function in the package that has an input argument called "tree" to automatically intake the tree from the dispRity object. This is especially useful for disparity metrics that requires calculations based on a phylogenetic tree (e.g. ancestral.dist or projections.tree) and if phylogeny (or phylogenie*s*) are going to be an important part of your analyses.

Trees are attached to dispRity object as soon as they are called in any function of the package (e.g. as an argument in chrono.subsets or in dispRity) and are stored in my_dispRity_object$tree. You can always manually attach, detach or modify the tree parts of a dispRity object using the utility functions get.tree (to access the trees), remove.tree (to remove it) and add.tree (to... add trees!). The only requirement for this to work is that the labels in the tree must match the ones in the data. If the tree has node labels, their node labels must also match the data. Similarly if the data has entries for node labels, they must be present in the tree.

Here is a quick demo on how attaching trees to dispRity objects can work and make your life easy: for example here we will measure how the sum of branch length changes through time when time slicing through some demo data with a acctran split time slice model (see more info here).

## Loading some demo data:
## An ordinated matrix with node and tip labels
data(BeckLee_mat99)
## The corresponding tree with tip and node labels
data(BeckLee_tree)
## A list of tips ages for the fossil data
data(BeckLee_ages)

## Time slicing through the tree using the equal split algorithm
time_slices <- chrono.subsets(data   = BeckLee_mat99,
                              tree   = BeckLee_tree,
                              FADLAD = BeckLee_ages,
                              method = "continuous",
                              model  = "acctran",
                              time   = 15)

## We can visualise the resulting trait space with the phylogeny
## (using the specific argument as follows)
plot(time_slices, type = "preview",
     specific.args = list(tree = TRUE))
## Note that some nodes are never selected thus explaining the branches not reaching them.

And we can then measure disparity as the sum of the edge length at each time slice on the bootstrapped data:

## Measuring the sum of the edge length per slice
sum_edge_length <- dispRity(boot.matrix(time_slices), metric = c(sum, edge.length.tree))
## Summarising and plotting
summary(sum_edge_length)
plot(sum_edge_length)

Of course this can be done with multiple trees and be combined with an approach using multiple matrices (see here)!

Disparity of variance-covariance matrices (covar) {#covar}

Variance-covariance matrices are sometimes a useful way to summarise multidimensional data. In fact, you can express the variation in your multidimensional dataset directly in terms of how your trait covary rather than simply the positions of your elements in the trait space. Furthermore, variance-covariance matrices can be estimated from multidimensional in sometimes more useful ways that simply looking at the the data in your trait space. This can be done by describing your data as hierarchical models like generalised linear mixed effect models (glmm). For example, you might have a multidimensional dataset where your observations have a nested structure (e.g. they are part of the same phylogeny). You can then analyse this data using a glmm with something like my_data ~ observations + phylogeny + redisduals. For more info on these models start here. For more details on running these models, I suggest using the MCMCglmm package (@MCMCglmm) from @hadfield2010 (but see also @mulTree). For an example use of this code, see @guillerme2023innovation.

Creating a `dispRity` object with a `$covar` component

Once you have a trait space and variance-covariance matrices output from the MCMCglmm model, you can use the function MCMCglmm.subsets to create a "dispRity" object that contains the classic "dispRity" data (the matrix, the subsets, etc...) but also a the new $covar element:

## Loading the charadriiformes data
data(charadriiformes)

Here we using precaculated variance-covariance matrices from the charadriiformes dataset that contains a set of posteriors from a MCMCglmm model. The model here was data ~ traits + clade specific phylogenetic effect + global phylogenetic effect + residuals. We can retrieve the model information using the MCMCglmm utilities tools, namely the MCMCglmm.levels function to directly extract the terms names as used in the model and then build our "dispRity" object with the correct data, the posteriors and the correct term names:

## The term names
model_terms <- MCMCglmm.levels(charadriiformes$posteriors)[1:4]
## Note that we're ignoring the 5th term of the model that's just the normal residuals

## The dispRity object
MCMCglmm.subsets(data = charadriiformes$data,
                 posteriors = charadriiformes$posteriors,
                 group = model_terms)

As you can see this creates a normal dispRity object with the information you are now familiar with. However, we can be more fancy and provide more understandable names for the groups and provide the underlying phylogenetic structure used:

## A fancier dispRity object
my_covar <- MCMCglmm.subsets(data = charadriiformes$data,
                             posteriors = charadriiformes$posteriors,
                             group = model_terms,
                             tree = charadriiformes$tree,
                             rename.groups = c(levels(charadriiformes$data$clade), "phylogeny"))
## Note that the group names is contained in the clade column of the charadriiformes dataset as factors

Visualising covar objects

One useful thing to do with these objects is then to visualise them in 2D. Here we can use the covar.plot function (that has many different options that just plot.dispRity for plotting covar objects) to plot the trait space, the 95% confidence interval ellipses of the variance-covariance matrices and the major axes from these ellipses. See the ?covar.plot help page for all the options available:

par(mfrow = c(2,2))
## The traitspace
covar.plot(my_covar, col = c("orange", "darkgreen", "blue"), main = "Trait space")
## The traitspace's variance-covariance mean ellipses
covar.plot(my_covar, col = c("orange", "darkgreen", "blue", "grey"), main = "Mean VCV ellipses",
           points = FALSE, ellipses = mean) 
## The traitspace's variance-covariance mean ellipses
covar.plot(my_covar, col = c("orange", "darkgreen", "blue", "grey"), main = "Mean major axes",
           points = FALSE, major.axes = mean)
## A bit of everything
covar.plot(my_covar, col = c("orange", "darkgreen", "blue", "grey"), main = "Ten random VCV matrices",
           points = TRUE, major.axes = TRUE, points.cex = 1/3, n = 10, ellipses = TRUE, legend = TRUE)

Disparity analyses with a `$covar` component

You can then calculate disparity on the "dispRity" object like shown previously. For example, you can get the variances of the groups that where used in the model by using the normal dispRity function:

summary(dispRity(my_covar, metric = variances))

However this is not applied on the variance-covariance matrices from the posteriors of the MCMCglmm. To do that, you need to modify the metric to be recognised as a "covar" metric using the as.covar function. This function transforms any disparity metric (or disparity metric style function) to be applied to the $covar part of a "dispRity" object. Basically this $covar part is a list containing, for each posterior sample $VCV, the variance-covariance matrix and $loc, it's optional location in the traitspace.

## The first variance covariance matrix for the "gulls" group
my_covar$covar[["gulls"]][[1]]

And this is how as.covar modifies the disparity metric:

## Using the variances function on a VCV matrix
variances(my_covar$covar[["gulls"]][[1]]$VCV)

## The same but using it as a covar metric
as.covar(variances)(my_covar$covar[["gulls"]][[1]])

## The same but applied to the dispRity function
summary(dispRity(my_covar, metric = as.covar(variances)))

Disparity and distances

There are two ways to use distances in dispRity, either with your input data being directly a distance matrix or with your disparity metric involving some kind of distance calculations.

Disparity data is a distance

If your disparity data is a distance matrix, you can use the option dist.data = TRUE in dispRity to make sure that all the operations done on your data take into account the fact that your disparity data has distance properties. For example, if you bootstrap the data, this will automatically bootstrap both rows AND columns (i.e. so that the bootstrapped matrices are still distances). This also improves speed on some calculations if you use disparity metrics directly implemented in the package by avoiding recalculating distances (the full list can be seen in ?dispRity.metric - they are usually the metrics with dist in their name).

Subsets

By default, the dispRity package does not treat any matrix as a distance matrix. It will however try to guess whether your input data is a distance matrix or not. This means that if you input a distance matrix, you might get a warning letting you know the input matrix might not be treated correctly (e.g. when bootstrapping or subsetting). For the functions dispRity, custom.subsets and chrono.subsets you can simply toggle the option dist.data = TRUE to make sure you treat your input data as a distance matrix throughout your analysis.

## Creating a distance matrix
distance_data <- as.matrix(dist(BeckLee_mat50))

## Measuring the diagonal of the distance matrix
dispRity(distance_data, metric = diag, dist.data = TRUE)

If you use a pipeline of any of these functions, you only need to specify it once and the data will be treated as a distance matrix throughout.

## Creating a distance matrix
distance_data <- as.matrix(dist(BeckLee_mat50))

## Creating two subsets specifying that the data is a distance matrix
subsets <- custom.subsets(distance_data, group = list(c(1:5), c(6:10)), dist.data = TRUE)
## Measuring disparity treating the data as distance matrices
dispRity(subsets, metric = diag)

## Measuring disparity treating the data as a normal matrix (toggling the option to FALSE)
dispRity(subsets, metric = diag, dist.data = FALSE)
## Note that a warning appears but the function still runs

Bootstrapping

The function boot.matrix also can deal with distance matrices by bootstrapping both rows and columns in a linked way (e.g. if a bootstrap pseudo-replicate draws the values 1, 2, and 5, it will select both columns 1, 2, and 5 and rows 1, 2, and 5 - keeping the distance structure of the data). You can do that by using the boot.by = "dist" function that will bootstrap the data in a distance matrix fashion:

## Measuring the diagonal of a bootstrapped matrix
boot.matrix(distance_data, boot.by = "dist")

Similarly to the dispRity, custom.subsets and chrono.subsets function above, the option to treat the input data as a distance matrix is recorded and recycled so there is no need to specify it each time.

Disparity metric is a distance

On the other hand if your data is not a distance matrix but you are using a metric that uses some kind of distance calculations, you can use the option dist.helper to greatly speed up calculations. dist.helper can be either a pre-calculated distance matrix (or a list of distance matrices) or, better yet, a function to calculate distance matrices, like stats::dist or vegan::vegdist. This option directly stores the distance matrix separately in the RAM and allows the disparity metric to directly access it at every disparity calculation iteration, making it much faster. Note that if you provide a function for dist.helper, you can also provide any un-ambiguous optional argument to that function, for example method = "euclidean".

If you use a disparity metric implemented in dispRity, the dist.helper option is correctly loaded onto the RAM regardless of the argument you provide (a matrix, a list of matrix or any function to calculate a distance matrix). On the other hand, if you use your own function for the disparity metric, make sure that dist.helper exactly matches the internal distance calculation function. For example if you use the already implemented pairwise.dist metric all the following options will be using dist.helper optimally:

## Using the dist function from stats (specifying it comes from stats)
dispRity(my_data, metric = pairwise.dist, dist.helper = stats::dist)

## Using the dist function from vegdist function (without specifying its origin)
dispRity(my_data, metric = pairwise.dist, dist.helper = vegdist)

## Using some pre-calculated distance with a generic function
my_distance_matrix <- dist(my_distance_data)
dispRity(my_data, metric = pairwise.dist, dist.helper = my_distance_matrix)

## Using some pre-calculated distance with a user function defined elsewhere
my_distance_matrix <- my.personalised.function(my_distance_data)
dispRity(my_data, metric = pairwise.dist, dist.helper = my_distance_matrix)

However, if you use a homemade metric for calculating distances like this:

## a personalised distance function
my.sum.of.dist <- function(matrix) {
  return(sum(dist(matrix)))
}

The dist.helper will only work if you specify the function using the same syntax as in the user function:

## The following uses the helper correctly (as in saves a lot of calculation time)
dispRity(my_data, metric = my.sum.of.dist, dist.helper = dist)

## These ones however, work but don't use the dist.helper (don't save time)
## The dist.helper is not a function
dispRity(my_data, metric = my.sum.of.dist, dist.helper = dist(my_data))
## The dist.helper is not the correct function (should be dist)
dispRity(my_data, metric = my.sum.of.dist, dist.helper = vegdist)
## The dist.helper is not the correct function (should be just dist)
dispRity(my_data, metric = my.sum.of.dist, dist.helper = stats::dist)

TGuillerme/dispRity documentation built on Dec. 21, 2024, 4:05 a.m.

rdrr.io home R language documentation Run R code online

CRAN packages Bioconductor packages R-Forge packages GitHub packages

Note that we can't provide technical support on individual packages. You should contact the package authors for that.

TGuillerme/dispRity Measuring Disparity