funspace: An R package to build, analyse and plot functional trait spaces

Functional trait space analyses are pivotal to describe and compare organisms' functional diversity across the tree of life. Yet, there is no single application that streamlines the many sometimes‐troublesome steps needed to build and analyse functional trait spaces.

Subsequent research has attempted to characterize such axes using functional traits for fish (Winemiller & Rose, 1992), corals (Darling et al., 2012), and vascular plants (e.g.Diaz et al., 2004), among other groups of organisms.The implication of having a handful of axes describing species' strategies implies that the diversity of organisms' form and function must be constrained by trade-offs among traits (Pianka, 1970 for animals; Greenslade, 1983 for insects; Westoby et al., 2002 for vascular plants).Nowadays, we know that pervasive trade-offs between traits constrain organisms' form and function along (relatively) few trait dimensions that are often independent, thus limiting possible trait combinations within two -or highly dimensional spaces, also called functional trait spaces, across the tree of life (Bueno et al., 2023;Carmona, Bueno, et al., 2021;Carmona, Tamme, et al., 2021;Chelli et al., 2024;Díaz et al., 2016;Junker et al., 2023;Mouillot et al., 2021;Westoby et al., 2021;Winemiller et al., 2015).Expressing the diversity of organisms' form and function within functional trait spaces of a reduced number of dimensions allows for robust generalization and quantification of species' strategic axes of trait variation (see Mouillot et al., 2021).Thus, functional trait spaces represent indispensable tools to describe organisms' functional diversity across the tree of life.
Different approaches allow building multivariate functional trait spaces.The most common quantitative tool is Principal Component Analysis (PCA), as it allows reducing a trait dataset to few independent trait dimensions defined by the inherent relationships between traits.However, other approaches for dimensionality reduction that allow considering both quantitative and categorical trait data, such as Principal Coordinates Analysis (PCoA) or nonmetric multidimensional scaling (NMDS) are also widely used.Recent methodological advances have permitted to explore additional features of functional trait spaces besides the trait dimensions that define them.
In particular, kernel density analyses (Carmona et al., 2016(Carmona et al., , 2019;;Duong, 2007) have revealed that species occupy functional trait spaces differentially, resulting in species clumping around some trait combinations that are much more frequent than others (Carmona, Bueno, et al., 2021;Díaz et al., 2016;Puglielli et al., 2021 for vascular plants; Carmona, Tamme, et al., 2021;Cox et al., 2021;Toussaint et al., 2021 for different animal groups).As a result, kernel density estimates are nowadays widely employed to estimate different properties about the way in which different species and groups of species occupy a functional trait space, including aspects such as the amount of space occupied (functional richness) or the degree to which species in a group have different trait combinations (functional divergence), among other indices of functional diversity (see Mammola et al., 2021 for a review).
A first pitfall arising when using PCA (but also distance-based ordinations approaches depending on the dissimilarity metric considered) is that missing values are not allowed in the dataset.This is especially the case of large-scale analyses, where the unprecedented trait availability provided by global databases contrasts with the use of disparate trait frameworks across studies and disciplines, resulting in large databases including many traits but with missing information for many of the species (e.g. the TRY plant trait database, Kattge et al., 2020).This has pushed research towards developing imputation methods (e.g.Penone et al., 2014;Schrodt et al., 2015;Stekhoven & Bühlmann, 2012) that mostly use machine-learning techniques to impute missing trait values based on trait-trait correlations sometimes improved by accounting for species' phylogenetic information (Penone et al., 2014).A key point to consider is that missing trait information is not distributed randomly within datasets.
For example, across mammals, big species with large range areas are generally better informed than small species or species with small ranges (González-Suárez et al., 2012); similar biases have been described for plants (Carmona, Bueno, et al., 2021;Sandel et al., 2015).
Recent research is increasingly suggesting that imputation can largely correct these biases, so that functional diversity patterns inferred from imputed datasets are much closer to the real ones than those that would be estimated using only species with complete trait information (Penone et al., 2014;Stewart et al., 2023).
Even when there are no missing data, another pitfall is deciding how many dimensions to retain to maximize the information contained in the dataset while minimizing information redundancy.
In the case of PCA, this is usually done by retaining the first two Principal Components (PCs) since they capture most of the variance in a dataset.However, if the intrinsic dimensionality of the considered trait dataset is higher, this approach might lead to the loss of biologically relevant information (Laughlin, 2014).Several methods have been proposed to determine the number dimensions (Mouillot et al., 2021;Peres-Neto et al., 2005).In the context of PCA, the parallel analysis method proposed by Horn (1965) is one of the most precise methods in identifying the correct number of PCs to retain (Peres-Neto et al., 2005), and one of the most widely used.This method contrasts eigenvalues produced through a PCA on n random data sets of uncorrelated variables with the same dimension of the original dataset to produce eigenvalues for components that are adjusted for the sample error-induced inflation.Apart from some large cross-species studies (e.g.Carmona, Bueno, et al., 2021;Carmona, Tamme, et al., 2021;Díaz et al., 2016;Guillemot et al., 2022;Toussaint et al., 2021), the two described pitfalls of using PCA are hardly considered when building functional trait spaces.We argue that this might mostly depend on lack of knowledge of the existence of specific procedures, and likely, on the potential difficulties when implementing them.However, the most common practices when building functional trait spaces remain choosing the first two PCs and/or dropping observations with missing values, with the inevitable loss of biological information.
Another aspect that deserves attention is how to use functional trait spaces as a ground to test the relationship between species' trait strategies and additional variables.This includes analyses such as exploring the relationship between species' ecological strategies (i.e.their position in the functional trait space) and climatic features at a species' habitat, or to explore how extinction risk (Carmona, Tamme, et al., 2021), as well as any other adaptive syndrome (Chelli et al., 2024;Pavanetto et al., 2023;Puglielli et al., 2022), is related to species traits.However, analysing these multidimensional relationships is inherently complex (Villéger et al., 2011).One layer of complexity is provided by the need to use the multiple dimensions defining the functional trait space as predictors in models attempting to link multidimensional species strategies to other ecological dimensions (e.g.climate).Another problem is that more than often we do not have any previous knowledge on the functional form (e.g.linearity) of the relationship linking the target dimension to species' strategies.In addition, visualizing such relationships is not straightforward, and we often rely on separately analysing relationships between response variables and single trait dimensions (e.g. a single PC).These difficulties can be overcome by considering functional trait spaces as a set of coordinates that we can use to build maps of any response variable of interest.General additive models (GAMs), used to estimate smooth functional relationships between predictor variables and a response (Pedersen et al., 2019), provide a solution to do that.Accordingly, they have been consistently used for fitting and mapping, for instance, species distribution models' predictions in spaces defined by geographical coordinates (e.g.Naimi & Araújo, 2016).GAMs allow expressing the models' predictor as a multidimensional smoother, allowing to use the functional trait space dimensions as a single predictor (see Carmona, Tamme, et al., 2021;Chelli et al., 2024;Pavanetto et al., 2023;Puglielli et al., 2022 for examples).In the case of functional trait spaces, GAMs provide a sound and flexible modelling solution because they operate using piecewise functions adapting to the local conditions within the functional trait space, so that GAMs behaviour in a particular portion of the data point cloud does not overall alter the global behaviour of the model (Pedersen et al., 2019;Venables & Dichmont, 2004).
While this comes at the cost of the interpretability of the coefficients, spline regressions are more easily interpreted from a graphical point of view rather than through the values of their coefficients (Venables & Dichmont, 2004).Thus, GAMs become particularly useful to visualize and test patterns of how target variables vary within functional trait spaces.
Despite functional trait spaces are widely used in disparate research fields spanning ecology, plant science, animal ecology, evolutionary ecology -for example, searching the term 'trait space' in Web of Science returned 11,440 documents -there is no single R package for building, exploring, mapping and plotting functional trait spaces.To fill this gap, we propose funspace, an R package to handle functional trait space analyses using bivariate or multivariate trait data.Due to the very little difference between using funspace with pairs of traits or considering a multivariate ordination as input (Figure 1), here we present the package functionalities in the context of PCA-based analyses, and we refer to funspace documentation and examples for other cases, including raw trait data or other ordination methods (e.g.NMDS).The package consists of three main interconnected modules (Figure 1; Table 1): 1. Building and exploring module: it consists of multiple functions to build the functional trait space from bivariate or multivariate input data using multivariate kernel density methods, and to analyse the main features of a functional trait space (e.g.functional diversity indexes, testing against null models).
If specified, all these analyses can be iterated across levels of a grouping variable (e.g.within groups such as families or populations, etc.).
2. Mapping module: it consists of one function that uses GAMs to statistically test for the link between a target variable at the species level (such as extinction risk) and the position of species within the functional trait space.The function handles multiple groups as well.
3. Plotting module: it consists of a generic plot function that can receive as input the objects built in either of the previous modules.In case the input is the object built in module 1 (building and exploring), the function prints a bivariate (or pairs of dimensions in case of multivariate) functional trait space displaying bivariate kernel density estimates (for single or multiple groups).If the input is an object built in module 2 (mapping module), then the function prints a heatmap depicting how a target variable is distributed within the functional trait space (for single or multiple groups).All plots can be customized and new features can readily be added after plotting.

| MODULE 1: BU ILD ING AND E XPLORING A FUN C TI ONAL TR AIT S PACE
funspace includes three functions to build and explore a functional trait space: funspaceDim(), impute(), and funspace() (Table 1; Figure 1).Here we only illustrate the main features of the functions when the input is a PCA object and refer to the funspace documentation for full details on the functions' arguments and input data.
However, note that the same exact procedures can be applied to any bivariate trait data (i.e.pairs of raw traits, or scores of species on a different type of ordination).An example of how to use the functions constituting module 1 is reported in Box 1. funspaceDim() allows the user to identify the number of dimensions that are needed to build a functional trait space.The identified dimensions are those that minimize redundancy while maximizing the information contained in the trait data based on the method proposed by Horn (1965) (see Introduction).Briefly, this method contrasts the eigenvalues produced through PCAs run on (30 * (number of variables)) random datasets with the same number of variables and observations as the input dataset.Eigenvalues >1 are retained in the adjustment.funspaceDim() implements the Horn test via the paran() function of the R package paran (Dinno, 2018).A trait matrix is the only input needed to run the function.funspaceDim() returns the number of dimensions to be retained from the subsequent PCA (Figure 1) and prints this information in the R console.Note that funspaceDim() can be used only when the user wants to pass a PCA object as argument of funspace().In the case of PCoA as the input object, the dimensions are iteratively identified by the analysis and depend on the choice of the distance measure used (e.g.Euclidean distance) together with the number of input variables.For NMDS, the user can specify the number of dimensions (k argument of vegan::monoMDS(), see below) before running the analysis, and the best value of k needs to be determined using the model stress, a measure of the model's goodness of fit obtained by comparing distances in the original dissimilarity matrix to the fitted one in the ordination space.
If the input trait matrix (either bivariate or multivariate) contains missing data, the user can fill the gaps using the impute() function.
By default, impute() fills the gaps in the dataset using information on trait-trait relationships using random forest models as implemented in the missForest R package (Stekhoven & Bühlmann, 2012).
Optionally, a phylogenetic tree (an object of class phylo) can be included to improve the imputation procedure by accounting for species phylogenetic relatedness together with the information on trait-trait relationships (Penone et al., 2014).Briefly, a given number of phylogenetic eigenvectors derived from the tree (specified by the argument nEigen, default is 10) are added to the original trait matrix, and the resulting dataset is then used in the above-mentioned random forest-based procedure.The addingSpecies argument allows the user to add to the phylogeny those species that are only present in the trait matrix.If TRUE (default is FALSE), the phytools::add.species.to.genus() function (Revell, 2012) is used to add species to the root of the genus if there are congeneric species in the tree.
Those species that cannot be added to the phylogeny are imputed without taking phylogenetic information into account.Users who want to make use of the multiple options available in add.species.
to.genus() should modify their phylogenetic tree beforehand.The impute() output is a list containing the imputed version of the trait matrix and its non-imputed counterpart.
After having defined the number of dimensions that define the functional trait space, and having ran a PCA, PCoA or NMDS, the user can build the functional trait space using funspace(), the F I G U R E 1 funspace workflow by package module.Functions classification according to package modules is highlighted by colour coding.Input data are shown in grey.In case the input trait matrix has missing trait data, the user might want to gap-fill the bivariate or multivariate trait matrix using the impute() function, with the option to account for phylogenetic information in the imputation process.The second step is running an ordination (PCA, PCoA, NMDS) before using the other funspace functions.In the case of a PCA object, the user can select the number of principal components to be retained in the analysis, by using a the multivariate trait matrix to funspaceDim().The doubleheaded arrows indicate that the dimensions retained using fuspaceDim() are those that should be extracted from the subsequent PCA as well as those used to build the trait space using funspace().In cases where the user decides to identify the trait space dimensions using funspaceDim(), and to pass the target PCs to funspace() instead of the PCA object, at this point the workflow follows the same path as that of any bivariate trait matrix.That is, the user can directly use a bivariate trait matrix (including coordinates obtained from any ordination method) to run funspace().Because funspace is strongly focused in providing graphical representations of functional spaces, in cases when the functional space has higher dimensionalities, two dimensions at a time have to be selected within the funspace() function.Finally, an object TPD::TPDs() or TPD::TPDMeans() can be passed as argument of the funspace() function, and we refer to Carmona et al. (2016Carmona et al. ( , 2019) ) for building TPDs objects.Once the funspace() object has been defined, the user can plot the output (i.e.bivariate functional trait space with kernel density estimates) using the plot() function.In case the user wants to test the observed functional trait space against a multivariate normal or a uniform distribution, the funspace() object can be passed to funspaceNull().funspaceGAM() allows assessing how an additional target variable (e.g.temperature, extinction risk, an additional trait non included in the trait space, etc.) is distributed within the functional trait space previously obtained with funspace().The resulting map of the target variable within the functional space can then be represented graphically as a heatmap using plot().All the modules can handle multiple groups.The package also includes a generic summary() function that can be used to print the output of funspace(), funspaceNull(), and funspaceGAM().Function descriptions are reported in Table 1.A full worked example is shown in Boxes 1-3.core function of the package.By default, the function needs only one object to run: either an ordination object (either a PCA obtained using the base::princomp() function, or a PCoA obtained using the vegan::capscale() function, or a NMDS obtained using the vegan::monoMDS() or vegan::metaMDS() functions), or a TPDs object obtained with the TPD::TPDs() or the TPD::TPDsMean() functions, or a data frame including species coordinates in any other type of ordination (as well as any multivariate trait data) (Figure 1).By default, funspace() uses the first two dimensions of the provided space, but the user can plot any pairs of dimensions by modifying the PCs argument.We recommend users to set the pairs of PCs based on the output of funspaceDim().With this input, funspace() estimates the probability of occurrence of trait combinations within the space defined by pairs of dimensions using kernel density estimation with unconstrained bandwidth selectors combining the functionalities available in the R packages ks (Duong, 2007) and TPD (Carmona et al., 2019).The grid size for computing multivariate kernel density estimates is defined by the n_divisions argument that sets the total number of divisions of each dimension or trait; since the resulting functional space is bidimensional, the total number of cells in which the functional trait space is divided is equal to n_divisions.Note that increasing n_divisions results in longer time for computing the functional trait space (this can be an issue particularly when the funspaceNull() is used), but it increases the smoothness of the space edges in subsequent plotting (see Module 3).When the input is an ordination object, funspace() uses by default the range of the selected axes to build the grid limits in which estimating the multivariate kernel density estimates, but user-defined ranges can be used for calculations as well.The probability threshold used for multivariate kernel density estimates (i.e. the boundaries of the trait probability density function) can be set using the threshold argument (default = 0.999 corresponding 99.9% of the trait probability density function being retained).If the dataset includes a categorical variable (for example containing the family to which each species belongs), the same procedure can be applied to each level of the grouping variable by specifying it in the group.vecargument.In this case, by default, funspace() constraints group bandwidths for plotting to the global bandwidth, to avoid kernel density estimates of each group to exceed the borders of the functional trait space defined using the whole dataset.However, group-specific bandwidths can be generated by specifying fixed.bw= FALSE.
Other than the information that is later used for plotting, and when the input data is a PCA object, funspace() also returns a table reporting the loadings of the traits in the provided PCA (eigenvectors multiplied by the square root of eigenvalues) (Table 2).In this sense, it is important to note that the 'loadings' returned by base::princomp() are actually eigenvectors, which do not provide information about the amount of variance contained within each PC.The squared loadings are used to estimate the proportion of the original variance of each trait that is explained by each selected component (specified in the PCs argument) as well as by the two selected components (Table 2, Box 1).This information can help understanding how well each trait is represented by the selected functional space.In addition, and independently of the data input (ordination object, TPDs object, or matrix), funspace() always returns a table reporting the functional richness and functional divergence of the dataset (Mason et al., 2005) (Table 3, Box 1).Functional richness represents the amount of functional trait space (i.e. the area in a two-dimensional context) occupied by the set of species.Functional divergence quantifies how much the TPD function is distributed towards the edges of the functional trait space occupied by an assemblage.These indexes are calculated using the trait probability density approach (Carmona et al., 2016(Carmona et al., , 2019) ) across the whole dataset, and, if specified, per each level of a grouping variable.The user can retrieve all the funspace() output(s) using the summary() function (Box 1).
Finally, the function funspaceNull() allows the user to test how the area (i.e.functional richness) of the observed functional trait space differs from that of a null model.Two null models are implemented.(i) TA B L E 1 funspace functions by module.

Module Function Description
Building & exploring funspaceDim() Identifies the number of dimensions that are needed to build a functional trait space.

impute()
Imputes missing data in a trait matrix using trait-trait relationships and, if specified, phylogenetic information.
funspace() Builds a functional trait space including multivariate kernel density at different quantiles.
Multi-group option available.

funspaceNull()
Tests for statistical difference between the observed functional trait space area versus theoretical areas generated using either a multivariate normal or a uniform distribution.
Mapping funspaceGAM() Fits GAMs using a target variable as the response variable and functional trait space axes as a two-dimensional smoother predictor.Returns GAM model predictions within a functional trait space.Multi-group option available.
Plotting plot() If a funspace() object is passed as argument, the function prints a functional trait space with multivariate kernel density estimates.If a funspaceGAM() object is passed as argument, the function returns a GAM map (i.e.heatmap) of how a target variable distributes within the functional trait space.Multi-group option available.
Building & exploring and mapping summary() Gives a summary of the main features and results of objects created with the funspace(), funspaceNull(), and funspaceGAM() functions.

BOX 1 Using funspace module 1
We use one of the funspace example dataset to illustrate the package workflow.The dataset is a subset of Carmona, Bueno, et al. (2021) dataset, containing trait information for the traits defining the Global Spectrum of Plant Form and Function (GSPFF, Díaz et al., 2016, namely plant height, seed mass, specific stem density, leaf area, leaf nitrogen content on a mass basis and specific leaf area) for >10,000 vascular plant species.Trait data was compiled from the TRY database (Kattge et al., 2020).
# Load the package library(funspace) # Load example data (traits already log10-transformed and scaled) data.aux<− funspace::GSPFF_missing The dataset contains missing trait information, so we use the impute() function, coupled with phylogenetic information to impute missing trait data.The phylogenetic tree was retrieved using the V.Phylomaler R package (Jin & Qian, 2019) and is contained in the phylo object that is loaded with funspace.Note that, when using user-specified data, the trait dataset provided must include species names as row names; the same naming convention should be used in the phylogenetic tree (if provided).
# Load the example tree (it must be an object of class phylo) funspace example datasets also include a data frame with taxonomic information for the species in our GSPFF subset.We will use this information to illustrate how to deal with groups when using funspace Module 1.

PCA_subset <− princomp(GSPFF_subset)
A grouping variable can be specified in funspace() by passing a vector including the variable to the group.vecargument.
# Building the functional trait space (using the first two PCs) including groups The first null model generates data following a bivariate normal distribution (null.distribution= 'multnorm') with the same mean and variance-covariance structure as the original data.This null model corresponds to the expectation that some trait combinations, specifically those towards the centre of the space, are more likely than others, and the resulting space has the approximate shape of an ellipse (Carmona, Bueno, et al., 2021).(ii) Another null model generates a dataset with variables following a uniform distribution (null.distribution= 'uniform'), creating a functional trait space where all trait combinations within the range of the observed functional trait space axes are equally possible, and the space assumes the approximate shape of a rectangle (Díaz et al., 2016).Given a funspace() object as input, funspa-ceNull() generates a vector of null model areas across a user-defined number of iterations.Then, using the function as.randtest from the ade4 R package (Dray & Dufour, 2007), it tests for the statistical difference between the observed value of functional richness and the null expectation (considering in both cases the probability threshold value used in the funspace() object input).Hypothesis testing options follow the as.randtest() specifications.To get meaningful comparisons with the original functional trait space, funspaceNull() builds the trait probability density functions at each iteration using the same grid size and bandwidth for computing kernel density estimates as the funspace() object.funspaceNull() returns: the observed functional richness, the average functional richness of the null model across iterations, and the p-value and standardized effect size associated to the hypothesis test.This information can be retrieved using the summary() function.Note that funspaceNull() is meant to test the global space against null models, so it does not handle multiple groups.

| MODULE 2: MAPPING A FUN C TI ONAL TR AIT S PACE
funspace includes the function funspaceGAM() that automatizes GAM modelling steps needed for mapping patterns of a target summary(trait_space_global) summary(trait_space_families) In this example, the summary() function applied to the object trait_space_global returns the proportion of variance explained for each trait by each of the selected components of the PCA (Table 2) and the functional diversity indexes describing the global space (Table 3).In the multi-group case (i.e.trait_space_families object), functional indexes are returned per each level of the grouping variable as well (Table 3).Note that functional diversity indexes (for single or multiple groups) are returned independently of the funspace() input.
BOX 1 (Continued) From this output we can see that the first principal component is mainly explaining plant height, specific stem density and seed mass, whereas the second component is more strongly related to leaf traits (leaf area, leaf nitrogen content and specific leaf area).Finally, the Overall _ explained column reports the percentage of variance explained by the considered space (i.e. the sum of the variance explained by individual components).In this functional space, the quality of the representation of plant height, specific stem density and specific leaf area is better than that of leaf area and leaf nitrogen content.Note that the function output is slightly different since the table has been reorganized for presentation purposes.).The quantile threshold at which the indexes are calculated is also returned in the output.The table is returned as part of the summary() output.In this example we can see that functional richness and functional divergence are larger for the Fabaceae family than for the other families, reflecting the fact that legumes include both herbaceous and tree species (see Figure 3), whereas species within the other families are generally more similar between them.Note that the function output is slightly different since the table has been reorganized for presentation purposes.

