A kernel integral method to remove biases in estimating trait turnover

Trait diversity, including trait turnover, that differentiates the roles of species and communities according to their functions, is a fundamental component of biodiversity. Accurately capturing trait diversity is crucial to better understand and predict community assembly, as well as the consequences of global change on community resilience. Existing methods to compute trait turnover have limitations. Trait space approaches based on minimum convex polygons only consider species with extreme trait values. Tree‐based approaches using dendrograms consider all species but distort trait distance between species. More recent trait space methods using complex polytopes try to harmonise the advantages of both methods, but their current implementation has mathematical flaws. We propose a new kernel integral method (KIM) to compute trait turnover, based on the integration of kernel density estimators (KDEs) rather than using polytopes. We explore how this approach and the computational aspects of the KDE computation can influence the estimates of trait turnover. The novel method is compared with existing ones using justified theoretical expectations for a large number of simulations in which the number of species and the distribution of their traits is controlled for. The practical application of KIM is then demonstrated using data on plant species introduced to the Pacific Islands of French Polynesia. Analyses on simulated data show that KIM generates results better aligned with theoretical expectations than other methods and is less sensitive to the total number of species. Analyses for French Polynesia data also show that different methods can lead to different conclusions about trait turnover and that the choice of method should be carefully considered based on the research question. The mathematical properties of methods for computing trait turnover are crucial to consider because they can have important effects on the results, and therefore lead to different conclusions. The novel KIM method provided here generates values that better reflect the distribution of species in trait space than other methods. We therefore recommend using KIM in studies on trait turnover. In contrast, tree‐based approaches should be kept for phylogenetic diversity, as phylogenetic trees will then reflect the speciation process.


