shiftingspace.Rmd

---
title: "Shifting spaces: which disparity or dissimilarity measurement best summarise occupancy in multidimensional spaces?"
author: "Thomas Guillerme, Mark N. Puttick, Ariel E. Marcy, Vera Weisbecker"
bibliography: references.bib
csl: mee.csl
date: "`r Sys.Date()`"
output:
  html_document:
    fig_width: 8
    fig_height: 8
    keep_tex: true
  self_contained: true
---


<!-- Keywords: disparity, dissimilarity, ecology, evolution, multidimensionality, statistics -->


<!-- When toggling to html, toggle snippet `compilation_html` below to eval = TRUE -->

<!-- To build the line numbered version (for submission) use the following yaml header snippet
output:
  pdf_document:
    fig_width: 8
    fig_height: 8
    keep_tex: true
  self_contained: false

Instead of

output:
  html_document:
    fig_width: 8
    fig_height: 8
---

Switch the text highlighting by replacing

<@@@font style="color:black">
to:
\@@@textcolor{blue}{

and 

<@@@/font>
to:
}@@@%col_tag\n

And then edit the }%col_tags in latex

----
and edit the save .tex file to add the following

\usepackage{lineno} %in the package header

\linespread{1.6} %in the document header

\modulolinenumbers[1] % just after the \begin{document} tag
\linenumbers
 -->
<!-- 

# Response to reviewers

We a really grateful for the reviewers refreshing, positive and constructive comments and we've addressed them all as detailed below.

## Reviewer 1 (Stefano Mammola)

In a conceptually similar study in the context of species distribution modelling (DOI: 10.1111/2041-210X.12397), Qiao and colleagues termed this general idea the “no silver bullet” paradigm (making a metaphor with the mythology of lycanthropes). This is certainly a very important point also in trait-based science. 

> This is indeed a really similar idea with similar conclusions. We've now added a reference to niche modeling in the introduction:

"This can also be extended to more complex ecological concepts such as niche modelling [@qiao2015]." l.@@@

1) While reading the text, it was not entirely clear to me what exactly the “moms” tool is. Is this an R package, or just a set of functions? Would be helpful if you could briefly specify.

> We've added an extra section explaining what the `moms` tool is. We've also added a sampling (see reviewer 2) and simulation feature allowing users to replicate the results of this paper with several clicks or exploring the properties of their own metrics:

"Therefore, we propose the [`moms`](https://tguillerme.shinyapps.io/moms/) shiny app to allow workers to help them choose their set of space occupancy measurements (and test the caveats mentioned above).
`moms` is an online graphical user interface to help analyse multidimensional data.
It allows users to upload their dataset of interest (or simulate one with specific parameters) and measure space occupancy using a variety of implemented measures (namely, but not only, the ones used in this study).
Furthermore, the package allows simulation of shifts in trait space occupancy as also presented in this paper to test whether some measures capture specific changes in space.
However, `moms` is not a tool for analysing multidimensional data _per se_ but rather for helping workers to chose the space occupancy measure most appropriated to their data and question.
To run multidimensional analysis, we suggest using dedicated `R` packages (such as - but not limited to: @oksanen2007vegan, @momocs, @bat2015, @disprity)." l.@@@

2)  I suspect that your proposed classification may somewhat overlap with the general classification scheme of diversity indexes by Pavoine and Bonsall (2011; Biological review), later expanded by Tucker et al. (2017; Biological Reviews) in the context of phylogenetic metrics (but it’s actually the same). In a nutshell, they grouped metrics/indexes in three main category of ‘Richness’, ‘Divergence’ and ‘Regularity’ components. In their view, these three concepts should capture the primary mathematical operation inherent to each metric, namely:
i) the ‘richness’ dimension encompasses indexes representing the sum of difference among taxa (sum);
ii) the ‘divergence’ dimension encompasses indexes representing functional (or phylogenetic) dissimilarity, reflecting the average difference among taxa (mean); and
iii) the  ‘regularity’ dimension encompasses indexes representing functional (or phylogenetic) variance, reflecting how regular the difference among taxa in an assemblages are (variance).
I see some overlap with your method, but also differences. It would be worth briefly acknowledging this and perhaps exploring the main conceptual differences and why your classification is an advance.

>We were not aware of these references and we are grateful for the reviewer pointing them out. There is so much overlapping work on ecology and evolution but so few visibility between them. Adding these references and concepts is definitely improving the main idea of this paper!

Note that this classification bears some similarities with @tucker2017 classifying phylogenetic diversity measurements into richness, divergence and regularity categories.
However, while @tucker2017 based their classification on the mathematical operation inherent to each metrics (the sum for richness, the mean for divergence and the variance for regularity), our three broad classifications are based on their geometric properties regardless of the formula of each metric (e.g. the size of a space can be calculated using a sum, mean or/and variance). l.@@@

3) The approach you choose for simulating changes in the multidimensional space (Figure 2) is  very interesting. While I was thinking about it, I indulged with an idea (but please ignore if not appropriate):

I thought it would be useful to try linking the different changes in the multidimensional space to actual example of biological processes (e.g., giving few examples in the description in the method or in the discussion). This, because similar changes in the trait space may occur in the real world due to different processes and, by providing examples grounded in the real world you may enhance the appealing of the ms to a broader audience.
For example, from the perspective of a conservation biologist, you may argue that ‘limit’ and ‘position’ change in the trait space may occur when there is the destruction of habitats with narrow environmental conditions that filters for few species possessing specific traits. I’m primarily a cave biologist: this mental bias made me think that, for example, if you open a quarry that destroys a cave (not an unfrequent event! DOI: 10.1093/biosci/biz064), you may end up removing from the trait space all points from a specific position, namely those traits clearly selected by the permanent darkness of the subterranean world. If trait 1 is a gradient of eye size and trait 2 a gradient of body pigment, you would remove trait contribution of specialized invertebrates with no eyes and no pigment. The ‘density’ change may occur, e.g., due to global wildlife trade, when there is the exploitation of only species with specific traits within the total functional tree of life (see DOI: 10.1126/science.aav5327). For example, the fishing industry exploits, within any given species, primarily fishes of larger size, thus selecting for specific densities in the size/weight trait space composed of the multiple species. And so on (these are just random examples to illustrate my point).

> We've added this reviewer's example (although not all to also balance with examples from macroevolution):
    
This type of change could be due to habitat destruction [e.g. @mammola2019b] or to mass extinctions [e.g. @wright2017]. l.@@@

This type of change could be due to accelerated rates of evolution [@close2015] or to differences in modes of life in macroevolution [e.g. @healy2019]. l.@@@

This type of change could be due changes in evolutionary trajectories [@endler2005] or to differences in ecosystem compositions [e.g. @jones2015]. l.@@@


### Minor Line comments

-L34: “studied” instead of “study”

> We fixed this typo.

L35-37: I think merging together these two sentences would put more emphasis on your idea... E.g. “While different fields use a different set of terms for such approaches (Table 1), they actually focus on the same mathematical objects: matrices with columns representing an original or transformed trait value and rows representing observations (taxon, field site, etc.; Guillerme 2018). 

> Merged these sentences.

-L63: “are are”

> We fixed this typo.

-Table 1: 
‘functional space’ instead of ‘function-space’? 
Morphospace I think is also widely used in ecology. 
Please add “etc.” at the end of the text in the cell “statistic” intersected “ecology”. Wouldn’t the “hypervolume” mentioned in the same cell more fit in the “Matrix (n x d)” row?

> We fixed these cells in the table. However, we did not add "hypervolume" to the "Mathematics" column. We use the term "statistic" here as a statistic (i.e. a measure) which does include the hypervolume as used in ecology. We've specified this.

-330–332: True, but there are at least two recent examples of attempt in this sense, the very recent package TPD (DOI: 10.1002/ecy.2876) and new functions in BAT (doi: 10.1101/2020.01.25.919373; shameless self-promotion and still a preprint). Both approaches covers metric in the Richness, Divergence, and Regularity domains (which should be a possible equivalent to your Size, Density, and Position view).

> We've added these references on line @@@.

-Figure 2: why the dots are in black and white in the first inset?

> The black dots are the 50% randomly removed points. We've specified this in the caption.


## Reviewer: 2 (Neil Brocklehurst) 

A general suggestion to start with. The authors state in their conclusions, very diplomatically, that "...no measure is better than the next one...". While I agree with the sentiment in the context of their conclusions (that different metrics show different things, and you should choose a variety of metrics relevant to the issue under study), would it not be accurate to say that your results indicate some metrics are worse than others? Three leapt out at me in table 5: Sum of Ranges, Minimum spanning tree average distances, and to a lesser extent average nearest neighbour distances,  seem to produce changes in value following random removals comparable to those seen in not one, but all three varieties of separation (limit, density and position) . Would it be fair to say that these metrics are, shall we say, problematic?