Set of species
dimension within the functional trait space (Table 1, Figure 1).An example of how to use funspaceGAM() is reported in Box 2.
funspaceGAM() fits a GAM with a single response variable (defined by the y argument) and a two-dimensional smoother defined by the functional trait space axes as the bivariate predictor.GAM specifications are set to the default ones (see mgcv::gam() function, Wood, 2017), only the family argument can be specified by the user.Functional trait space axes are inherited by the funspace() object specified in the funspace argument of funspaceGAM().Given a funspace() object, the function automatically defines the functional trait space boundaries and grid size for generating GAM model predictions by inheriting the threshold and the grid size information from the funspace() object.Note that, the greater the grid size is, the longer the GAM computational time will be.
If a grouping variable is specified when running funspace() (see Module 1, Box 1), GAMs are also fitted per each level of the grouping variable and only within the portion of the space occupied by the data points of each subgroup (see details of fixed_bw argument of funspace() for different ways of estimating the bandwidth for each group).To avoid fitting models with not enough data, if there are fewer observations in a particular group than specified in the argument minObs (defaults to 30), the model is not fitted for that group, and a warning message is returned to inform the user.GAM summary statistics for all data pooled or, if specified, for each level of the specified grouping variable, can be retrieved using the summary() function.The output is not shown in the example below because the summary is equivalent to that of the function mgcv::summary().