| INTRODUC TI ON
Biodiversity is a complex concept and can most easily be quantified by distinguishing three complementary and interrelated facets: taxonomic diversity based on a site-by-species matrix that captures the compositional properties of a community; phylogenetic diversity that captures the evolutionary relatedness among community members, using phylogenetic distance between species alongside the site-by-species matrix; and trait diversity that describes a community according to the traits of its resident species, using a species-by-trait matrix alongside the site-by-species matrix (Devictor et al., 2010).The study of functional traits has been advocated as fundamental to better understand and quantify community assembly (McGill et al., 2006), as well as the impact of global change on community resilience and on the ecosystem services that biodiversity provides (Gross et al., 2017).For example, through comparison with null models and by relating traits to environmental gradients and to each other, trait diversity can provide information about the assembly processes structuring an ecological community (Ackerly & Cornwell, 2007), including biotic interactions between species (Laureto et al., 2015).It also enables the estimation of components of ecosystem function, such as nutrient use and storage, or ecosystem productivity (Cadotte et al., 2011;Hillebrand & Matthiessen, 2009).
In addition to the decomposition of biodiversity into taxonomic, trait and phylogenetic components, unravelling how biodiversity is organised requires an understanding of how assemblages of species are more or less similar to one another at different places and times, that is, turnover (Anderson et al., 2011).To do so, beta (β) diversity provides a direct link between biodiversity at the regional (gamma-γ-diversity) and local (alpha-α-diversity) scales (Anderson et al., 2011;Chao et al., 2005Chao et al., , 2019)).In particular, taxonomic β diversity has been shown to be important for assessing the effects of conservation actions (Socolar et al., 2016), for example, for estimating the effect of the spatial distribution of protected areas and their subdivision into multiple subareas on species diversity (Deane et al., 2022), or for extrapolating regional species richness from limited data (Kunin et al., 2018).Although using β diversity to describe trait turnover has received less attention than to describe species turnover (taxonomic β diversity), it has been applied to measure change in trait diversity across communities and regions (e.g.Carmona et al., 2012;Loiseau et al., 2017;Siefert et al., 2013;Swenson et al., 2012;Villéger et al., 2013).
As a valuable and increasingly measured biodiversity facet, there are multiple important steps to consider when estimating trait turnover over space or time.First, the precise choice of traits can substantially influence the outcome (Petchey & Gaston, 2006).
Second, despite recent initiatives to collate large amounts of data for multiple traits across species (e.g.Kattge et al., 2020;Middleton-Welling et al., 2020;Tobias et al., 2022), trait data are still missing for many species and types of traits across taxonomic groups.Finally, and also the focus of this work, different mathematical methods exist to compute trait diversity and turnover that differ in outcome and therefore in the conclusions drawn about biodiversity (Loiseau et al., 2017;Sobral et al., 2016;Villéger et al., 2017).A systematic comparison of these methods can help identify informative and robust methods and establish standards for quantifying trait turnover.
There are two main categories of methods for calculating trait β diversity: (i) methods based on the concept of trait space (referred to here as the 'trait space approach', and (ii) methods that use dendrograms (referred to here as the 'tree-based approach').The trait space approach is based on a multidimensional space whose axes are determined by the traits included in the analyses.Axes can correspond directly to the original traits or can be derived from these traits through ordinations to reduce dimensionality.A particular species or population is typically represented as a single point in this trait space, which represents the average trait value over the species or the specific population considered, and a polytope is computed as the trait envelope of the set of points representing the species present in a community or assemblage.The minimum convex polytope (MCP), a convex hull, that encompasses all species of a community in the trait space (Figure 1), has traditionally been used in these analyses (Loiseau et al., 2017).As the MCP only captures information about the species with extreme trait values in a community, it is sensitive to outliers and ignores how species are distributed in the trait space, which can be crucial to delineate the functional roles of different species within an ecosystem (Mouillot et al., 2021).Although other hull methods can be used to compute the trait envelope (e.g.Irl et al., 2017), they are typically computationally intensive and have seldom been applied to β diversity analyses.
The tree-based approach consists of computing all pairwise distances between species based on a set of traits, typically using the Gower distance to incorporate both continuous and discrete traits.
As for the trait space approaches, a trait value can represent the average value over the whole species or a specific population.A clustering algorithm is then applied to these distances to generate a dendrogram, from which measures of β diversity can be computed.
Although the tree-based approach considers all species in the computation of trait turnover, the dendrogram splits into successive branches, and using the length of the branches connecting two species as a measure of distance distorts the original trait distance In contrast, tree-based approaches should be kept for phylogenetic diversity, as phylogenetic trees will then reflect the speciation process.

K E Y W O R D S
beta diversity, convex hull, French Polynesia, hypervolume, kernel, trait turnover, traits between them compared to the distance obtained through ordination in the trait space (Maire et al., 2015).In addition, the choice of the clustering algorithm for generating the dendrogram will inevitably influence the outcome (Loiseau et al., 2017).
The convex hull of trait space and the tree-based approach therefore make different computational trade-offs, and the appropriateness of the two approaches for measuring trait β diversity has been debated (Loiseau et al., 2017).In response to this debate and to incorporate information from all species, Mammola and Cardoso (2020) proposed another trait space approach where polytopes are defined by applying a threshold to the kernel density estimation (KDE; Figure 1; see details in Methods below).The resulting polytope is typically not convex, and its shape better reflects the distribution of species in the trait space.Although it has the potential to provide a more accurate estimate of trait diversity than the other two methods, this has not been formally assessed.The computational aspects when computing kernel densities have also largely been overlooked.As we demonstrate here, the appropriate contribution of all species in the community to β diversity is key to a robust measure of trait turnover.
Here, we propose a new trait space method, the kernel integral method (KIM), for computing trait β diversity that is based directly on the integration of the KDE rather than on the polytope.How the computational aspects of the KDE can influence the estimates of trait β diversity with different methods is examined, as well as how systematic bias can arise in the computation of trait turnover due to species distributions in the trait space of a community.A set of theoretical examples is used to compare methods, for which the behaviour of the trait β diversity metric can be justified, and therefore, biases can be identified.We further apply the KIM method to compute non-native plant trait turnover across islands and archipelagos of the Pacific Islands of French Polynesia and compare results with the other methods, to show how the theoretical differences between methods unveiled through the simulations can generate different results and lead to different conclusions when analysing real data.

| Convex hull
Computing trait turnover between two communities using the convex hull methods simply consists in computing (i) the minimum convex polytopes (MCP) encompassing the points representing average trait values for a species or a specific population for each community, and (ii) the hypervolumes of the intersection and the union of these two MCPs (Figure 1a,e).It is then possible to compute a range of β diversity indices based on these four values.Here, following Mammola and Cardoso (2020), we used the Jaccard dissimilarity index J (Jaccard, 1908) and the Williams replacement index W (Williams, 1996), defined as: F I G U R E 1 Summary of the trait space approaches for two communities with different species in a two-dimensional trait space.(a, e) The convex hull remains the same irrespective of the additional species present in community 2, resulting in the same outcome when computing β diversity metrics.(b, f) The KDH (kernel density hypervolume) method generates a polytope for each community, whose shape will vary with the absence or presence of the additional species in community 2 and is often non-convex.As a result, the outcome of the Jaccard dissimilarity or the Williams replacement formulas will differ.(c, d) KDEs (kernel density estimators) corresponding to the polytopes in (b).(g, h) KDEs corresponding to the polytopes in (f).
where A and B are the MCPs of the two communities, and V(A) and V(B) their hypervolumes.V(A∩B) and V(A∪B) are therefore the hypervolumes of the intersection and union of the two MCPs.In the analyses, the MCPs and the indices were computed using the hull.build()and hull.beta() functions from the BAT R package V.2.8.1 (Cardoso et al., 2015(Cardoso et al., , 2022)).The Williams replacement index evaluates the contribution of trait replacement to trait β diversity (Carvalho et al., 2012(Carvalho et al., , 2013)), and the difference between Jaccard and Williams indices quantifies how the trait richness difference between communities contributes to β diversity.Although there are other approaches and indices that can decompose β diversity into turnover replacement components, the relevance of these approaches is still debated (Baselga, 2010;Baselga & Leprieur, 2015;Cardoso et al., 2014;Carvalho et al., 2012Carvalho et al., , 2013)).
This debate is beyond the scope of this manuscript, and, to compare methods the decomposition used by Mammola & Cardoso (2020) was followed (see sections on kernel density hypervolumes [KDH] below), readily available from the BAT R package (Cardoso et al., 2015(Cardoso et al., , 2022)).
The main issue with the convex hull method is that it is insensitive to the addition or removal of species within the MCP in the trait space (Figure 1).A corollary is that it is sensitive to outliers, as they will define the MCP.The convex hull method is therefore highly biased in the sense that the different indices computed with this method are entirely determined by the species defining the extreme trait values in the community.