> We agree with the reviewer and have toned down or "diplomacy" level. We now also state concern about certain types of metrics that are highly sensitive to outliers or to the number of dimensions:

We insist that although no measure is objectively better than the next one, some can be more problematic than other in specific contexts. 
For example, the results for the Sum of Ranges, Minimum spanning tree average distances, and to a lesser extent average nearest neighbour distances produced results in the reduced space often similar to the randomly reduced spaces (Table 5).
This does not make them "bad" measures but rather heavily context dependent.
Regardless, we believe that workers should identify the most appropriate measures based on their trait space properties as well as their specific biological question.
We believe this could be fostered by following these several suggestions: l.@@@


Following on from this: its worth testing, or at the very least acknowledging, that some metrics probably respond to incomplete sampling worse than others? Or are certain targeted separations (by limit, density or position) harder to detect with incomplete sampling? A test of this would be reasonably simple; after carrying out the separations by limit, density or position, randomly delete different numbers of data points to see how this affects the results. I realise this would be a large number of extra analyses to ask for and may be considered beyond the immediate scope of the paper, so I’ll leave the choice of whether or not to do the actual analyses up to the authors, but I think it should at least be discussed in the text. 

> We agree again with this reviewer's point and find this a really interesting point. We've updated the results in the supplementary material to highlight what happens when removing 80% and 20% of the elements (rather than 50%).
We've also added a sampling option to the latest version of the `moms` app allowing to test the effect of sampling on the metric of choice and, finally, we've added the following mentions to sampling issues in the caveat section:

Furthermore, we did not take into account the effect of sampling on space occupancy measurements (but see additional results with 80% and 20% space reduction in the supplementary materials 4).
In fact, sampling has been previously shown to have an effect on measurements depending on range or volumes (e.g. the sum of ranges or the hypervolume @ciampaglio2001).
This effect is especially expected to be acerbated in macroevolutionary studies when using the fossil record [@brocklehurst2013] but can be tackled using rarefaction and bootstrapping techniques [@disprity]. l.@@@

Some more specific comments:
Fig 2, panels E and F, and similar in table 5: the dots and lines representing the changes in occupancy for blue, orange and random groups. I assume the dot is the median and the solid and dashed lines are quartiles and ranges, but specify in caption, please

> We already specified that "the dots represent the median space occupancy values across all simulations for each scenario of trait space change (Table 2), the solid and dashed line respectively the 50% and 95% confidence intervals" (Caption Figure 2). We've now also added this to the Table 5 caption.

Table 5 headers are inconsistent with the terms used in the text and figure 5. In the text and figure 5 you describe Limit and Density parameters for changes, but the headers talk about size change and arrangement change. Maybe this will create confusion (particularly the ‘arrangement’ header, which might not naturally be associated in peoples minds with density)?

> Nice spotting, this comes from a classic version mashup. We've fixed this in the table header and double checked throughout the manuscript.

And finally some typos and wording comments:

First sentence of the abstrach: “Multidimensional analysis of traits are now a common in ecology and evolution…”; should be either “…now common in…” or “…now a common method in…”

> We remove "a".

Abstract, line 8: “…subspace of this space…”; maybe say subset instead of subspace? The two ‘spaces’ next two each other makes slightly awkward reading

> We've changed "subspace" to "subset".

Intro line 48: “ciampaglio2001”; Space missing between name and year, name needs capital letter
Methods, page 8, line 153: “hopkins2017”; Space missing between name and year, name needs capital letter

> We fixed the references.

Page 20, line 317: “spacesmight” missing space between words

> We've added a space.
 -->


```{r header, echo = FALSE, results = 'hide', message = FALSE, warning = FALSE}
## Repeatability note:
## This whole paper is entirely reproducible and compiles as a single document. The data, figures
## and tables are all generated through the code snippets in this file. Note that the first time
## that this paper is compiled, it will generate the data which will take substantial time
## as the shiftingspace_supplementary.Rmd file completion takes ~10 minutes. Subsequent compilations
## will be much faster!

## Loading the packages
if(!require(devtools)) install.packages("devtools")
if(!require(knitr)) install.packages("knitr"); library(knitr)
if(!require(rmarkdown)) install.packages("rmarkdown"); library(rmarkdown)
if(!require(ape)) install.packages("ape"); library(ape)
if(!require(dispRity)) install.packages("dispRity"); library(dispRity)
if(packageVersion("dispRity") < "1.2.4") {
    ## dispRity must be above v1.2.3
    devtools::install_github("TGuillerme/dispRity"); library(dispRity)
}
if(!require(moms)) devtools::install_github("TGuillerme/moms"); library(moms)

## Setting the datapath
DATA_PATH <- "../data/processed/"

## Create data directory
if(!dir.exists(DATA_PATH)) {
    dir.create(path = DATA_PATH, showWarnings = FALSE)
}

## Setting the default parameters for the space plots
defaults <- list(pch = 20,
                 xlim = c(-3, 3),
                 ylim = c(-3, 3),
                 col1 = "grey",
                 col2 = "black",
                 xlab = "Trait",
                 ylab = "Trait",
                 cex = 1)
## Generating the default palette
default.palette <- function(n) {
    hues <- seq(15, 375, length = n + 1)
    grDevices::hcl(h <- hues, l = 65, c = 100)[1:n]
}

## Checking whether the data exists
if(!all(c("remove_05.Rda") %in% list.files(DATA_PATH))) {
    ## Run The supplementary material (simulation) ~ 30 minutes
    rmarkdown::render("shiftingspace_supplementary_simulation.Rmd", "html_document")    
}
if(!all(c("empirical_results.Rda") %in% list.files(DATA_PATH))) {
    ## Run The supplementary material ~ 15 minutes
    rmarkdown::render("shiftingspace_supplementary_empirical.Rmd", "html_document")
}
```

```{r compilation_html, echo = FALSE, eval = TRUE}
## Changing defaults
body(plot.id)[[2]] <- substitute(type <- ".png")
```

# Abstract


Multidimensional analysis of traits are now common in ecology and evolution and are based on trait spaces in which each dimension summarises the observed trait combination (a morphospace or an ecospace).
Observations of interest will typically occupy a <font style="color:black">subset</font> of this space, and researchers will calculate one or more measures to quantify how organisms inhabit that space.
In macroevolution and ecology these measures <font style="color:black">called disparity or dissimilarity metrics and are</font> generalised as space occupancy measures.
Researchers use these measures to investigate how space occupancy changes through time, in relation to other groups of organisms, and in response to global environmental changes.
However, the mathematical and biological meaning of most space occupancy measures is vague with the majority of widely-used measures lacking formal description.

Here we propose a broad classification of space occupancy measures into three categories that capture changes in size, density, or position.
We study the behaviour of 25 measures to changes in trait space size, density and position on simulated and empirical datasets.
We find that no measure describes all of trait space aspects but that some are better at capturing certain aspects. 
Our results confirm the three broad categories (size, density and position) and allow us to relate changes in any of these categories to biological phenomena.