| MODULE 3: PLOT TING A FUN C TIONAL TR AIT S PACE
funspace includes a generic plot() function that allows the user to plot objects of the funspace class created with the funspace() and funspaceGAM() functions (see Box 3 for a worked example).
By default, the only argument needed for plotting is an object of class funspace that is assigned to the x argument of the plot() function.If the plotting input was created with the funspace() function, the plotting function returns a functional trait space where coloured areas represent the density of occupation of the functional space by the considered species (i.e. the trait probability density function as described in Carmona et al., 2016Carmona et al., , 2019; Figure 2a).
When the input is the result of the funspaceGAM() function, the plot() function prints a heatmap in the functional space depicting the predicted values for the response variable within the external boundaries of the trait probability density function (Figure 2b).Note that, in both cases, the resulting plot(s) might appear more or less smooth depending on the number of n_divisions set when building the functional trait space using the funspace() function (Box 1).Larger values of n_divisions will result in smoother plots, but also in higher computational times.An important argument of the plot() function is the type argument.By default, this argument is set to 'global', and it allows the user to plot the functional trait space for all data pooled (Figure 2a,b).However, if type is set to 'groups', and given that a grouping vector has been already specified when building the funspace() object (Box 1), the plot() function creates a plot for each level of the grouping variable (Figure 3, see the fixed_bw argument of funspace()).
Many graphical parameters can be set for plotting.For example, if the input object comes from funspace(), when the argument quant.plot is set to TRUE (default is FALSE), contour lines are drawn at quantiles of the trait probability density function specified by the argument quant (by default these are the threshold argument that was used in funspace() and the 0.50 and 0.25 quantiles).If quant.
plot is set to TRUE when the input comes from the funspaceGAM() function, then the contour lines indicate the quantiles of the values of the response variable predicted by the GAM model.When quantile lines are displayed, their features (e.g.type, colour) can be set using specific graphical parameters (e.g.quant.lty,quant.col).
Independently of the input object, the gradient palette (with a number of colours specified by the ncolors argument, which defaults to 100) can be easily modified by specifying a vector including the two or more colours in the colours argument.Plot axes limits can be adjusted by modifying the base plot graphical parameters xlim and ylim.When the input used to create the funspace() object is a PCA, the user can also decide whether to plot arrows representing the trait loadings in the PCA by setting the argument arrows to TRUE