| Tree-based method
The tree-based method consists in computing a dendrogram from the trait distance between all species in the species pool (i.e. the entire list of species over all included sites, not just those occurring in the pair of sites for each calculation of trait turnover; Figure 2a-note that as for the trait space approaches, one may decide to use trait values at the population rather than species level).Multiple clustering algorithms can be used to generate the dendrogram.Here, the approach in Loiseau et al. (2017) was followed using the unweighted pair group method with arithmetic mean (UPGMA) algorithm, and the hclust() function from the stats R package (R Core Team, 2022).
This has been shown to best conserve distances between species compared with the original distances in the trait space.
For each site, a subtree including only the species present is generated by trimming the overall tree (Figure 2b,c).It is then possible to compute the trees corresponding to the union and the intersection of the two subtrees (Figure 2d,e).Equations 1 and 2 can then be adapted to compute the Jaccard and Williams indices, by using the (1) total length of remaining branches.Importantly, the subtrees must be computed from the original tree generated from the entire species pool, not those computed from only the residing species.This conserves the internal branches in the union and the intersection of the two subtrees, even if these internal branches do not lead to any present species (see branch h in Figure 2e).Therefore, for the example of Figure 2, Equations 1 and 2 become: where A and B are the trees computed for the two communities, and L(A) and L(B) the sum of all their branches.L(A∩B) and L(A∪B) are therefore the length of the branches of the intersection and union of the two trees.The equations including lower-case letters (a, …, h) show how to apply this formula for the trees presented in Figure 2.
The tree-based method offers the advantage over the convex hull method that all species will be accounted for when computing the β diversity indices, therefore removing the bias towards species with extreme trait values.However, the clustering algorithms often generate branch lengths between species in the dendrogram that differ from the original distances in the trait space, which will necessarily influence the value of any β diversity index (Loiseau et al., 2017).

| Kernel density hypervolumes
Mammola and Cardoso (2020) introduced the use of KDH for computing indices of species turnover.This approach generates polytopes from KDEs.The polytopes are often nonconvex (and can even be disjunct) and can be seen as a trait envelope around the species points in the trait space.The recommended method is based on a Gaussian estimator of the KDE (Mammola & Cardoso, 2020) and follows a series of four steps (see Blonder et al., 2018 for further details): (i) Points are drawn randomly within a hypersphere around each species point in the trait space (as for the convex hull approach, each point can also represent average trait values at the population level if data are available); (ii) these points are resampled to uniform density; (iii) a KDE is computed from these points; (iv) a threshold (typically 95%) is applied to truncate the KDE and define the polytope, from which hypervolumes can be computed.The indices are then computed as per Equations 1 and 2, where A and B represent polytopes instead of MCPs.Here, the kernel.beta() function from the BAT R package V.2.8.1 (Cardoso et al., 2015(Cardoso et al., , 2022) ) was used to apply the KDH method.
This method, although more computationally intensive than the convex hull method, allows the distribution of species points in the trait space to contribute to defining the polytopes and therefore the hypervolumes used in the computation of the turnover indices (Figure 1b-d,f-h).As a result, the KDH method is less sensitive to outliers.
The KDH method nonetheless has some caveats.First, the choice of the threshold used to construct the polytope will necessarily influence the components of the β diversity indices, and therefore the final values.Second, by resampling random points to uniform density, some information about the distribution of species points in the trait space is lost.More specifically, a bias is introduced by giving a lower weight to species with traits similar to each other, for which the density of random points is higher before resampling.Finally, in the current implementation of the method in the BAT R package V.2.8.1 (Cardoso et al., 2015(Cardoso et al., , 2022)), the radius of the hyperspheres within which random points are drawn around the species points and the bandwidth used during the computation of the KDE (the bandwidth is a parameter that determines how smooth the KDE will be) are determined based on the species point distribution of each community separately, using the estimate_bandwidth() function from the hypervolume R package.Following this implementation, the more similar species are to each other, the smaller the volume of a hypersphere and the closer random points will be to each other, and consequently the KDE will show a steeper gradient (Figures B1-B24, compare columns 3 and 4 with columns 6 and 7).In other words, a systematic bias is introduced when using community-specific bandwidths for assessing trait turnover between two communities, due to the difference of how tight species are packed in the trait space of each community.For a β diversity index to be unbiased we argue that all species should have the same weight and contribute equally to the KDE when relative abundance and intraspecific trait variation are not concerned.

| A kernel integral method
To solve the issues associated with the KDH method, we propose a novel computational method to compute trait turnover in the trait space.This KIM method computes β diversity indices from the kernels themselves, therefore removing the influence of the threshold used to generate the polytopes, and uses different kernels than those used in the KDH method.The KIM method consists of using only steps (i) and (iii) from the KDH method.(i) Points (typically 1000, but the number can be adjusted to account for species abundance, for example; see Appendix E) are drawn randomly within a hyperellipsoid-an ellipsoid in n dimensions-around each species point in the trait space.Under a simple setting, the diameters of the hyperellipsoid can be the same for all species along different trait axes, which would assume that all species or populations have the same trait variability, and that the variability is the same for all traits.As species or populations can have different trait variability (i.e.intraspecific trait variability can differ between species or populations), one can define species-specific diameters.Similarly, different traits can have a different variability for the same species, resulting in elongated hyperellipsoids (see Appendix E for illustrations).In the following, we used the same intraspecific trait variability for all traits and all species for simplicity, and hyperellipsoids are actually hyperspheres.(iii) A KDE is computed from these points, and rescaled between [0,1]. (3) From the KDE, we then propose the following equations to compute the Jaccard dissimilarity index and the Williams replacement index: where KDE A and KDE B are the KDEs for communities A and B, and ∫ KDE A is the integral of the KDE for community A over all dimensions of the trait space.This is similar in essence to the index of niche overlap proposed by Mouillot et al. (2005).In practice, since KDEs are computed as multidimensional matrices, an integral is simply computed as the sum of all elements of the matrix.The minimum and the maximum of two KDEs are analogous to the intersection and the union of the polytope in the KDH method (Figure 3).
Kernel integration method overcomes the limitations and biases of the KDH method.First, there is no need to define a threshold: If the KDE is estimated over a large enough area or volume, the local kernel density will approach zero and the integral will therefore converge.Second, the radius within which the random points are drawn is the same for all communities (but a suitable value must be chosen, which can be adjusted to account for intraspecific trait variability).
Finally, because there is no resampling to uniform density, the distribution of species points in the trait space will be reflected more accurately in the KDE.