Because the choice of space occupancy measures is specific to the data and question, we introduced [`moms`](https://tguillerme.shinyapps.io/moms/), a tool to both visualise and capture changes in space occupancy for any measurement.
[`moms`](https://tguillerme.shinyapps.io/moms/) is designed to help workers choose the right space occupancy measures, given the properties of their trait space and their biological question.
By providing guidelines and common vocabulary for space occupancy analysis, we hope to help bridging the gap in multidimensional research between ecology and evolution.

# Introduction

Groups of species and environments share specific, recognisable, correlated characteristics: guilds or biomes with shared phenotypic, physiological, phylogenetic or behavioural traits.
Organisms or environments should therefore be studied as a set of traits rather than some specific traits in isolation [@donohue2013; @hopkins2017].
Biologists increasingly been using ordination techniques [see @legendre2012 for a summary] to create multidimensional trait spaces to either explore properties of data or test hypotheses [e.g. @oksanen2007vegan; @blonder2018; @disprity].
For example, in palaeobiology, @wright2017 used trait spaces to study how groups of species' characteristics change through time; in ecology, @jones2015 <font style="color:black"> studied </font> evidence of competition by looking at trait overlap between two populations.
<font style="color:black"> While different fields use a different set of terms for such approaches (Table 1), they actually focus on the same mathematical objects: matrices with columns representing an original or transformed trait value and rows representing observations [taxon, field site, etc.; @disprity]. </font>

Mathematics | Ecology | Macroevolution | This paper
--------------------|--------------------|--------------------|------------
Matrix ($n \times d$) with a structural relation between rows and columns | <font style="color:black"> Functional space, morphospace </font>, etc. | Morphospace, traitspace, etc. | trait space
Rows (*n*) | Taxa, field sites, environments, etc. | Taxa, specimen, populations, etc. | observations
Columns (*d*) | Traits, Ordination scores, distances, etc. | Traits, ordination scores, distances, etc. | dimensions
Matrix subset ($m \times d$; $m \leq n$) | Treatments, phylogenetic group (clade), etc. | Clades, geological stratum, etc. | group
Statistic <font style="color:black"> (i.e. a measure) </font>  | Dissimilarity index or metric, hypervolume, functional diversity, <font style="color:black"> etc. </font> | Disparity metric or index | space occupancy measure
Multidimensional analysis | Dissimilarity analysis, trait analysis, etc. | Disparity analysis, disparity-through-time, etc. | multidimensional analysis
Table 1: <font style="color:black">Different terms are used for equivalent measures in</font> mathematics, ecology and macroevolution.

Ecologists and evolutionary biologists often use trait spaces with respect to the same fundamental questions:
are groups occupying the same amount of trait space?
Do some groups contain more species than others in the same amount of trait space?
Are some specific factors correlated with different patterns of trait space occupancy?
Because of the multidimensional nature of these trait spaces, it is often not possible to study them using bi- or tri-variate techniques [@diaz2016; @hopkins2017; @mammola2019].
Studying the occupancy of trait spaces is done using disparity indices in macroevolution [@wills2001; @hopkins2017; @disprity] or comparing hypervolumes in ecology [@donohue2013; @diaz2016; @blonder2018; @mammola2019].
Despite the commonalities between the measures used in ecology and evolution (which are often metric but don't necessarily need to be), surprisingly little work has been published on their behaviour [but see @ciampaglio2001; @villéger2008; @mammola2019].

Different occupancy measures capture different aspects of  trait space [@ciampaglio2001; @villéger2008; @mammola2019].
<font style="color:black">This</font> may be widely-known, but to our knowledge it is infrequently mentioned in peer-reviewed papers.
First, space occupancy measures are often named as the biological aspect they are describing ("disparity", "functional diversity") rather than what they are measuring (e.g. the product of ranges), which obscures the differences and similarities between studies.
Second, in many studies in ecology and evolution, authors have focused on measuring the size of the trait space [e.g. ellipsoid volume @donohue2013; hypervolume @diaz2016; product of variance @wright2017; Procrustes variance @marcy2016].
However, the size of the trait space only represents one aspects of occupancy, disregarding <font style="color:black">other measures</font> such as the density [@geiger2008] or position [@wills2001; @ciampaglio2001].
For example, if two groups have the same size, this can support certain biological conclusions.
Yet, an alternative aspect of space occupancy may indicate that the groups' position are different, leading to a different biological conclusion (e.g. the groups are equally diverse but occupy different niches).
Using measures that only capture one aspect of the trait space may restrain the potential of multidimensional analysis [@villéger2008].

Here we propose a broad classification of space occupancy measures as used across ecology and evolution and study their power to detect changes in trait space occupancy in simulated and empirical data.
<font style="color:black">Note this does not account whether or not it is possible for a space to be occupied (e.g., some spaces may represent biologically impossible shapes); this, however, may be important in some cases, such as testing whether a region is infinite or not.</font>
We provide an assessment of each broad type of space occupancy measures along with a unified terminology to foster communication between ecology and evolution. 
Unsurprisingly, we found no one measure describes all changes <font style="color:black">in space</font> and that the results from each measures are dependent on the characteristics of the space and the hypotheses.

<font style="color:black">There can be an infinite number of measures and that it is thus impossible to propose a comprehensive analysis for all the measures properties respective to how they measure changes in trait space.
We therefore propose [`moms`](https://tguillerme.shinyapps.io/moms/), a tool for researchers to design, experiment and visualise their own space occupancy measure tailored for their project.
The tool will help researchers understand the "null" behaviour of the measures of interest.</font>

## Space occupancy measures

In this paper, we define trait spaces as any matrix where rows are observations and columns are traits, where both observations and traits are structurally related (e.g. there is a phylogenetic relation between observations - and traits, etc.).
These traits can widely vary in number and types: they could be coded as discrete [e.g. presence or absence of a bone; @beck2014; @wright2017], continuous measurements [e.g. leaf area; @diaz2016] or more sophisticated measures [e.g. landmark position; @marcy2016; Fourier ellipses; @momocs].
Traits can also be measured by using relative observations [e.g. community compositions; @jones2015] or distance between observations [e.g. @close2015].
However, regardless of the methodology used to build a trait space, three broad occupancy measures can be used: the **size** which approximates the amount of space occupied, the **density** which approximates the distribution in space and the **position** which approximates the location in space [Fig. 1; @villéger2008].
Of course any combination of these three aspects is always possible.

```{r fig_measures_types, echo = FALSE, fig.height = 3, fig.width = 9, results = 'hide', fig.cap = paste("different type of information captured by space occupancy measures: (A) size, (B) density and (C) position.") }

set.seed(11)
## Plot space function (utility shortcut)
## The elements
elements <- 10

## Trait space
trait_space <- space.maker(elements, 2, distribution = rnorm)

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

## The size
plot(trait_space, pch = defaults$pch, main = "A - Size\n(e.g. Sum of ranges)", xlab = "Trait 1", ylab = "Trait 2")
## The range lines
lines(x = range(trait_space[,1]), y = rep(mean(trait_space[,2]), 2), col = defaults$col1)
lines(x = rep(mean(trait_space[,1]), 2), y = range(trait_space[,2]), col = defaults$col1)

## The density
plot(trait_space, pch = defaults$pch, main = "B - Density\n(e.g. pairwise distances)", xlab = "Trait 1", ylab = "")
## Plotting the pairwise lines
pair.line <- function(line, trait_space, defaults) {
    lines(x = trait_space[line,1], y = trait_space[line,2], col = defaults$col1)
}
apply(combn(1:elements, 2), 2, pair.line, trait_space, defaults)
points(trait_space, pch = defaults$pch)

## The position
plot(trait_space, pch = defaults$pch, main = "C - Position\n(e.g. distance from centre)", xlab = "Trait 1", ylab = "")
## Plotting the center of the space
arrows(x0 = 0, y0 = 0, x1 = trait_space[,1], y1 = trait_space[,2], code = 0, col = defaults$col1)
points(trait_space, pch = defaults$pch)
points(0,0, pch = 13, cex = 2)

par(op)
```

#### 1. Size

Size captures the spread of a group in the trait space.
They can be interpreted as the amount of the trait space that is occupied by observations.
Typically, larger values for such measures indicate the presence of more extreme trait combinations.
For example, if group A is bigger than B, the observations in A achieve more extreme trait combinations than in B.
This type of measure is widely used in both ecology [e.g. the hypervolume; @blonder2018] and in evolution [e.g. the sum or product of ranges or variances; @wills2001].

Although size measures are suitable indicators of a group's trait space occupancy, they are limited to comparing the range of trait-combinations between groups.
Size measures do not take into account the distribution of the observations within a group and can often be insensitive to unoccupied "holes" in the trait space (overstimating the size; @blonder2018).
They can make it difficult to determine whether all the observations are on the edge of the group's distribution or whether the size is simply driven by outliers.


#### 2. Density

Density gives an indication of the quantity of observations in the trait space.
They can be interpreted as the distribution of the observations *within* a group in the trait space.
Groups with higher density contain more observations (i.e. more observations per approximation of size) that will tend to be more similar to each other.
For example, if group A is greater is size than group B and both have the same density (observations are equally distant within each group), similar mechanisms could be driving both groups' trait space occupancy. 
Indeed, this pattern could suggest that A is older and has had more time to achieve more extreme trait combinations under essentially the same process as younger, smaller group B [@endler2005].
Note that density based measures can be sensitive to sampling.
Density measures are less common compared to size measures, but they are still used in both ecology [e.g. the minimum spanning tree length; @oksanen2007vegan] and evolution [e.g. the average pairwise distance; @geiger2008]. 

#### 3. Position

Position captures where a group lies in trait space.
They can be interpreted as where a group lies in the trait space either relative to the space itself or relative to another group.
For example, if group A has a different position than group B, A will have a different trait-combination than in B.

Position measures may be harder to interpret in multidimensional spaces.
In a 2D space, two groups can be equally distant from a fixed point but in different parts of the space (left, right, up, or down - with the amount of parts of space increasing with dimensions).
However, when thinking about unidimensional data, this measure is obvious: two groups A or B could have the same variance (size) with the same number of observations (density) but could have two different means and thus be in different positions.
These measures are used in ecology to compare the position of two groups relative to each other [@mammola2019].

<font style="color:black"> Note that this classification into size, density and position bears some similarities with @tucker2017 classifying phylogenetic diversity measurements into richness, divergence and regularity categories.
However, while @tucker2017 based their classification on the mathematical operation inherent to each metrics (the sum for richness, the mean for divergence and the variance for regularity), our three broad classifications are based on their geometric properties regardless of the formula of each metric (e.g. the size of a space can be calculated using a sum, mean or/and variance). </font>

## No measure to rule them all: benefits of considering multiple measures

The use of multiple measurements to assess trait space occupancy provides a more detailed characterisation of occupancy changes.
If the question is to look at how space occupancy changes in response to mass extinction, using a single space occupancy measure can miss part of the picture: a change in size could be decoupled from a change in position or density in trait space.
For example, the Cretaceous-Paleogene extinction (66 million years ago) shows an increase in size of the mammalian trait space [adaptive radiation; @halliday2015] but more specific questions can be answered by looking at other aspects of trait space occupancy: does the radiation expand on previously existing morphologies [elaboration, increase in density; @endler2005] or does it explore new regions of the trait space [innovation, change in position; @endler2005]?
Similarly, in ecology, if two groups have the same trait space size, <font style="color:black">the differences in density within these two groups is potentially illuminating:</font> different selection pressure can lead to different density within equally sized groups. 
<font style="color:black">This can also be extended to more complex ecological concepts such as niche modelling [@qiao2015].</font>

Here, we provide the first interdisciplinary review of 25 space occupancy measures that uses the broad classification of measures into size, density and position to capture pattern changes in trait space.
We assess the behaviour of measures using simulations and six interdisciplinary empirical datasets covering a wide range of potential data types and biological questions.
We also introduce a tool for measuring occupancy in multidimensional space ([`moms`](https://tguillerme.shinyapps.io/moms/)), which is a user-friendly, open-source, graphical interface to allow the tailored testing of measurement behaviour for any use case.
[`moms`](https://tguillerme.shinyapps.io/moms/) will allow workers to comprehensively assess the properties of their trait space and the measures associated with their specific biological question.






<!--  -->
<!--  -->
<!--  -->
<!-- METHODS SECTION -->
<!--  -->
<!--  -->
<!--  -->








# Methods

```{r metrics_list, echo = FALSE}
## Loading the list of measures
source("list.of.metrics.R")
```

We tested how 25 space occupancy measures relate to each other, are affected by modifications of traits space and affect group comparisons in empirical data:

1.  We simulated 13 different spaces with different sets of parameters;
2.  We transformed these spaces by removing 50% of the observations following four different scenarios corresponding to different empirical scenarios: randomly, by <font style="color:black">size</font> (e.g. expansion or reduction of niches), by density (e.g. different degrees of competition within a guild), and by position (e.g. ecological niche shift).
3.  We measured occupancy on the resulting transformed spaces using `r metric.names()` different space occupancy measures;
4.  We applied the same space occupancy measures to six empirical datasets (covering a range of disciplines and a range of dataset properties).

Note that the paper contains the results for only `r metric.names()` measures which were selected as representative of common measures covering the size, density and position trait space aspects.
The results for an additional 17 measures is available in the supplementary material 4.

## Generating spaces

We generated trait spaces using the following combinations of size, distributions, variance and correlation:

| space name | size  | distribution(s)                | dimensions variance | correlation |
|------------|-------|--------------------------------|---------------------|-------------|
| 3D uniform   |200*3  | Uniform (min = -0.5, max = 0.5)| Equal            | None        |
| 15D uniform  |200*15 | Uniform                        | Equal            | None        |
| 50D uniform  |200*50 | Uniform                        | Equal            | None        |
| 150D uniform |200*150| Uniform                        | Equal            | None        |
| 50D uniform correlated |200*50 | Uniform                        | Equal            | Random (between 0.1 and 0.9) |
| 3D normal    |200*3  | Normal  (mean = 0, sd = 1)     | Equal            | None        |
| 15D normal   |200*15 | Normal                         | Equal            | None        |
| 50D normal   |200*50 | Normal                         | Equal            | None        |
| 150D normal  |200*150| Normal                         | Equal            | None        |
| 50D normal correlated  |200*50 | Normal                         | Equal            | Random (between 0.1 and 0.9) |
| 50D with random distributions     |200*50 | Normal, Uniform, Lognormal (meanlog = 0, sdlog = 1)| Equal            | None        |
| 50D PCA-like   |200*50 | Normal                         | Multiplicative           | None        |
| 50D PCO-like   |200*50 | Normal                         | Additive              | None        |
Table 2: different simulated space distribution.
 *Name* of the simulated space; *dimensions* of the matrix (row\*columns); *distribution(s)* of the data on each dimensions (for the 'Random', dimensions are randomly chosen between Normal, Uniform or Lognormal); *dimension variance*: distribution of the variance between dimensions (when equal, the dimensions have the same variance); *correlation* between dimensions. 

The differences in trait space sizes (200 elemeents for 3, 15, 50 or 150 dimensions) reflects the range found in literature [e.g. @hopkins2017; @mammola2019].
We used a range of distributions (uniform, normal or a random combination of uniform, normal and lognormal) to test the effect of observation distributions on the measurements.
We used different levels of variance for each dimensions in the spaces by making the variance on each dimension either equal ($\sigma_{D1} \simeq \sigma_{D2} \simeq \sigma_{Di}$) or decreasing ($\sigma_{D1} < \sigma_{D2} < \sigma_{Di}$) with the decreasing factor being either multiplicative (using the cumulative product of the inverse of the number of dimensions: $\prod_i^d(1/d)$) or additive ($\sum_i^d(1/d)$).
Both reductions of variance are used to illustrate the properties of ordinations where the variance decreases per dimensions [lognormal in principal components analysis - PCA; e.g. @marcy2016; @healy2019; and normal win Multidimensional Scaling - MDS, PCO or PCoA; e.g. @close2015; @wright2017].
Finally, we added a correlation parameter to illustrate the effect of co-linearity between traits (especially in non-ordinated trait spaces).
We repeated the simulation of each trait space 20 times (resulting in 260 spaces).

## Spatial occupancy measures

We then  calculated `r metric.names()` different measures  on the resulting transformed spaces, including a new one, the average displacement, which we expect to be influenced by changes in trait space position.

| Name | Definition | Captures | Source | Notes  |
|--------------------|------------------------------|----------|----------|------------------------------|
| Average Euclidean distance from centroid | $\frac{\sqrt{\sum_{i}^{n}{({k}_{n}-Centroid_{k})^2}}}{d}$ | Size | @laliberté2010 | the functional dispersion (FDis - without abundance) |
| Sum of variances | $\sum_{i}^{d}{\sigma^{2}{k_i}}$ | Size | @foote1992 | common measure used in palaeobiology [@ciampaglio2001; @wills2001] |
| Sum of ranges | $\sum_{i}^{d}{\|\text{max}(d_{i})-\text{min}(d_{i})\|}$ | Size | @foote1992 | more sensitive to outliers than the sum of variances |
| Ellipsoid volume | $\frac{\pi^{d/2}}{\Gamma(\frac{d}{2}+1)}\displaystyle\prod_{i}^{d} (\lambda_{i}^{0.5})$ | Size | @donohue2013 | less sensitive to outliers than the convex hull hypervolume [@diaz2016; @blonder2018] |
| Minimum spanning tree average distance | $\frac{\sum(\text{branch length})}{n}$ | Density | @sedgewick1990 | similar to the unscaled functional evenness [@villéger2008] |
| Minimum spanning tree distances evenness | $\frac{\sum\text{min}\left(\frac{\text{branch length}}{\sum\text{branch length}}\right)-\frac{1}{n-1}}{1-\frac{1}{n-1}}$ | Density | @villéger2008 | the functional evenness without weighted abundance [FEve; @villéger2008] |
| Average nearest neighbour distance | $\sqrt{\sum_{i}^{n}{min({q}_{i}-p_{i})^2}})\times \frac{1}{n}$ | Density | @foote1992 | the density of pairs of observations |
| Average displacement | $\frac{\sqrt{\sum_{i}^{n}{({k}_{n})^2}}}{\sqrt{\sum_{i}^{n}{({k}_{n}-Centroid_{k})^2}}}$ | Position | This paper | the ratio between the observations' position from their centroid and the centre of the trait space (coordinates: 0, 0, 0, ...). A value of 1 indicates that the observations' centroid is the centre of the trait space |
Table 3: List of measures with *n* being the number of observations, *d* the total number of dimensions, *k* any specific row in the matrix, *Centroid* being their mean and $\sigma^{2}$ their variance. $\Gamma$ is the Gamma distribution and $\lambda_{i}$ the eigenvalue of each dimension and ${q}_{i}$ and $p_{i}$ are any pairs of coordinates.

We selected these `r metric.names()` space occupancy measures to illustrate how they capture different aspects of space occupancy (not as an expression of our preference).
These measures are specific to Euclidean and isotropic trait spaces (which is not necessary for all measures).
The supplementary material 4 contains the same analysis as described below, performed on 17 measures.
Furthermore, [`moms`](https://tguillerme.shinyapps.io/moms/) allows exploration into the effect of many more measures as well as the customisation of measures by combining them or using user-designed functions.

<!--  -->
<!--  -->
<!--  -->
<!-- COMPARING METRICS -->
<!--  -->
<!--  -->
<!--  -->


## Measure comparisons

We compared the space occupancy measures correlations across all simulations between each pair of measures to assess their captured signal [@villéger2008; @laliberté2010].
We used the measures on the full 13 trait spaces described above.
We then scaled the results and measured the pairwise Pearson correlation to test whether measures were capturing a similar signals or not using the `psych` package [@psych].

<!--  -->
<!--  -->
<!--  -->
<!-- SHIFTING SPACES -->
<!--  -->
<!--  -->
<!--  -->

## Changing space {#changing-spaces}

To assess how the measures responded to changes within trait spaces, we removed 50% of observations each time using the following algorithms:

* **Randomly:** by randomly removing 50% of observations (Fig. 2-A).
This reflects a "null" biological model of changes in trait space: the case when observations are removed regardless of their intrinsic characteristics.
For example, if diversity is reduced by 50% but the space size remains the same, there is a decoupling between diversity and space occupancy [@ruta2013].
Our selected measures are expected to not be affected by this change.

* **<font style="color:black">Size</font>:** by removing observations within a distance from the centre of the trait space lower or greater than a radius $\rho$ (where $\rho$ is chosen such that 50% observations are selected) generating two limit removals: *maximum* and *minimum* (respectively in orange and blue; Fig. 2-B).
This can reflect a strict selection model where observations with trait values below or above a threshold are removed leading to an expansion or a contraction of the trait space.
<font style="color:black">This type of change could be due to habitat destruction [e.g. @mammola2019b] or to mass extinctions [e.g. @wright2017].</font>
Size measures are expected to be most affected by this change.

<!-- [[@@@NOTE - I think the explanation for density is slightly unclear. I think it may also be slightly confusing to have density as the name of a measure and a simulation procedure (same for 'position').]] -->

* **Density:** by removing any pairs of point with a distance $D$ from each other where (where $D$ is chosen such that 50% observations are selected) generating two density removals: *high* and *low* (respectively in orange and blue; Fig. 2-C).
This can reflect changes within groups in the trait space due to ecological factors [e.g. niche repulsion resulting in lower density; @grant2006].
<font style="color:black">This type of change could be due to accelerated rates of evolution [@close2015] or to differences in modes of life in macroevolution [e.g. @healy2019].</font>
Density measures are expected to be most affected by this change.

* **Position:** by removing points similarly as for **<font style="color:black">Size</font>** but using the distance from the furthest point from the centre generating two position removals: *positive* and *negative* (respectively in orange and blue; Fig. 2-D).
This can reflect global changes in trait space (e.g. if an entire group remaining diverse but occupying a different niche).
<font style="color:black">This type of change could be due changes in evolutionary trajectories [@endler2005] or to differences in ecosystem compositions [e.g. @jones2015].</font>
Position measures are expected to be most affected by this change.

The algorithm to select $\rho$ or $D$ is described in the Supplementary material 1.


```{r fig_reduce_space, echo = FALSE, fig.height = 12, fig.width = 8, results = 'hide', fig.cap = paste("different type of space reduction. Each panel displays two groups of 50% of the data points each. Each group (orange and blue) are generated using the following algorithm: A - randomly (the removed elements are displayed in black and the analysed ones in grey); B - by size (maximum and minimum limit); C - by density (high and low); and D - by position (positive and negative). Panel E et F represents two typical display of the reduction results displayed in Table 5: the dots represent the median space occupancy values across all simulations for each scenario of trait space change (Table 2), the solid and dashed line respectively the 50% and 95% confidence intervals. Results in grey are the random 50% reduction (panel A). Results in blue and orange represent the opposite scenarios from panels B, C, and D. The displayed value is the amount of overlap (Bhattacharrya Coefficient) between the blue or orange distributions and the grey one. Panel E and F shows respectively the \"ideal\" and \"worst\" results for any type of measures, where the space occupancy measurement respectively manages or fails to captures a specific type of reduction (i.e. size, position or density; Table 5).")}

set.seed(42)

## Change colors
defaults$col1 <- "blue"
defaults$col2 <- "orange"

op <- par(mfrow = c(3,2), bty = "n")
## The elements
elements <- 300
## The amount to remove
remove <- 0.5

## Trait space
trait_space <- space.maker(elements, 2, distribution = rnorm)

## Random removal colours
defaults$col1 <- "black"
defaults$col2 <- "grey"

## The randomly removed space
random_rm <- reduce.space(trait_space, type = "random", remove = 0.5)
## Plotting the reduction
plot.space(trait_space, random_rm, main = "A - Random removal", defaults)

## Back to normal colours
defaults$col1 <- "blue"
defaults$col2 <- "orange"

## The limit removal
limit_rm <- reduce.space(trait_space, type = "limit", remove = 0.5)
## Plotting the reduction
plot.space(trait_space, limit_rm, main = "B - Limit", defaults)

## The density removal
density_rm <- reduce.space(trait_space, type = "density", remove = 0.5)
## Plotting the reduction
plot.space(trait_space, density_rm, main = "C - Density", defaults)

## The displacement removal
displace_rm <- reduce.space(trait_space, type = "displacement", remove = 0.5)
## Moving the space to the upper right corner (positive)
plot.space(trait_space, displace_rm, main = "D - Position", defaults)


## Test example (positive)
set.seed(42)
colours <- c("grey", "orange", "blue")
par(bty = "n")
## Plot size
plot(NULL, ylim = c(0.8, 3.2), pch = 19, xlim = c(-1,1), xlab = "Centred and scaled space occupancy", ylab = "", xaxt = "n", yaxt = "n", main = "E - Space reduction \"ideal\" results example (Table 5)")
## Adding lines
abline(v = 0, lty = 2, col = colours, lwd = 2)
## Adding the x axis
axis(1, at = c(-1, -0.5, 0, 0.5, 1), labels = TRUE, tick = TRUE, col.ticks = "black", col = "black", lwd = 2)
## Adding the values
values <- cbind(rnorm(30, mean = 0, sd = 0.15),
                rnorm(30, mean = -0.75, sd = 0.15),
                rnorm(30, mean = 0.75, sd = 0.2))
quantile_vals <- apply(values, 2, quantile, probs = c(0.025, 0.250, 0.750, 0.975))
centtend_vals <- apply(values, 2, median)

## Loop through the lines
for(column in 1:3) {
    ## Get the x y values
    line_x_vals <- quantile_vals[, column]
    line_y_vals <- rep(column, 2)

    ## Add the lines
    n_cis <- 4
    for(ci in 1:(n_cis/2)) {
        lines(x = line_x_vals[c(ci, n_cis-(ci-1))], y = line_y_vals, lty = (n_cis/2 - ci + 1), lwd = ci * 1.5 * 2, col = colours[column])
    }
}
## Add the points
points(x = centtend_vals, y = 1:ncol(values), pch = 19, col = colours, cex = 1.5 + 2)

bc_1 <- dispRity::bhatt.coeff(values[,2], values[,1])
bc_2 <- dispRity::bhatt.coeff(values[,3], values[,1])

text(x = 0.2, y = 3,  labels = round(bc_2, 3), cex = 1.4)
text(x = 0.2, y = 2.9,  labels = paste0("\n(Probability of overlap between\nthe blue and the grey distributions)"), cex = 0.8)
text(x = -0.2, y = 2,  labels = round(bc_1, 3), cex = 1.4)
text(x = - 0.2, y = 1.9,  labels = paste0("\n(Probability of overlap between\nthe orange and the grey distributions)"), cex = 0.8)

## Test example (negative)
set.seed(42)
colours <- c("grey", "orange", "blue")
par(bty = "n")
## Plot size
plot(NULL, ylim = c(0.8, 3.2), pch = 19, xlim = c(-1,1), xlab = "Centred and scaled space occupancy", ylab = "", xaxt = "n", yaxt = "n", main = "F -Space reduction \"worst\" results example (Table 5)")
## Adding lines
abline(v = 0, lty = 2, col = colours, lwd = 2)
## Adding the x axis
axis(1, at = c(-1, -0.5, 0, 0.5, 1), labels = TRUE, tick = TRUE, col.ticks = "black", col = "black", lwd = 2)
## Adding the values
values <- cbind(rnorm(30, mean = 0, sd = 0.15),
                rnorm(30, mean = 0, sd = 0.17),
                rnorm(30, mean = 0, sd = 0.2))
quantile_vals <- apply(values, 2, quantile, probs = c(0.025, 0.250, 0.750, 0.975))
centtend_vals <- apply(values, 2, median)

## Loop through the lines
for(column in 1:3) {
    ## Get the x y values
    line_x_vals <- quantile_vals[, column]
    line_y_vals <- rep(column, 2)

    ## Add the lines
    n_cis <- 4
    for(ci in 1:(n_cis/2)) {
        lines(x = line_x_vals[c(ci, n_cis-(ci-1))], y = line_y_vals, lty = (n_cis/2 - ci + 1), lwd = ci * 1.5 * 2, col = colours[column])
    }
}
## Add the points
points(x = centtend_vals, y = 1:ncol(values), pch = 19, col = colours, cex = 1.5 + 2)

bc_1 <- dispRity::bhatt.coeff(values[,2], values[,1])
bc_2 <- dispRity::bhatt.coeff(values[,3], values[,1])

text(x = 0.5, y = 3,  labels = round(bc_2, 3), cex = 1.4)
text(x = -0.6, y = 2,  labels = round(bc_1, 3), cex = 1.4)

```

Because occupancy measures are dependent on the space, we scaled and centred them between -1 and 1 to make them comparable (by subtracting the observed occupancy without reduction to all the measures of the reduced spaces and then divided it by the maximum observed occupancy).
A value of 0 indicates no effect of the space reduction and $>0$ and $<0$ respectively indicates an increase or decrease in the measure value.
We then measured the amount of overlap between the non-random removals (<font style="color:black">size</font>, density and position) and the random removals using the Bhattacharrya Coefficient [@bhattacharyya1943].

### Measuring the effect of space and dimensionality

Distribution differences and the number of dimensions can have an effect on the measure results.
For example, in a normally distributed space, an increase in density can often lead to a decrease in size (though this is not necessarily true if the space is log-normal or uniform).
High dimensional spaces (>10) are subject to the "curse of multidimensionality" [@cursedimensionality]: data becomes sparser with increasing number of dimensions.
This can have two main consequences: 1) the probability of overlap between two groups decreases as a product of the number of dimensions; and 2) the amount of samples needed to "fill" the spaces increases exponentially [see this interactive illustration by Toph Tucker](https://observablehq.com/@tophtucker/theres-plenty-of-room-in-the-corners).
The "curse" can make the interpretation of high dimensional data counter-intuitive.
For example if a group expands in multiple dimensions (i.e. increase in size), the actual hypervolume ($\prod_{i}^{d} range_{Di}$) can decrease (Fig. 3 and Tables 6, 7).

We measured the effect of space distribution and dimensionality using an ANOVA ($occupancy \sim distribution$ and $occupancy \sim dimensions$) by using all spaces with 50 dimensions and the uniform and normal spaces with equal variance and no correlation with 3, 15, 50, 100 and 150 dimensions (Table 2) for testing respectively the effect of distribution and dimensions.
The results of the ANOVAs (F and *p*-values) are reported in Table 5 (full results in supplementary material 3).


<!--  -->
<!--  -->
<!--  -->
<!-- EMPIRICAL EXAMPLES -->
<!--  -->
<!--  -->
<!--  -->


## Empirical examples

We analysed the effect of the different space occupancy measures on six different empirical studies covering a range of fields that employ trait space analyses.
For each of these studies we generated trait spaces from the data published with the papers.
We divided each trait spaces into two biologically-relevant groups and tested whether the measures differentiated the groups in different ways.
Both the grouping and the questions were based on a simplified version of the topics of these papers (with no intention to re-analyse the data and questions).
The procedures to generate the data and the groups varies between studies and is detailed in the supplementary materials 2.


study | field | taxonomic group | traits | trait space | size | groups | question |
----------|-----------|-----------|--------------|--------------|----------|--------------|----------------
@beck2014 | Palaeontology | Mammalia | discrete morphological phylogenetic data | Ordination of a distance matrix (PCO) | 106*105 | 52 crown vs. 54 stem | Are crown mammals more disparate than stem mammals?|
@wright2017 | Palaeontology | Crinoidea | discrete morphological phylogenetic data | Ordination of a distance matrix (PCO) | 42*41 | 16 before vs. 23 after | Is there a difference in disparity before and after the Ordovician mass extinction?|
@marcy2016 | Evolution | Rodentia | skull 2D landmark coordinates  | Ordination of a Procrustes Superimposition (PCA) | 454*134 | 225 *Megascapheus* vs. 229 *Thomomys* | Are two genera of gopher morphologically distinct? |
@hopkins2016 | Evolution | Trilobita | 3D landmark coordinates | Ordination of a Procrustes Superimposition (PCA) | 46*46 | 36 adults vs. 10 juveniles | Are juvenile trilobites a subset of adult ones <font style="color:black">in trait space</font>? |
@jones2015 | Ecology | Plantae | Communities species compositions | Ordination of a Jaccard distance matrix (PCO) | 48*47 | 24 aspens vs. 24 grasslands | Is there a difference in species composition between aspens and grasslands? |
@healy2019 | Ecology | Animalia | Life history traits | Ordination of continuous traits (PCA) | 285*6 | 83 ecthotherms vs. 202 endotherms | Do endotherms have more diversified life history strategies than ectotherms? |
Table 4: details of the six empirical trait spaces.

For each empirical trait space we bootstrapped each group 500 times [@disprity] and applied the `r metric.names()` space occupancy measure to each pairs of groups.
We then compared the means of each groups using the Bhattacharrya Coefficient [@bhattacharyya1943].

<!--  -->
<!--  -->
<!--  -->
<!-- RESULTS SECTION -->
<!--  -->
<!--  -->
<!--  -->








# Results

```{r loading_results, echo = FALSE}
## Loading the results
remove_05 <- load.results("remove_05")
```

```{r running_tests, echo = FALSE}
## Anova function
anova.fun <- function(data) {return(aov(glm(disparity ~ factor, data = data)))}

## Running the tests
all_test <- test.simulation(remove_05, test = anova.fun, scale = TRUE)
all_dim_test <- test.simulation(remove_05, test = anova.fun, scale = TRUE,
                      factors = c("uniform3", "uniform15", "uniform50", "uniform100", "uniform150",
                                  "normal3", "normal15", "normal50", "normal100", "normal150"))
space_test <- test.simulation(remove_05, test = anova.fun, scale = TRUE,
                              factors = c("uniform50", "uniform50c", "normal50", "normal50c",
                                          "random50", "pca_like", "pco_like"))
```

## Measure comparisons 


```{r fig_measure_correlation, echo = FALSE, fig.height = 8, fig.width = 8, results = 'hide', fig.cap = paste("pairwise correlation between the scaled measures. Numbers on the upper right corner are the Pearson correlations. The red line are linear regressions (with the confidence intervals in grey). Av.: average; dist.: distance; min.: minimum; span.: spanning.")}

## Adding proper names
results_pairwise <- remove_05
names(results_pairwise[[1]]) <- metric_names
## Shortening the names for plotting
names(results_pairwise[[1]])[1] <- "Av. dist.\nfrom centroid"
names(results_pairwise[[1]])[5] <- "Min. span. tree\nav. dist."
names(results_pairwise[[1]])[6] <- "Min. span. tree\ndist. evenness"
names(results_pairwise[[1]])[7] <- "Av. nearest\nneighbour dist."
names(results_pairwise[[1]])[8] <- "Av. displacements"

## Plotting the pairwise results
pairwise.plot(results_pairwise, scale = TRUE, type = "base", plot = "cor")
```

<font style="color:black">Most measures of space were positively correlated (Pearson correlation of 0.99 for the average Euclidean distance from centroid and sum of variance or 0.97 for the average nearest neighbour distance and minimum spanning tree average length; Fig. 3).
The remaining measures were either somewhat correlated or had a negative pairwise distribution </font> (ranging from 0.66 for the sum of variances and the ellipsoid volume to -0.09 between the average displacement and the average Euclidean distance from centroid; Fig. 3).
All measures but the ellipsoid volume were normally (or nearly normally) distributed (Fig. 3).

## Space shifting

```{r fable_results, fig.show='hide', echo=FALSE, fig.height=3, fig.width=3}
## The measures names (shortened vector)
name <- metric_names

## Making a list of parameters for each mini plot
plot.param <- list(scaler = 3,
                   bg.col = "black",
                   col = c("grey", "orange", "blue"),
                   quantiles = c(95, 50),
                   cent.tend = median,
                   pch = 19,
                   metric.max = length(metrics_list),
                   cex = 2)

## Looping through each mini plot for every measure
for(metric in 1:length(metrics_list)) {
    generate.fable.plot(data = remove_05, metric = metric, what = "limits", plot.param = plot.param, overlap = TRUE)
    generate.fable.plot(data = remove_05, metric = metric, what = "densit", plot.param = plot.param, overlap = TRUE)
    generate.fable.plot(data = remove_05, metric = metric, what = "displa", plot.param = plot.param, overlap = TRUE)
}

```

```{r, echo = FALSE, eval = FALSE}
## Function for printing the table in the R console
print.fable <- function(n_measures, byrow, ncol, test) {
  ## Make the plot.id table 
  ids <- matrix(1:(n_measures*ncol), ncol = ncol, byrow = byrow)

  for(one_measure in 1:n_measures) {
    ## Print all the rows one by one
    text <- c(paste0("`r name[", one_measure, "]`"), paste0("`r plot.id(", ids[one_measure, ], ")`"))
    ## Add some tests?
    if(test) {
      text <- c(text, paste0("`r s.test(", one_measure, ", \"s\")`"), paste0("`r s.test(", one_measure, ", \"r\")`"))
    }
    ## Print the line
    cat(paste(text, collapse = " | "))
    cat("|\n")
  }
}
## Getting the fable to copy paste under the table header below.
print.fable(length(metrics_list), byrow = TRUE, ncol = 3, test = TRUE)
```

Measure      | Size change  | Density change | Position change | Distribution effect | Dimensions effect |
:-----------|----------------|-----------------|----------------|---------------------|-------------------|
`r name[1]` | `r plot.id(1)` | `r plot.id(2)` | `r plot.id(3)` | `r s.test(1, "s")` | `r s.test(1, "r")`|
`r name[2]` | `r plot.id(4)` | `r plot.id(5)` | `r plot.id(6)` | `r s.test(2, "s")` | `r s.test(2, "r")`|
`r name[3]` | `r plot.id(7)` | `r plot.id(8)` | `r plot.id(9)` | `r s.test(3, "s")` | `r s.test(3, "r")`|
`r name[4]` | `r plot.id(10)` | `r plot.id(11)` | `r plot.id(12)` | `r s.test(4, "s")` | `r s.test(4, "r")`|
`r name[5]` | `r plot.id(13)` | `r plot.id(14)` | `r plot.id(15)` | `r s.test(5, "s")` | `r s.test(5, "r")`|
`r name[6]` | `r plot.id(16)` | `r plot.id(17)` | `r plot.id(18)` | `r s.test(6, "s")` | `r s.test(6, "r")`|
`r name[7]` | `r plot.id(19)` | `r plot.id(20)` | `r plot.id(21)` | `r s.test(7, "s")` | `r s.test(7, "r")`|
`r name[8]` | `r plot.id(22)` | `r plot.id(23)` | `r plot.id(24)` | `r s.test(8, "s")` | `r s.test(8, "r")`|
Table 5: Results of the effect of space reduction, space dimension distributions and dimensions number of the different space occupancy measures. The dots represent the median space occupancy values across all simulations for each scenario of trait space change (Table 2), the solid and dashed line respectively the 50% and 95% confidence intervals. See Fig. 2 for details on the interpretation of the figures distributions and values. F-values for distribution effect and dimensions effect represents respectively the effect of the ANOVAs space occupancy ~ distributions and space occupancy ~ dimension represent the ratio of sum squared difference within and between groups (the higher, the more the factor has an effect on the measure) and associated _p_-values (0 '\*\*\*' 0.001 '\*\*' 0.01 '\*' 0.05 '.' 0.1 '' 1). This figure illustrates how different measures can be influenced by different aspects of changes in the trait space. E.g. the Average Euclidean distance from centroid (row 1) captures mainly changes in size (column 1), but also captures changes in density (column 2) but does not capture changes in position (column 3).

As expected, some different measures capture different aspects of space occupancy.
However, it can be hard to predict the behaviour of each measure when 50% of the observations are removed.
We observe a clear decrease in <font style="color:black">the median measure value</font> in less than a third of the space reductions (10/36).

In terms of change in size, only the average Euclidean distance from centroid and the sum of variances seem to capture a clear change in both directions.
<!-- However, the increase in size does not correspond to an *actual* increase in volume in the trait space (i.e. the size from the blue observations in Fig. 2-B is equivalent to the one in Fig. 2-A). -->
In terms of change in density, only the minimum spanning tree average distance and the average nearest neighbour distance seem to capture a clear change in both directions.
And in terms of change in position, only the average displacement metric seems to capture a clear change in direction (albeit not in both directions).
This is not surprising, since the notion of positions becomes more and more complex to appreciate as dimensionality increases (i.e. beyond left/right, up/down and front/back).

## Empirical example

```{r test_empirical, echo = FALSE}
## Loading the results
empirical_results <- load.results("empirical_results")

## Testing the differences for each distributions
bhatt.coeff.safe <- function(x, y, tol = 1e-16, ...) {
  if(all(x < tol) | all(y < tol)) {
    bhatt.coeff(1, 0)
  } else {
    bhatt.coeff(x, y, ...)
  }
} 
disparity_test <- lapply(empirical_results, lapply, test.dispRity, test = bhatt.coeff.safe)
```

```{r fable_results_empirical, fig.show='hide', echo=FALSE, fig.height=3, fig.width=3}
## Plotting parameters
plot.param <- list(cex = 2,
                   col = c("#F7B27E", "#BFE4E3"),
                   border = c("#F65205", "#3E9CBA"),
                   na.cex = 3,
                   scaler = 3)

## Looping through each mini plot for every measure
data_names <- c("Beck and Lee 2014", "Wright 2017", "Marcy et al. 2016",
                "Hopkins and Pearson 2016", "Jones et al. 2015", "Healy et al. 2019")

for(dataset in 1:length(data_names)){
    for(measure in 1:length(metrics_list)){
        ## Plotting the results
        generate.fable.empirical(data = empirical_results[[dataset]][[measure]],
                                 test = disparity_test[[dataset]][[measure]],
                                 precision = 1e-5, plot.param, dataset = dataset)
    }
}
```

```{r change_plotid_chain, echo = FALSE, eval = TRUE}
## Changing defaults
body(plot.id)[[3]] <- substitute(chain <- "fable_results_empirical")
```

```{r, echo = FALSE, eval = FALSE}
print.fable(length(metrics_list), byrow = FALSE, ncol = 6, test = FALSE)
```

Measure     | Beck and Lee 2014  | Wright 2017 | Marcy et al. 2016 | Hopkins and Pearson 2016 | Jones et al. 2015 | Healy et al. 2019  |
:-----------|----------------|-----------------|------------------|----------------|------------------|----------------|
Comparisons (orange *vs.* blue) | crown *vs.* stem mammals morphologies | crinoids morphologies before *vs.* after the end-Ordovician extinction | *Megascapheus* *vs.* *Thomomys* skull shapes | adults *vs.* juveniles trilobites cephalon shapes | aspens *vs.* grasslands communities compositions | ecthotherms *vs.* endotherms life history traits |
`r name[1]` | `r plot.id(1)` | `r plot.id(9)` | `r plot.id(17)` | `r plot.id(25)` | `r plot.id(33)` | `r plot.id(41)`|
`r name[2]` | `r plot.id(2)` | `r plot.id(10)` | `r plot.id(18)` | `r plot.id(26)` | `r plot.id(34)` | `r plot.id(42)`|
`r name[3]` | `r plot.id(3)` | `r plot.id(11)` | `r plot.id(19)` | `r plot.id(27)` | `r plot.id(35)` | `r plot.id(43)`|
`r name[4]` | `r plot.id(4)` | `r plot.id(12)` | `r plot.id(20)` | `r plot.id(28)` | `r plot.id(36)` | `r plot.id(44)`|
`r name[5]` | `r plot.id(5)` | `r plot.id(13)` | `r plot.id(21)` | `r plot.id(29)` | `r plot.id(37)` | `r plot.id(45)`|
`r name[6]` | `r plot.id(6)` | `r plot.id(14)` | `r plot.id(22)` | `r plot.id(30)` | `r plot.id(38)` | `r plot.id(46)`|
`r name[7]` | `r plot.id(7)` | `r plot.id(15)` | `r plot.id(23)` | `r plot.id(31)` | `r plot.id(39)` | `r plot.id(47)`|
`r name[8]` | `r plot.id(8)` | `r plot.id(16)` | `r plot.id(24)` | `r plot.id(32)` | `r plot.id(40)` | `r plot.id(48)`|
Table 6: Comparisons of pairs of groups in different empirical trait spaces. NAs are used for cases where space occupancy could not be measured due to the curse of multidimensionality. The displayed values are the amount of overlap between both groups (Bhattacharrya Coefficient).

<font style="color:black">As with the as for the simulations, there is no measure that summarises all the aspects of distributions for empirical data.</font>
For all `r metric.names()` measures (except the ellipsoid volume) we see either one group or the other having a bigger mean than the other and no consistent case where a group has a bigger mean than the other for all the measures.
For example, in the @beck2014's dataset, there is a clear difference in size using the average Euclidean distance from centroid or the sum of variances (overlaps of respectively 0.175 and 0.159) but no overlap when measuring the size using the sum of ranges (0.966).
However, for the @hopkins2016's dataset, this pattern is reversed (no clear differences for the average Euclidean distance from centroid or the sum of variances - 0.701 and 0.865 respectively - but a clear difference for the sum of ranges (0).
For each dataset, the absolute differences between each groups is not consistent depending on the measures.
For example, in @hopkins2016's dataset, the orange group's mean is clearly higher than the blue one when measuring the sum of ranges (0) and the inverse is true when measuring the average displacement (0).


<!--  -->
<!--  -->
<!--  -->
<!-- DISCUSSION SECTION -->
<!--  -->
<!--  -->
<!--  -->








# Discussion

Here we tested 25 measures of trait space occupancy on simulated and empirical datasets to assess how each measure captures changes in trait space size, density and position.
Our results show that the correlation between measures can vary both within and between measure categories (Fig. 3), highlighting the importance of understanding the measure classification for the interpretation of results.
Our simulations show that different measures capture different types of trait space change (Table 5), meaning that the use of multiple measures is important for comprehensive interpretation of trait space occupancy.
We also show that the choice of measure impacts the interpretation of group differences in empirical datasets (Table 6).


#### Measures comparisons

Measures within the same category of trait space occupancy (size, density or position) do not have the same level of correlation with each other.
For example, the average Euclidean distance from centroid (size) is highly correlated to the sum of variances (size - correlation of 0.99) and somewhat correlated with the minimum spanning tree average distance (density - correlation of 0.66) but poorly with the ellipsoid volume (size - correlation of 0.17) and the minimum spanning tree distances evenness (density - correlation of -0.05).
Furthermore, the fact that we have such a range of correlations for normal distributions suggests that each measure can capture different summaries of space occupancy ranging from obvious differences (for measures not strongly correlated) to subtle ones (for measures strongly correlated).


#### Space shifting

Most measures capture no changes in space occupancy for the "null" (random) space reduction (in grey in Table 5).
This is a desirable behaviour for space occupancy measures since it will likely avoid false positive errors in studies that estimate biological processes from space occupancy patterns [e.g. convergence @marcy2016, life history traits @healy2019].
However, the average nearest neighbour distance and the sum of ranges have a respectively positive and negative "null" median.
<font style="color:black">In itself this is not necessarily a negative</font> property but it should be kept in mind that even random processes can increase or decrease these measures' values.

For changes in size, the sum of variances and the average Euclidean distance from centroid are good descriptors (Table 5).
However, as illustrated in the 2D examples in Fig. 2-B only the blue change results (Table 5) should not result in a direct change in overall size because the trait space is merely "hollowed" out.
That said, "hollowing" is harder to conceptualise in many dimensions and the measures can still be interpreted for comparing groups (orange has a smaller volume than blue).

The average nearest neigbhour distance and the minimum spanning tree average distance consistently detect changes in density with more precision for low density trait spaces (in blue in Table 5).
However, we can observe some degree of correlation between the changes in density and the changes in size for most measure picking either signal.
This could be due to the use of normally distributed spaces where a change in density often leads to a change in size.
This is not necessarily the case with empirical data.

Regarding the changes in position, only the average displacement measure seems able to distinguish between a random change and a displacement of the trait space (Table 5).
However, the average displacement measure does not distinguish between positive or negative displacement: this might be due to the inherent complexity of *position* in a multidimensional trait space.

#### Empirical examples

Although most differences are fairly consistent within each dataset with one group having a higher space occupancy score than the other for multiple measures, this difference can be more or less pronounced within each dataset (ranging from no to nearly full overlap - BC $\in(0;0.995)$) and sometimes even reversed.
This indicates that opposite conclusions can be drawn from a dataset depending on which space occupancy measure is considered.
For example, in @wright2017, crinoids after the Ordovician mass extinction have a higher median measure value for all measures but for the average displacement.
These differences depending on the measures are also more pronounced in the empirical datasets where the observations per group are unequal [@hopkins2016; @healy2019].

### Caveats

While our simulations are useful to illustrate the behaviour of diverse space occupancy measures, they have several caveats.
First, the simulated observations in the trait spaces are independent.
This is not the case in biology where observations can be spatially [@jones2015] or phylogenetically correlated [e.g. @beck2014].
Second, the algorithm used to reduce the trait <font style="color:black"> spaces might </font> not always accurately reflect changes.
This might favour some specific measures over others, in particular for the changes in density that modify the nearest neighbour density rather than changing the global density.
This algorithmic choice was made in order to not confound changes in density along with changes in size.
However, the results presented here probably capture the general behaviour of each measure since results are consistent between the simulated and empirical analysis.

<font style="color:black"> Furthermore, we did not take into account the effect of sampling on space occupancy measurements (but see additional results with 80% and 20% space reduction in the supplementary materials 4).
In fact, sampling has been previously shown to have an effect on measurements depending on range or volumes (e.g. the sum of ranges or the hypervolume @ciampaglio2001).
This effect is especially expected to be acerbated in macroevolutionary studies when using the fossil record [@brocklehurst2013] but can be tackled using rarefaction and bootstrapping techniques [@disprity]. </font>

### Using `moms` to choose the appropriate measurements

<font style="color:black">Therefore, we propose the [`moms`](https://tguillerme.shinyapps.io/moms/) shiny app to allow workers to help them choose their set of space occupancy measurements (and test the caveats mentioned above).
`moms` is an online graphical user interface to help analyse multidimensional data.
It allows users to upload their dataset of interest (or simulate one with specific parameters) and measure space occupancy using a variety of implemented measures (namely, but not only, the ones used in this study).
Furthermore, the package allows simulation of shifts in trait space occupancy as also presented in this paper to test whether some measures capture specific changes in space.
However, `moms` is not a tool for analysing multidimensional data _per se_ but rather for helping workers to chose the space occupancy measure most appropriated to their data and question.
To run multidimensional analysis, we suggest using dedicated `R` packages (such as - but not limited to: @oksanen2007vegan, @momocs, @bat2015, @disprity). </font>


### Conclusions

<font style="color:black">We insist that although no measure is objectively better than the next one, some can be more problematic than other in specific contexts. 
For example, the results for the Sum of Ranges, Minimum spanning tree average distances, and to a lesser extent average nearest neighbour distances produced results in the reduced space often similar to the randomly reduced spaces (Table 5).
This does not make them "bad" measures but rather heavily context dependent.
Regardless, we believe that workers should identify the most appropriate measures based on their trait space properties as well as their specific biological question.
We believe this could be fostered by following these several suggestions: </font>

First, we suggest using multiple measures to tackle different aspects of the trait space. 
This follows the same logical thinking that the mean might not be sufficient to describe a distribution (e.g. the variance might be a good additional descriptor).
Although using multiple measures is not uncommon in macroevolutionary studies [e.g. @halliday2015] or in ecology [@mammola2019], they often do no cover more than one of the three categories of trait space measures <font style="color:black">(but see the recent work of @carmona2019 and @mammola2020).</font>

Second, we suggest selecting the measures that best address the biological question at hand.
If one studies an adaptive radiation in a group of organisms, it is worth thinking what would be the expected null model: would the group's size increase (radiation in all directions), would it increase in density (niche specialisation) or would it shift in position (radiation into a new set of niches)?

Third, we suggest not naming measures after the biological aspect they describe which can be vague (e.g. "disparity" or "functional dispersion") but rather after what they are measuring and why (e.g. "we used sum of ranges to measure the space size").
We believe this will support both a clearer understanding of what *is* measured as well as better communication between ecology and evolution research where measures can be similar but have different names.

Multidimensional analyses have been acknowledged as essential tools in modern biology but they can often be counter-intuitive [@cursedimensionality].
It is thus crucial to accurately describe patterns in multidimensional trait spaces to be able to link them to biological processes.
When summarising trait spaces, it is important to remember that a pattern captured by a specific space occupancy measure is often dependent on the properties of the space and of the particular biological question of interest.
We believe that having a clearer understanding of both the properties of the trait space and the associated space occupancy measures (e.g. using [`moms`](https://tguillerme.shinyapps.io/moms/)) as well as using novel space occupancy measures to answer specific questions will be of great use to study biological processes in a multidimensional world.


# Acknowledgements

We thank Natalie Jones and Kevin Healy for helping with the empirical datasets and Stefano Mammola and Neil Brocklehurst for their positive and encouraging reviews.
We acknowledge funding from the Australian Research Council DP170103227 and FT180100634 awarded to VW.

# Authors contributions

TG, MNP, AEM and VW designed the project.
TG and AEM collected the empirical dataset.
TG ran the analyses and designed the software.
TG, MNP, AEM and VW wrote the manuscript.

# Data Availability, repeatability and reproducibility

The raw empirical data is available from the original papers [@beck2014; @jones2015, @marcy2016; @hopkins2016; @wright2017; @healy2019].
The subsets of the empirical data used in this analysis are available on figshare 
[DOI: 10.6084/m9.figshare.9943181.v1](https://doi.org/10.6084/m9.figshare.9943181.v1).
The modified empirical data are available in the package accompanying this manuscript (`data(moms::demo_data)`).
This manuscript (including the figures, tables and supplementary material) is repeatable and reproducible by compiling the vignette of the [GitHub `moms R` package](https://github/TGuillerme/moms).

# Conflict of interest
None declared.

# References