BOX 2 Using funspace mapping module
For this example, we will consider the GSPFF dataset that contains information on six traits for the set of species with complete trait information.We will create a response variable that increases as we move further from the origin of coordinates of the space and add some random normally distributed noise to avoid a perfect fitting.(default is FALSE).Full details on plotting options (e.g. points display, etc.) can be found in the funspace documentation.

| DISCUSS ION
funspace provides users with an important tool for conducting functional trait-space analyses from raw trait data and a way to increase the reproducibility of such analyses.The interconnected modules of the package provide in fact all the necessary steps to streamline sometimes troublesome operations, including imputation of missing data (Stewart et al., 2023) or the definition and analysis of functional spaces using multivariate kernel density (Mammola et al., 2021).
funspace is explicitly based on the TPD framework (Carmona et al., 2016(Carmona et al., , 2019)), which is not optimized for plotting.Thus, funspace extends the TPD functionalities in terms of plotting, and if the functions are defined within the same trait range (by setting the trait_range argument), the results of the two packages are fully comparable.This allows integrating funspace() outputs with additional analyses not included in the package, such as TPD-based dissimilarity, redundancy, or uniqueness (Carmona et al., 2016).In particular, the main features of funspace include: (i) it requires minimum input from users, increasing its scope of application; (ii) it is optimized for graphical representation and visual exploration of functional trait spaces.In the following, we discuss how funspace integrates and expands the current landscape of functional diversity R packages.