| Test of the different methods on theoretical data and expectations
The theoretical advantages and caveats of each of the different methods are described above.Next, seven different methods were applied to explore how the differences between methods influence the results (Table 1): 1.A convex hull method (hereafter COVHULL).

A tree-based method (hereafter TREE).
3. The original kernel density hypervolume method with communityspecific bandwidths and uniform resampling (hereafter KDH V1). 4. A modified kernel density hypervolume method computed with the same bandwidth for each pair of communities and uniform resampling (hereafter KDH V2), to explore the influence of the bandwidth on the outcome.5.The kernel integral method using kernel densities estimated with community-specific bandwidths and uniform resampling (hereafter KIM V1), to explore the influence of the kernel-based vs the polytope-based formulas (Figure 3).
6.The kernel integral method using kernel densities estimated with the same bandwidth for both communities within a pair and uniform resampling (hereafter KIM V2), to further explore the influence of the kernel-based vs the polytope-based formulas.
7. The kernel integral method estimated with the same bandwidth for both communities within a pair and without uniform resampling (hereafter KIM V3).
For each method, Jaccard dissimilarity and Williams replacement were computed, as defined in Equations 1-6.The seven methods were then examined for how they behaved in a set of theoretical contexts for which qualitative predictions could be made for how an index of turnover should behave when capturing the trait differences between communities.
In total, 72 different pairs of communities were simulated (Figure 4 and Figures A1, B1-B24) and the 14 indices computed (the Jaccard and Williams indices for each of the seven methods) for each pair.For simplicity and computational efficiency, a trait space defined by two theoretical traits was used.Each community was first delimited by a MCP represented by four species Computation of trait turnover for two pairs of communities (one pair in (a) and one in (b)) following different approaches.In each graph (a, b), the curves fictional KDEs (kernel density estimators) for the two communities, in one dimension (for simplification, we assume their densities are 0 beyond intersecting the horizontal axis).Using the KIM (kernel integration method) formula, The Jaccard index is computed as one minus the area in grey divided by the striped area (Equation 5), and the value is different for the two pairs of communities in (a) and (b) (as is the value of Willams replacement index, Equation 6).(c) The horizontal segments represent one-dimensional polytopes (defined using the values where the KDEs intersect the horizontal axis for simplification), used to compute the Jaccard or Williams indices using the KDH (kernel density hypervolume) method (Equations 1 and 2).Contrary to KIM, the KDH method only generates a single value for each index for the two pairs of communities.
TA B L E 1 Summary of the methods used to compute trait turnover.There is one obvious difference between these theoretical communities and communities that would be analysed for real-world applications.Real-world communities belonging to the same ecological system (such as those described in the next section) will usually share species, resulting in many species points overlapping in the trait space (unless the approach is applied at the population rather than the species level, that is, each community is represented by its specific population).Here, we used independent random species distributions in the trait space for the two communities to have greater flexibility over the species distributions, to explore in detail how each of the seven methods would behave across a wide variety F I G U R E 4 One of the 50 instances of the nine theoretical configurations of pairs of communities for different sizes of the MCPs, defined by four extreme species within each community (i.e.corners of the squares, see Figure A1 for same size MCPs), using only 10 randomly drawn species points for clarity (point distributions for 40, 70 and 100 species were also generated).How the MCPs overlap ("nested", "overlapping" or "disjunct"), and how the rest of the species points within the community are distributed within the MCPs was varied ("different"-distributed in opposite corners of the MCPs-, "random"-randomly distributed within the MCPs-or "similar"-distributed either within the same small area, or in the closest corners of the MCPs).When the trait distributions of species within the MCPs follow a "different" or "similar" configuration, three out of the four extreme species represent outliers within the communities.For the "different" configuration, the outlier species are more similar to each other than the rest of the species between the two communities.For the "similar" configuration, the outlier species are more different to each other than the rest of the species between the two communities.

Method
of extreme configurations, and to better disentangle the implications of their computational specificities.F I G U R E 5 Differences in Jaccard dissimilarity between the seven methods summarised in Table 1, for MCPs of different sizes, for 10 and 100 species points (see Appendix C for the full set of results).The first row shows the qualitative differences in β values for trait turnover predicted under different simulated configurations of species points in the trait space when the MCPs of the two communities have different sizes.The position of the symbols on the y-axes are coarse approximations of the value an index should take for these configurations.Overall, β values, are expected to decrease as the overlap between the trait profiles of the two simulated communities increases.Jaccard dissimilarity accounts for differences in the spread of the trait profiles in the trait space (i.e.trait richness).For detailed explanations about the changes in β values, see Table A1 in Supporting Information.
The replacement component should decrease if the difference in area covered by the two sets of species points increases.These patterns should be especially clear for large numbers of species, that is, for high densities of species points.For few species and low species point density (e.g. when only 10 species points were randomly drawn in the trait space), these patterns are expected to be weak, because the stochastic element of the species point distributions may obscure the results.

F I G U R E 6
Changes in Williams replacement index for the seven methods summarised in Table 1, for MCPs of different sizes, for 10 and 100 species points (see Appendix C for the full set of results).The first row shows the qualitative differences in β values for trait turnover predicted under different simulated configurations of species points in the trait space when the MCPs of the two communities have different sizes.The position of the symbols on the y-axes are coarse approximations of the value an index should take for these configurations.Overall, β values, are expected to decrease as the overlap between the trait profiles of the two simulated communities increases.Williams replacement is independent of differences in trait richness.For detailed explanations about the changes in β values, see Table A1 in Supporting Information.

| Established non-native plants in French Polynesia
The theoretical pairwise configurations presented above are used because it is possible to make some qualitative predictions about For these 417 species, data on species woodiness (woody vs. herbaceous species), seed mass, plant height and specific leaf area (SLA) were extracted from multiple trait databases, including TRY (Kattge et al., 2011(Kattge et al., , 2020)), LEDA (Kleyer et al., 2008), PLANTATT (Hill et al., 2004), Austraits (Falster et al., 2021), BIEN (Maitner, 2022), EcoFlora (Fitter & Peat, 1994) and BROT (Tavşanoğlu & Pausas, 2018).Seed mass, plant height and SLA have been used to characterise different plant life strategies (Díaz et al., 2016;Westoby, 1998).When different databases contained different values, means were used for seed mass, plant height and SLA, and the most frequent category for woodiness.Data on plant woodiness were available for all 417 species.Trait data for seed mass and plant height were only available for 250 out of the 417 species.Data for seed mass, plant height and SLA were only available for 124 out of 417 species.Three sets of analyses were run: (i) a set for the 250 species with data on seed mass and plant height, (ii) a set for the 124 species with data on the three traits, and (iii) a set for the same 124 species, using data on seed mass and plant height only, to assess the robustness of the results to data availability and trait selection.The results are presented and discussed mainly for seed mass and plant height for the 250 species (see Figure D1 for the distribution of plant species in this two-dimensional trait space), as these represent more than half of the species and results should thus be less biased despite using only two traits.
Data from multiple open-access databases were used due to the lack of trait data that would be specific to the populations introduced to French Polynesia and trait values may be specific to the populations on which they were estimated.We therefore implicitly assumed that averaging trait values across databases would be representative of values at the species level, and that trait average was the same for the French Polynesia populations and for the whole species.Similarly, in the absence of data on trait variability, we considered that all traits had the same variability for all species, represented by hyperspheres with the same radii in the analyses.
Prior to analysis, seed mass, plant height and SLA were log-transformed and rescaled between [0,1] so that the traits would be more uniformly distributed in the trait space.Jaccard dissimilarity and the Williams replacement indices were then computed for all species together, and for woody and herbaceous species separately.
This split was made for a more comprehensive assessment of potential differences between methods, because woody and herbaceous plants are different in terms of growth form and function (Díaz et al., 2016).The indices were also computed for all French Polynesian islands, and for each archipelago separately.
Finally, the behaviours of the different indices were analysed using randomisation tests.The presence-absence matrices were randomised for all islands and for each archipelago by keeping species occupancy and island richness constant (i.e. the sim9 algorithm from [Gotelli, 2000]).The Jaccard dissimilarity and Williams replacement indices generated by the seven methods for the original matrices were compared with the indices computed over 10 randomised matrices for each original matrix (the number of randomisations was constrained by computation time).
The purpose of applying the methods to an empirical data set was specifically to examine how the biases and theoretical differences between the seven methods examined through simulations can generate different results and lead to different conclusions when analysing more complex, real data, and to assess each method's range of sensitivity.Because each archipelago contains multiple different islands whose combinations will fall across a large spectrum of trait profile configurations, it was not possible to define a priori expectations for this case study.The purpose is therefore not to determine which methods are in line or not with a priori expectations, contrary to the theoretical analyses, but rather to demonstrate the consequences of applying biased methods for research findings.

| Theoretical data
Overall, KDH and KIM tended to converge towards similar values and behaviours as the number of species points increased ), corresponding to theoretical expectations (Figures 5-7, Figures A2 and A3).In contrast, the convex hull and the tree-based methods generated indices of turnover different from both the other methods and from the theoretical expectations.
The main differences between observed and expected values for all methods, except the tree-based method, were for the contribution of replacement to overall turnover (computed as the ratio of the Williams replacement index to the Jaccard dissimilarity index), for MCPs of the same size in the nested / random and the nested / similar configurations, which was lower than the expected value of 1 (Figures A3 and C6).This is likely because these are the configurations for which the values of the Jaccard index are small, and small changes in Williams replacement index due to stochasticity in the distribution of species points in the trait space will be disproportionally large.
For all indices of turnover, the three KIM methods generated values above 0.5 and above other methods when the point distri- Adjusting the bandwidth to be common between communities in each pair in the computation of the kernels for the KDH and KIM methods (i.e.switching from V1 to V2) resulted in lower dissimilarity values, both for the Jaccard dissimilarity index and the Williams replacement index, for all configurations.This is because the radius of the hyperspheres and therefore the steepness of the kernels were the same for both communities in the V2 methods, increasing similarity.The effect of removing the resampling of random points to uniform density (i.e. from KIM V2 to KIM V3) often had a larger and more variable effect than adjusting the bandwidth.This is because the values generated by KIM V3 could be either larger, smaller or in between those of the KIM V1 and V2 methods.
Overall, results were consistent across the 50 replicates with little overlap between the interquartile intervals between methods (Figures C13-C20).Variability across replicates also increased as the number of species points decreased.KIM V3 tended to generate more variable values across replicates (but not for all configurations, see e.g. Figure C16), which is consistent with the fact that it was designed to capture more precisely the contribution of each species, and therefore differences in their distributions across replicates.