BOX 3 Using funspace plotting module
For this example, we use the trait_space_families and fit.gamobjects defined in Boxes 1 and 2 and generated using the funspace() and funspaceGAM() function, respectively.
# 1 -Plotting a functional space including a global (all data pooled) trait probability distribution.In this case, we will plot the trait_space_ families object.
The output is shown in Figure 2a.
# Plotting an object created with the funspaceGAM() function: we will use the fit.gamobject built in Box 2.
# To display the use of the 'type' argument, we also plot the functional trait space with an individual trait probability density function per each level of the grouping variable (Fabaceae, Lauraceae, Pinceae, Poaceae).The GAM output per each group is not shown, but it could be obtained replacing the trait_space_families object with the fit.gamone in the following code.The output is shown in Figure 3.
One advantage of funspace for users is its ability to generate functional spaces using only a matrix of functional features, eliminating the need for additional input data (i.e.abundance or presence/absence matrices) required by other packages that are mainly focused on community-level diversity indices.For example, both the BAT and betapart R packages (Baselga & Orme, 2012;Cardoso et al., 2015) require a community matrix to estimate the functional space.However, it should be noted that both packages are more focused on estimating alpha and beta diversity indices, including taxonomic, phylogenetic, and functional aspects, so their focus is different from that of our package.Moreover, and importantly, they do not include specific functions for visualizing or exploring functional spaces, nor any function to map external variables in such spaces.
Mapping different variables in a functional space is one of the main features of funspace, and this is done using the funspa-ceGAM() function.funspaceGAM() requires as input an object created by the core funspace() function and the target response variable, allowing for mapping target variables in any trait space in a seamless and intuitive manner.Finally, the resulting plots can be easily customized, and, by default, they are publication ready.To the best of our knowledge, only two R functions perform similar operations.These are the envfit() and ordisurf() functions included in the vegan R package (Oksanen et al., 2022) and they allow for fitting multiple regression or GAM with a similar purpose as funspaceGAM().However, envfit() and ordisurf() need the ordination axes as input.Thus, ordination axes need to be separately extracted from objects of ordination functions, requiring an additional step from the users.Moreover, those alternative functions, being more focused on modelling, provide default plots that are generally barebone, often requiring additional tuning from the users as well as additional packages to improve plotting (e.g.ggordiplots, Quensen et al., 2023).funspaceGAM() solves all these limitations by automatically retrieving ordination axes for modelling, thus providing a single-step solution for generating publication-ready plots, and associated statistical tests, of how target variables distribute within ordination analyses outputs.
Finally, the mFD R package (Magneville et al., 2022), a recently released extension of the previous FD package (Laliberté et al., 2014), offers a wide range of functions to calculate traitbased distances, construct multidimensional functional spaces, and compute various alpha and beta functional diversity indices.
In addition, mFD allows for the construction of functional spaces using only a trait matrix.However, and this is the main difference with funspace, mFD relies on convex hull estimation to build functional spaces, an approach that may have limitations with specific types of data.For example, convex hulls are extremely sensitive to outliers and are not able to account for the differential distribution of data points within trait spaces (i.e.multivariate The output shows that there are two hotspots (corresponding, from left to right, to herbaceous plants and angiosperm trees, see Díaz et al., 2016).The variance explained by each component and the loadings of the original traits are also shown.(b) Plot showing the GAM predicted values for a variable (see Box 3 for explanations) within the same functional space; the input was created using funspaceGAM().
The heatmap shows the predicted values for a variable with increasing values as we move away from the centre of coordinates of the functional space.Contour lines are the quantiles of the GAM predictions expressed in the same unit as the response variable.The code to reproduce each panel is available in Boxes 1-3.density) (Mammola et al., 2021).This limitation of convex hulls becomes especially important in cases when missing data need to be imputed (Stewart et al., 2023).Additionally, mFD lacks support for external variable mapping within the functional spaces, although this falls outside the package's intended scope.In sum, mFD and funspace are complementary packages to broaden the set of tools for users to explore the diversity of organismal form and function across the tree of life.