| Established non-native plants in French Polynesia
Raw values of Jaccard dissimilarity, of Williams replacement and of the contribution of replacement to turnover differed greatly between the different methods when applied to the case study data.
Maximum differences in values between methods were around 0.6 for Jaccard dissimilarity, 0.  D3).
Importantly, compared with the other methods, the KIM methods sometimes generated a different ranking between archipelagos for Jaccard dissimilarity.This is especially true for woody species, for which KIM V3 suggests that trait turnover was higher for the Gambier than for any other archipelagos, for all combinations of traits, and for both Jaccard and Williams replacement indices (Figure 8,Figures D2 and D3).By contrast, the other methods generated results that were more variable depending on the combination of traits and species used.
Randomisation of presence-absence matrices show that the KDH V2, KIM V2, KIM V3 and Tree methods tended to generate more consistent values for Jaccard dissimilarity compared with the convex hull, KDH V1 and KIM V1 methods (Figure D4).For Williams replacement, values were also more consistent across randomised matrices for KIM V3 than for the other methods, especially for herbaceous species.

| DISCUSS ION
Here, we compared existing and novel methods to compute trait turnover for simulated and empirical data, to illustrate how differences in the computational aspects of these methods reflect different aspects of trait diversity and can affect inferences made from trait diversity comparisons.

| Theoretical aspects of trait turnover computation
Comparing the seven methods using simulated data, for which we had complete control of the community trait profiles in the trait space, showed that the computational specifics of each method has significant effects on the value of trait β diversity.The two original approaches to assess trait turnover, the convex hull and tree-based methods, consistently generated results most different from theoretical expectations (Figures 5-7).Although using a trait space approach better conserves trait distance between species than a tree-based approach when all species are included (Maire et al., 2015), this property is broken when species are ignored in   The respective effects of these three computational aspects on trait turnover (β values) depend on the index used (Jaccard dissimilarity or Williams replacement) and on the configuration of the community trait profiles.For example, Jaccard dissimilarity is sensitive to the difference in bandwidth between communities (Figure 5).This is because using different bandwidths changes the shape of the KDEs (akin to making the distributions larger or narrower in Figure 3) and generates polytopes with different areas (akin to changing the lengths of A and B in Figure 3).Consequently, Jaccard dissimilarity values reflect this artificial difference in trait richness, whereas Williams replacement does not.In contrast, for Williams replacement, the use of polytopes or kernel integrals proved to be the most important factor (Figure 6).This is because kernel integrals better reflect small variations in the shape of the KDE (akin to changing the shape of the distributions and the overlapping area in Figure 3) and thus better capture trait replacement.Similarly, resampling also affected Williams replacement for communities with an 'overlapping' configuration (Figure 6), especially for species-poor communities.This is because, when compared to the more uniformly distributed trait profile of species-rich communities, each species has a greater effect on the shape of the trait profile in species-poor communities, and the idiosyncrasy in the position of different species can drastically change trait profiles if resampling is not applied.
Both trait envelope and kernel-based community trait profiles can have complementary uses, and the choice of an analytical approach will depend on the research or management question.On the one hand, species with extreme trait values defining a trait envelope for a given community can help capture the whole range of trait values of species that may potentially join the community.The trait envelope may therefore be an important piece of information to assess the risk of potential invaders to a region (e.g. to test the 'join the locals' vs. the 'try harder' hypotheses; Tecco et al., 2010), or an indicator of the loss of trait extremes.The CONVHULL method is appropriate for such applications.By contrast, capturing the trait distribution of all species in a community in the trait space provides a more comprehensive description of trait diversity and is necessary for identifying community assembly processes (Falster et al., 2017).
The distributional profile of species in trait space can also highlight gaps within the trait envelope, where introduced species with corresponding traits could have a better chance of establishing (i.e. the 'empty niche hypothesis ';MacArthur, 1970;Molofsky et al., 2022).
Changes in a community trait profile quantified in this way could therefore be used to reflect changes in ecosystem resilience or quantify trait redundancy (Hui et al., 2021;Mouillot et al., 2021).The results presented here demonstrate that the KIM V3 method is most informative and least biased for addressing trait diversity questions.

| Empirical test of methods using plant data from French Polynesia
The empirical testing of these methods provided further insight on the behaviour of the different methods for a mixture of configurations of species points in the trait space (Figure D1), and on how for KIM V3).These different conclusions, in addition to the different rankings generated by the different methods could be crucial for conservation decisions.For example, assuming management actions are influenced by species traits, low trait dissimilarity between islands, as indicated by KIM V3, would suggest that a similar management approach is appropriate across most islands, simplifying management and potentially improving management efficiency.In addition, the high contribution of replacement to turnover suggests that existing differences in community trait profiles are unlikely to be the result of differences in colonisation pressure, and may point towards either idiosyncratic or niche-driven factors.As only three traits were used in these analyses due to lack of data for other traits, further analyses would be needed to confirm these recommendations.In particular, other traits related to niche tolerance, competitive ability and dispersal could be important.
Factors to consider when selecting traits to include in trait turnover analyses are beyond the scope of this work, and frameworks exist for different systems and applications (e.g.Bremner et al., 2006;Lavorel & Garnier, 2002;Luck et al., 2012).It is nonetheless important to note that some factors have direct bearing on other potential sources of bias.(i) The integration and correlation between traits should be examined, not only for ecological and evolutionary reasons (Falster et al., 2017), but because non-independence of traits can bias the analysis (Pigliucci, 2003).(ii) Increasing the number of dimensions of the trait space (i.e.including many traits without using a reduction of dimensionality method) may lead to strange behaviours of the different indices due to differences in how the volumes of different shapes vary as the number of dimensions increases (Smith & Vamanamurthy, 1989).It is therefore advisable to apply a principal component analysis or similar method to restrict the analyses to three or less dimensions or syndromes (Salguero-Gómez et al., 2018).(iii) Trait distributions (mean and variance) may differ between multiple populations of the same species, and may not follow a unimodal distribution with homogenous variance for all traits.This would especially be likely for populations introduced in novel environments, such as the ones used in this study.However, although KIM allows to account for differences between populations by characterising each population by a species point and drawing the random point according to a distribution specific to each given population, such fine data at the population level is often not available, and in the absence of population-specific trait data, using a normal distribution centred around the trait mean at the species level can be used for the KDH and KIM methods.The radii of the hyperellipsoids within which random points are generated can be adjusted to reflect the variance of the different traits.Note that the radii can vary independently for different species (see Appendix E for illustrations).
Multiple populations belonging to the same species but differing in their trait distributions (e.g. for analyses at very large spatial scales) can also be incorporated into analyses by considering them as different operational taxonomic units specific to their residing sites.
Finally, future works could explore how it may be possible to generate random points according to different theoretical or empirical distributions, as trait variation around the mean may not be even and isotropic.This additional functionality could be implemented in the future in the open-sourced code accompanying this study for cases when such high-resolution trait data is available.

| CON CLUS IONS
The KIM method presented here computes trait β diversity by directly integrating KDEs.Out of the three different indices this method can generate, KIM V3 implements the same bandwidth for the paired communities without resampling random points to uniform distribution.KIM V3 generates values that better reflect the distribution of species in the trait space (i.e. the community trait profiles) than methods based on convex hulls or dendrograms, and also better than other methods based on KDEs.The approach is flexible, information rich and readily adapted to account for relative abundance between species and intraspecific trait variation, by using different numbers of random points and radii to generate the KDEs.
When used for the appropriate purpose, the convex hull method to inform on the trait envelope, the tree-based approaches for quantifying phylogenetic diversity, and the version of KIM using the same bandwidth and non-uniform point distribution (KIM V3), together provide a complementary set of metrics for understanding patterns of trait diversity and turnover.
Components of the tree-based approach for the computation of trait turnover between two sites whose species are part of a larger species pool.(a) Dendrogram for all species in the species pool.(b) Dendrogram for species present in site A, after trimming the global tree from a).(c) Dendrogram for species present in site B. (d) Dendrogram for species present in either site A or site B. Since species 3 does not occur at any of the two sites, branch d is not included.(e) Dendrogram for species present in both sites A and B. Although species 4 and 5 are present in only one of the two sites, branch h appears in both dendrograms, and is therefore conserved.
polytope (MCP) encompassing all species points in the trait space is computed for each community.The hypervolumes are used in the computation of the β diversity indices The distance between the original species points in the trait space is representative of how functionally different the species are Only the species with extreme trait values are used in the computation of the β diversity indices.Sensitive to outliers Cornwell et al. (2006) and Villéger et al. (2008, 2013) Tree-based method (TREE) A dendrogram is computed for the species pool based on the distance between all species based on their traits.The resulting tree is trimmed for each community, and the lengths of branches are used in the computation of the β diversity indices All species contribute to the computation of the β diversity indices.Not as sensitive to outliers as the convex hull method The distance between the species computed from the dendrograms is somewhat different from the original trait difference between species for the traits considered in the analyses Cardoso et al. (2014) Kernel density hypervolume (KDH) V1 A polytope is computed for each community based on a kernel density estimator (KDE).The hypervolumes are used in the computation of the β diversity indices The distance between the original species points in the trait space is representative of how functionally different the species are.The polytope accounts for how species are distributed in the trait space Sensitive to the threshold used to compute the MCP.The random points are resampled to uniform distribution, which diminishes the effect of the species point distributions in the trait space.A different bandwidth is used for each community to generate the random points and compute the KDE.As a result, the KDE is influenced by the distance between the species in the trait space Mammola and Cardoso (2020) Kernel density hypervolume (KDH) V2 As KDH V1, but using the same bandwidth for each pair of communities.As KDH V1, plus independence from how different species are from each other in each community As KDH V1, except for the issue related to using different bandwidths This article Kernel Integral Method (KIM) V1 The β diversity indices are directly computed from KDEs of species points in the trait space.The KDEs are based on the same procedure as KDH V1 As KDH V1, but with the advantage of being insensitive to the threshold used to compute the MCP This article Kernel Integral Method (KIM) V2 As KIM V1, but using the same bandwidth for each pair of communities As KDH V2, but with the advantage of being insensitive to the threshold used to compute the MCP This article Kernel Integral Method (KIM) V3 The β diversity indices are directly computed from kernel density estimators of species points in the trait space, but the KDEs are estimated directly from points randomly drawn around species points in the trait space, without resampling to uniform density As KDH V1 & V2, but reflecting more accurately the distribution of species points in the trait space.The radius of the hyperspheres still needs to be determined This article arranged as a square.These four species were distributed in the trait space so that the MCPs were either of different sizes (square side of lengths 4 and 2; Figure 4) or of the same size (square side of length 4; Figure A1).They were also either nested within each other, partially overlapping, or disjunct.For each of these configurations, we generated the remaining species within the community by randomly drawing species points within the MCPs according to three patterns: (i) the species points were located in a small area (square sides of length 1) in opposite corners of the MCPs (hereafter called the 'different' point distribution); (ii) the species points were randomly drawn within the MCPs (hereafter called the 'random' point distribution); (iii) the species points were located in a small area (squares of length 1) in the closest corners of the MCPs (hereafter called the 'similar' point distribution).We tested these 18 configurations for 10, 40, 70 and 100 species points, and performed analyses 50 times for each of the resulting 72 configurations (2 MCP size setups × 3 relative positions ×3 random point distributions × 4 sets of point numbers × 50 repeats = 72 configurations × 50 repeats).
Figures 5-7, FiguresA2 and A3, and their justification provided in TableA1.In summary, Jaccard dissimilarity should increase as most species points in the two communities move away from each other.
how an appropriate index of trait turnover should vary for different configurations Figures A1 and A2).This is not possible for real data with multiple pairwise comparisons between communities (thus the need for developing a new index).It is nonetheless important to explore how indices with different performances for F I G U R E 7 Changes in the contribution of replacement to overall turnover, computed as Williams replacement divided by Jaccard dissimilarity, for the seven methods summarised in Table1, for MCPs of different sizes (see Appendix C for the full set of results).The first row shows the qualitative differences in β values for trait turnover predicted under different simulated configurations of species points in the trait space when the MCPs of the two communities have different sizes.The position of the symbols on the y-axes are coarse approximations of the value an index should take for these configurations.Overall, β values, are expected to decrease as the overlap between the trait profiles of the two simulated communities increases.For detailed explanations about the changes in β values, see TableA1in Supporting Information.theoreticalconfigurations may lead to different conclusions when analysing real data.To do so, the trait diversity of plant species introduced to the Pacific islands of French Polynesia were examined, comparing trait turnover across islands and archipelagos using each method.Data were extracted from PacIFlora(Wohlwend et al., 2021a).For French Polynesia, PacIFlora contains data on naturalised non-native plant species across the 80 Pacific islands over five archipelagos: The Society Islands, the Gambier Islands, the Tuamotu Islands, the Tubuai Islands, and the Marquesas.Only the 417 naturalised species in PaciFlora appearing in the Appendix ofFourdrigniez and Meyer (2008) were used (excluding cultivated and endemic species).
Figures C1-C12).The KIM V3 method also tended to be less sensitive to the number of species points than the other KDH and KIM methods, with values and behaviours being similar from 10 to 100 species points.When communities had MCPs of the same size, the KIM V1 and V2 methods generated similar results to the KDH V1 and V2 methods, respectively, for all β diversity indices (Figures C4-C6 and C10-C12).However, when the MCPs had different sizes, the KIM methods tended to generate values more similar to each other than to the KDH methods (Figures 5-7, Figures C1-C3 and C7-C8).
2 for Williams replacement and 0.8 for the contribution of replacement to turnover (Figure 8, Figures D2 andD3).The KIM V3 and the KDH V2 methods generated the lowest Jaccard dissimilarity, and KIM V1 and TREE the highest.In contrast, KIM V3 consistently generated much higher values for the contribution of replacement to turnover than other methods, as expected because it better accounts for differences in species point distributions in the trait space.Results were similar for all combinations of traits and numbers of species used in the analyses(Figures D2 and the computation of trait turnover.Because CONVHULL only uses a subset of the species in the trait space, changing the distributions of species points within the MCP had no effect on the values of the Jaccard and Williams replacement indices.Results for the tree-based method also differed greatly from theoretical expectations compared to the KDH and KIM methods, but were closer than CONVHULL for the Jaccard index for the nested and overlapping configurations, and for Williams replacement index for the nested configuration.Overall, the tree-based method tended to either underestimate or overestimate dissimilarity in the trait profiles of communities in the MCPs, and tended to generate high values for F I G U R E 8 Application of seven trait turnover methods ( for using seed mass, plant height and SLA on 124 species, and seed mass and plant height on 124 species).Note that the order of the archipelagos is arbitrary, and lines between symbols are used as a visual aid and not to depict continuous change.theWilliams replacement index.Although the CONVHULL and treebased methods have been contrasted in the literature and have been shown to generate different results (e.g.Loiseau et al., 2017), neither of these two methods accurately reflects the trait profile of a community.In contrast, the other five trait space methods that were compared (KDH V1-V2 and KIM V1-V3) generated results more in line with theoretical expectations.These methods therefore offer a more consistent representation of the community trait profile in trait space, that is, they better capture the contribution of all species to the assessment of turnover.The computational aspects of these approaches to estimating trait turnover have nonetheless important effects on the generated β values, with potential consequences for inferences made about trait turnover in an assemblage or community.Specifically, three computational aspects of these methods were examined: (i) the use of polytopes versus kernel integrals (KDH V1 vs. KIM V1; Equations 1 and 2 vs. Equations 5 and 6); (ii) the use of the same or different bandwidths for each community in a pair when computing the KDE (KDH V1 vs. V2 and KIM V1 vs. V2); and (iii) the use of point resampling when computing the KDE (KIM V2 vs. V3).All three aspects proved to have important effects on the β diversity values calculated.Using kernel integrals, the same bandwidth and not resampling (i.e. using KIM V3) generated results most in line with theoretical expectations.
using different indices can lead to different conclusions.The KIM V3 method consistently generated much higher values for the contribution of replacement to turnover than other methods, even when randomising site-by-species matrices (Figure 8, Figures D2-D4), suggesting that the higher Jaccard dissimilarity values generated by the other methods may reflect an overestimation of the contribution of trait richness difference (the complement of replacement) to turnover.Importantly, depending on the method used, one could either conclude that most islands of an archipelago are very different from each other in terms of trait diversity (e.g.Jaccard dissimilarity values >0.5 for KIM V1 for the Gambier, Tuamotu and Society archipelagos), or very similar (Jaccard dissimilarity values mostly <0.2

Table 1 )
to empirical data on French Polynesian plant species.Jaccard dissimilarity index, Williams replacement index, and contribution of replacement to overall turnover, computed as Williams replacement divided by Jaccard dissimilarity.Results are for French Polynesia (FP) and its archipelagos, for all species, woody species and herbaceous species, using seed mass and plant height, for 250 out of 417 species (see FiguresD2 and D3