| CON CLUS IONS
funspace streamlines all the necessary steps to build, explore, map, and plot functional trait spaces, making such analyses readily accessible to any user interested in bivariate or multivariate functional trait analyses.Importantly, funspace automatizes most of the procedures needed to run such analyses, and for this reason can be of use to inexperienced R users as well.This package provides a F I G U R E 3 Plot of a funspace object created with the function funspace() including groups for four major plant families (Fabaceae, Lauraceae, Pinaceae, Poaceae).For each family, its corresponding trait probability distribution is represented within the functional space; the 99.9% probability quantile of the global trait probability distribution (i.e. the one including all species from the four families together) is shown to provide a common reference and make comparisons easier.Colours, interpretation of contour lines and explanation of axes is the same as in Figure 2a.The code to reproduce this figure is available in Boxes 1 and 3.

F
I G U R E 2 plot() output by input object.(a) Plot of an object created with the function funspace() obtained using all data points in the GSPFF dataset.Colours indicate the probabilistic distribution of trait combinations in the functional trait space defined by a PCA (red = high probability; yellow = low probability).Contour lines indicate 0.99, 0. 50, and 0.25 quantiles of the probability distribution.
TA B L E 2Note: A table summarizing the percentage of variance of each trait that is explained by each principal component (Comp.1 and Comp.2 in this example) and across components (Overall_explained).The table in this example is returned as part of the output of summary(trait _ space _ global), corresponding to the global spectrum of plant form and function using imputed trait information (see Module 1 and Box 1).