Trajectory analysis in community ecology

Trajectory analysis in community ecology. Ecological


INTRODUCTION
Ecologists have long been interested in how community dynamics operate in different kinds of ecosystems.The speed and direction of community changes, as well as the biotic and abiotic factors and processes that drive them, are central to research in community ecology (Pickett et al. 1987, McEwan et al. 2011, Vellend 2016).For example, an overarching question in ecology is how do communities respond to disturbances, depending on disturbance intensity, historical legacies, and environmental factors (Pickett and White 1985, Allison 2004, Reyer et al. 2015).A related question revolves around the existence of a single stable state under given environmental conditions, as opposed to the existence of multiple equilibria driven by priority effects (Chase 2003, Smith 2012).More generally, community ecologists are often interested in discerning the relative contribution of drift, selection, speciation, and dispersal processes to observed community dynamics and the resulting diversity patterns (Vellend 2016).Addressing such questions and goals, or others, requires both longterm experimental/observational data sets, and insightful ways of representing, comparing, and analyzing the dynamics of ecological communities.
Several statistical frameworks have been proposed and used to test hypotheses of community dynamics: (1) questions about temporal divergence/convergence of a set of communities can be addressed by testing for between-survey differences in multivariate dispersion (Anderson 2006) among other approaches (Fukami et al. 2005, Walters andCoen 2006); (2) temporal shifts in average composition in a set of monitored communities can be tested using general procedures such as distance-based or permutational MANOVA (Legendre and Anderson 1999, Anderson 2001, 2017), or specific tests based on null models (Schaefer et al. 2005); (3) spacetime interactions and complex spatiotemporal structures can be tested in different ways (Angeler et al. 2009, Legendre et al. 2010, Legendre and Gauthier 2014); (4) multivariate autoregressive models and joint dynamic species distribution models derived from latent variable modeling can be used to explicitly model species interactions and community dynamics (Hampton et al. 2013, Thorson et al. 2016, Ovaskainen et al. 2017).
Complementing statistical frameworks for hypothesis testing, the dynamics of ecological communities have traditionally been represented using scatter plots (e.g., ordination diagrams), where the plotted points correspond to the community states observed at sites at different sampling times, and temporal community trajectories are indicated using lines, arrows, or numbers (Austin 1977, Fukami et al. 2005, Trexler et al. 2005, Magalhães et al. 2007, Matthews et al. 2013).In accordance with their visualization, the formal (i.e., quantitative) study of community dynamics can be based on the geometric analysis of trajectories.For example, the site-to-site geometric comparison of trajectories can be used to group sites that exhibit the same dynamic patterns and to assess the amount of variation in community dynamics among a set of sites.Although they did not present it this way, the "second-stage community analysis" of Clarke et al. (2006) is an example of this geometry-based perspective.Clarke et al. (2006) compared pairs of trajectories by calculating the Mantel correlation between the dissimilarity matrices representing them.The pairwise comparison of all trajectories provides a new space that summarizes their similarity and can be used to explore the factors underlying differences in community dynamics (Clarke et al. 2006).Other approaches can also be included in this geometric perspective because they analyze dissimilarity values between states of a community trajectory (in other words, they analyze the shape of the trajectory) with the aim to distinguish patterns of community dynamics.Examples include "time-lag analysis," a method to detect patterns of temporal variation by comparing community dissimilarity and the time difference between surveys (Collins et al. 2000, Thibault et al. 2004, Takeuchi et al. 2010), as well as methods that analyze the time series of dissimilarity values with respect to a baseline condition (Philippi et al. 1998, Bagchi et al. 2017).Although these approaches are useful, to our knowledge, the full potential of geometrically based methods remains to be explored, especially by expanding the ways community trajectories are described and compared.Such advances would enable the development of comprehensive frameworks to describe, summarize, and interpret the spatial variation of community dynamics.
The analysis of community trajectories has an analogue in the movement of objects through space.In physics and geography, a spatial trajectory is a set of positional information for a moving object, ordered in time (e.g., trajectories of animals, humans, or hurricanes).Ecologists have been traditionally interested in the movement of animals and plant propagules in space (Turchin 1998); and a wealth of statistical methods and null models are available for this purpose (Hooten et al. 2017).Moreover, during the last decade, engineers have also further developed procedures for the analysis of movement in geographic space (Lee et al. 2007, Zheng 2015, Besse et al. 2016), and the discipline that deals with large volumes of spatial trajectory data has been called "trajectory data mining" (Buchin et al. 2013, Zheng 2015).Adapting tools designed for the study of movement to the study of community dynamics can offer novel ways of visualizing, analyzing, comparing, and ultimately understanding them.Indeed, movement metrics developed in animal ecology have been already used to categorize community dynamics into fundamental trajectories related to ecological resilience (Bagchi et al. 2017).This transfer of methodological approaches, however, needs to be done with caution, because of the difference between the geographic space in which spatial trajectories occur and the multidimensional, often non-Euclidean, spaces typical of community ecology.
With the aim of expanding the existing statistical toolbox for spatiotemporal community analysis, here, we present a community trajectory analysis (CTA) framework based on the geometric analysis and comparison of community trajectories.Our CTA framework includes (1) definition of a community trajectory in a multivariate space as formal representation of the dynamics in a community, (2) geometric properties of a community trajectory, (3) projection of a community state onto a community trajectory, (4) convergence/divergence between a pair of community trajectories, (5) geometric resemblance between a pair of community trajectories, and (6) spatial variation in community dynamics.
We study the behavior of our CTA framework using simulated community dynamics, which allows us to show how the results of the analysis may be influenced by sampling decisions.We then illustrate CTA in practice using two example applications on real spatiotemporal community data sets differing in what component of dynamics is of interest: compositional dynamics in a species-rich tropical forest and a disturbance-related structural dynamics in a species-poor forest.Finally, we discuss the advantages and limitations of CTA in relation to previous work, its domain of application, and ways in which it can be improved or extended.

Definition of a community trajectory
Let o 1 , o 2 , . .., o n be an ordered set of n community observations (n > 1) and t 1 , t 2 , . .., t n be a corresponding set of ordered times (i.e., t 1 < t 2 < ÁÁÁ< t n ).For all i in {1, 2, . .., n}, let x i contain the coordinates, or community state, corresponding to o i in a chosen multidimensional space Ω (e.g., of species composition).We define community trajectory as the sequence T = {(x 1 , t 1 ), . .., (x n , t n )}, where n is the size of the trajectory.Alternatively, the geometry of T can also be formalized using a set of n À 1 directed segments {s 1 , . .., s nÀ1 }, where s i = {x i , x i+1 } is a segment with endpoints (community states) x i and x i+1 (Fig. 1).A trajectory will often correspond to repeated observations of the same sampling unit (e.g., surveys of a permanent plot), but not necessarily.For example, a trajectory may also represent a forest chronosequence, where o 1 , o 2 , . .., o n are observations made at different sites and t 1 , t 2 , . .., t n correspond to forest ages (Foster and Tilman 2000).A trajectory may also represent an average dynamic pathway (e.g., each x i corresponding to the average composition of a set of observations).
On the multivariate space Ω.-Unlike spatial trajectories that represent the movement of objects in two-dimensional or three-dimensional Euclidean space, community trajectories represent dynamics in spaces that usually have many dimensions and may not be Euclidean.In the past, the geometric study of community trajectories has been sometimes conducted on Ω being defined by two or three axes of an ordination (Matthews et al. 2013), but this has the drawback of discarding the information of trajectories that is contained in additional dimensions.
In our framework, we assume that Ω is defined by the resemblance between pairs of community observations, measured using a dissimilarity coefficient d.All analyses that follow are based on the dissimilarity values contained in a distance matrix D = [d].This approach is not limited in the number of dimensions taken into account.
The study of community dynamics can focus on temporal changes in a wide range of properties (species composition, number and size of individuals, trait distribution, etc.), and the choice of d (hence of Ω) will depend on the property of interest.For example, here we use the percentage difference (Bray-Curtis) coefficient (Bray and Curtis 1957) to define the space Ω in which we study compositional dynamics.Some operations of our framework require d to be a metric (i.e., a distance satisfying the triangle inequality).Not all the dissimilarity coefficients used in ecology are metric, but transformations exist that help achieve this property (Gower andLegendre 1986, Legendre andLegendre 2012).Commonly used transformations like taking the square root, however, may strongly modify the angles between consecutive segments, and the overall trajectory shape.In Principal Coordinates Analysis (PCoA; Gower 1966), the usual practice of corrections for negative eigenvalues (Lingoes 1971, Cailliez 1983) can also result in strong shape distortions.In these cases, metric multidimensional scaling (mMDS; Borg and Groenen 2005) will be a better alternative.Note, however, that any Euclidean representation of semi-metric dissimilarities will introduce larger distortions than are needed to ensure metricity only.Moreover, ordination methods modify dissimilarities in different ways depending on the subset of trajectories included in matrix D. While ordination is the correct tool to display trajectories, for the calculation of CTA statistics, we advocate local transformation of semi-metric dissimilarities in every triplet when the triangle inequality is required.This creates the inconsistency that the same distance value may be transformed differently depending on the triangle that is evaluated, but in our opinion has the advantage of introducing a much smaller distortion of overall trajectory shape compared to other approaches.

Geometric properties of a trajectory
These are basic properties that describe the shape of a given trajectory, and their values may be compared between trajectories if they have been defined on comparable Ω spaces (i.e., using the same dissimilarity coefficient and similar spatial and temporal resolution of sampling; see Simulation Study).
Length and speed.-Comparing the lengths of community segments and the speeds of trajectories may be interesting, for instance, in studies focusing on the resistance and resilience of communities (Halpern 1989, Allison 2004, Rydgren et al. 2004), or to distinguish between gradual and abrupt changes (Matthews et al. 2013).The FIG. 1. Example of a community trajectory T issued from making five community observations (o 1 , . .., o 5 ) at five ordered points in time (t 1 , . .., t 5 ).Community observations are represented using a corresponding set of community states (x 1 , . .., x 5 ) in a multidimensional space (here two principal axes are shown only).The trajectory can also be formalized in terms of four directed segments (s 1 , . .., s 4 ) in the same space.The x i are positions (coordinates) and the s i are vectors that have lengths (i.e., displacement), so that the sum of segment lengths is the total length of the trajectory pathway.The same trajectory can be represented in a three-dimensional plot including the time axis (continuous arrows).
length of a segment s i is given by the distance between its two endpoints in space Ω, that is L(s i ) = d(x i , x i+1 ).The total path length of a trajectory L(T) is simply the sum of the lengths of all segments that compose it: (1) Knowing the lengths of segments s 1 , . .., s nÀ1 and time points t 1 , . .., t n , one can calculate the speed of community change along segments, S(s i ), or the average speed along the trajectory, S(T): Both lengths and speeds along trajectories have units that depend on the choice of d.If the Euclidean distance is calculated on individual counts, then lengths can be interpreted as changes in the number of individuals, but other distances may not have interpretable units.
Direction.-Anotherimportant geometric feature of a trajectory T concerns its changes in direction.Unlike in the analysis of movement in space, where angles are often measured on the X-Y plane (Turchin 1998), community trajectories are defined in spaces of more than two dimensions and this complicates the study of direction.Let {x i , x j , x k } be a triplet of community states of T that are ordered in time, that is where t i < t j < t k .If the set of distances {d(x i , x j ), d(x i , x k ), d(x j , x k )} fulfills the triangle inequality, then angles can be measured on the (Euclidean) plane that contains the three states.We define the angle 0 ≤ h(x i , x j , x k ) ≤ 180°as the change in direction of vector x j x k with respect to vector x i x j in this plane.Values of 0°indicate that the three states are completely aligned, and there is no change in direction; whereas values of 180°indicate that the orientation of the two vectors is the same but they are opposite in sense (i.e., origins and endpoints are at opposite ends of the vectors).Unlike in spatial movement, it is not useful to distinguish between positive and negative angles (e.g., between clockwise and counterclockwise turns) here, because each triplet forms its own plane.Of particular interest are the angles between consecutive segments s i and s i+1 , that is h(s i , s i+1 ) = h(x i , x i+1 , x i+2 ), because knowing the n À 1 sequence of h(s i , s i+1 ) values of a trajectory allows abrupt changes in trajectory direction to be detected (Hughes 1990, Matthews et al. 2013).For example, a large angle (i.e., a sudden change in direction) will often occur between the segment including a disturbance event and the preceding one, due to elevated mortality rates and the recruitment of new individuals/species.
In addition to inspecting the angle between consecutive segments, it is interesting to evaluate the overall directionality of the trajectory (i.e., the degree to which the community is consistently following a particular direction in Ω), which may be useful, for example, to distinguish between communities subjected to stabilizing (non-directional) selection from those influenced by directional or disruptive selection (Matthews et al. [2013]; selection sensu Vellend [2016]).We advocate directionality statistics bounded between 0 and 1, where the maximum value corresponds to a straight trajectory.For two-dimensional spatial trajectories and angles measured from the same reference direction (e.g., north), the sample mean vector length of circular statistics (Pewsey et al. 2013) is a good candidate to measure overall directionality because it is bounded between 0 and 1 and its highest value is attained when all trajectory segments are in the same direction.Unfortunately, angles are defined differently in the case of community trajectories, so a different statistic must be found.In a straight trajectory h(x i , x j , x k ) = 0°for all r time-ordered triplets; and for non-straight trajectories, the maximum penalty should occur for h = 180°.Moreover, angles of triplets corresponding to distant surveys should have a larger influence in the determination of overall directionality than angles between consecutive surveys.Therefore, we suggest overall directionality of a trajectory T is measured using the following: where h ijk = h(x i , x j , x k ), w ikj = d(x i , x j ) + d(x j , x k ) and summation are over all r time-ordered triplets.0 ≤ DIR(T) ≤ 1 and one should expect DIR(T) to be close to 1 if compositional dynamics are directional as in primary or secondary succession.Philippi et al. (1998) presented several ways to measure "progressive change."One of them is to conduct a Mantel test of the correlation between state dissimilarity and difference in survey times, using all pairs of states in the trajectory (i.e., a type of multivariate seriation test).This approach, however, is sensitive not only to changes in trajectory direction but also to changes in trajectory speed.

Projection of a community state onto a trajectory
Let y be the coordinates (state) in space Ω of o A , the community that has been sampled at a given site A. Let now T = {(x 1 , t 1 ), . .., (x n , t n )} be a trajectory in Ω representing a known dynamic pathway such as a successional sequence.Assessing the distance between y and its projection onto T is ecologically relevant because, if this distance is small enough, the future dynamics of o A could follow T. Additionally, knowing the relative position of y along T would allow the assessment of how far a given disturbed community is along a given pathway known to occur under a particular model of community dynamics.We begin by considering the comparison of y with a single directed segment s i = {x i , x i+1 }.D ps , the distance between a community state y and a segment s i , is defined as (Besse et al. 2016;Fig.2a, b) where ps stands for "point-to-segment" and proj(y, s i ) is the orthogonal projection of y onto s i .If the orthogonal projection of y onto s i falls outside the segment (Fig. 2b), the second line in Eq. 4 defines the distance between y and s i as the distance between y and the closest of the two endpoints.To calculate d(y, proj(y, s i )), the triplet {d(y, x i ), d(y, x i+1 ), d(x i , x i+1 )} again needs to fulfill the triangle inequality.The same triplet allows R (y, s i ), the relative position of proj(y, s i ) within s i , to be calculated (R 2 [0,1]): The distance from a community state y to a trajectory T, D pt , is naturally defined as the minimum distance between y and the various segments of T (Besse et al. 2016;Fig. 2c) D pt ðy; TÞ ¼ minfD ps ðy; s 1 Þ; . ..; D ps ðy; s nÀ1 Þg (6) where pt stands for "point-to-trajectory."And the relative position of proj(y, s i ) within T is where R(y, T) is a normalized measure of the position of y when projected onto T.

Convergence/divergence between a pair of trajectories
The temporal divergence/convergence of a set of communities synchronously sampled can be addressed by testing for between-survey differences in multivariate dispersion (Anderson 2006).We are interested here in testing the convergence or divergence between a particular pair of trajectories T 1 ¼ fðx  3a).

FIG. 2. (a, b)
Illustration of the calculation of the distance between a community state y and a segment s i (dashed line; Eq. 4) and the relative position of y on s i (Eq.5) when the projection of y onto s i (a) belongs to the segment and (b) when it does not.(c) Illustration of the calculation of the distance between a community state y and a trajectory T = {s 1 , s 2 , s 3 , s 4 } (Eq.6) and the relative position of y on T (Eq.7), with dotted lines representing distances between y and particular segments and the dashed line representing the minimum among them.Red lines indicate the absolute position, which leads to the relative position when compared to the total length of the segment/trajectory.T 1 Þ; . ..; D pt ðx 2m ; T 1 Þg (Fig. 3b, c).
The comparison of the two trajectories in approach 1 is symmetric, that is the two trajectories converge/diverge or not simultaneously.However, in approach 2, the analysis is asymmetric, so that T 1 may approach T 2 while the reverse may not be true.It may even happen that T 1 approaches T 2 while T 2 separates from T 1 , for example if they are one behind the other or at right angles to one another.Approach 1 addresses the question of whether the states of the two trajectories become closer/farther over time, while approach 2 addresses the question of whether the states of a given observed trajectory become closer/farther to a reference trajectory over time.In both cases, a Mann-Kendall test (Mann 1945) can be used to test for a monotonically increasing (divergence or separation) or decreasing (convergence or approach) trend.

Geometric resemblance between a pair of trajectories
Assessing the resemblance between spatial trajectories can be done in many ways (Vlachos et al. 2002, Lee et al. 2007, Besse et al. 2016, Tripathi et al. 2016).Analogously, there are multiple ways of assessing the resemblance between community trajectories.For example, in second-stage community analysis (Clarke et al. 2006), the resemblance between a pair of trajectories is measured by calculating the nonparametric Spearman correlation between their corresponding dissimilarity matrices.In the following two subsections, we develop a geometrically based approach to trajectory resemblance that takes into account the shape, size, direction, and position of trajectories but is not sensitive to differences in speed.While accounting for speed would be preferable, a purely geometric approach is easier to define and may produce similar results if the time intervals between consecutive surveys are more or less constant.The last subsection deals with the issue of excluding differences in position while keeping the other components of trajectory resemblance.
Dissimilarity between two segments.-Letus first consider D HS , the Hausdorff distance between two segments s i = {x i , x i+1 } and s j = {x j , x j+1 } (see general definition of the Hausdorff distance in Appendix S1; Besse et al. 2016; see Fig. 4a) An important limitation of D HS for the purpose of comparing community dynamics is that it does not account for the fact that directed segments have a direction (from the initial endpoint to the final endpoint).Because the direction sense is important in ecology (e.g., population growth should be distinguished from mortality), D HS will underestimate the difference between the two directed segments if they have opposite directions.To account for the direction of segments, Lee et al. (2007) proposed to explicitly incorporate the angle between segments in the definition of the distance.Alternatively, we define that a segment s i has a direction opposite to s j if the projection of the final endpoint of s i (i.e., x i+1 ) on s j has a relative position (Eq.5) smaller than that of the projection of the origin of s i (i.e., x i ) on s j , that is if R(x i+1 , s j ) < R(x i , s j ) (Fig. 4b).Visually, s i has a direction opposite to s j if the angle between them is larger than 90°, but this angle is not directly measured because the segments s i and s j have different origins.As such, we define the directed segment dissimilarity, D DS , using a modification of Eq. 8  FIG. 4. (a,b) The calculation of resemblance between a pair of directed segments s i and s j .In panel a, D HS (Eq.8) and D DS (Eq.9) are equal and the distance value results from finding the maximum (black dashed line) among the four endpoint-to-segment distances (dotted lines).When the direction of segment s j is reversed (b), to calculate D DS (Eq.9), the distance between the final endpoint of one segment and the other segment is modified (Eq.10).Here, the red dashed paths indicate the modified distances D 0 ps ðx iþ1 ; s j Þ and D 0 ps ðx jþ1 ; s i Þ. where: The distance between x i+1 and s j is modified if R(x i+1 , s j ) < R(x i , s j ).The idea underlying the modification is that if s i has a direction opposite to s j the distance from the final endpoint of s i (i.e., from x i+1 ) to s j should be larger than the distance from the initial endpoint of s i (i.e., from x i ).If R(x j+1 , s i ) < R(x j , s i ), then D 0 ps ðx jþ1 ; s i Þ is defined analogously.Ecologically, these modifications are intended to make dissimilarity between population growth and mortality processes (or between constant composition and complete species replacement) larger than dissimilarity between different degrees of population growth (or between colonization events by different species while keeping existing ones).The red dashed paths in Fig. 4b illustrate the calculation of D 0 ps ðx iþ1 ; s j Þ and D 0 ps ðx jþ1 ; s i Þ in these cases, and Fig. 4c-h shows values of D HS and D DS under different situations (note the increase in D DS from 4c to g).
Dissimilarity between two trajectories.-LetT 1 and T 2 be two trajectories to be compared.If they represent the dynamics of two sites that have been surveyed synchronously (i.e., if n = m and t 11 ¼ t 21 ; t 12 ¼ t 22 ; . ..; t 1n ¼ t 2n ), a straightforward way of comparing them is to calculate the mean across surveys of community dissimilarity between the two sites (i.e., the mean of the sequence fdðx 11 ; x 21 Þ; dðx 21 ; x 22 Þ; . ..; dðx 1n ; x 2n Þg).However, we are interested here in the comparison of dynamics beyond this specific case.As for segments, the resemblance between trajectories T 1 and T 2 can be assessed using the Hausdorff distance, but this approach and others are inappropriate because segment directions are not taken into account (Appendix S1).Instead, we advocate what we call the directed segment path dissimilarity (D DSP ) where, analogously to the distance between a community state and a trajectory (Eq.6), the distance of segment s 1i to trajectory T 2 is D DSP (Eq.11) takes into account the direction of trajectories because it is based on D DS (Eq.9), which explicitly accounts for the direction and sense of segments.D DSP is not symmetric, which precludes the application of most multivariate analysis methods unless it is symmetrized, for example by averaging D DSP (T 1 , T 2 ) and D DSP (T 2 , T 1 ): Calculating this average results in the loss of information about asymmetry in trajectory resemblance.For example, if the pathway of T 1 is a subset of the pathway in T 2 , then D DSP (T 1 , T 2 ) = 0 because one can perfectly map T 1 onto T 2 , but D DSP (T 2 , T 1 ) > 0 because T 2 includes segments that are lacking in T 1 .This distinction is lost when calculating D SDSP .Although D SDSP will be preferred to D DSP in most cases, asymmetric resemblance between trajectories may be useful, for example, to measure the degree to which the pathway of a given monitored community is similar to the (lengthier) pathway expected under a particular model of community dynamics.
Dissimilarity between centered trajectories.-Sofar, we assumed that trajectory resemblance should take into account shape, size, direction, and position.However, if d measures differences in community composition, a large amount of the variation in Ω may be due to spatial compositional differences that are time invariant (i.e., species with different abundances at different sites but whose abundance is not changing through time) and that may mask temporal changes in composition (i.e., species that are decreasing or increasing).In this context, an approach excluding differences in trajectory position may be appropriate to assess resemblance in community dynamics, for example to focus on the effect of local extinction and colonization.Such an approach can be achieved by centering all trajectories prior to the calculation of trajectory distances.Starting from the matrix D = [d] of compositional dissimilarities between states, trajectory centering can be conducted in a way that parallels the derivation of residual plots in PERMANOVA, by specifying a model matrix coding for trajectories and using linear algebra (see details in Anderson [2017]).The centering procedure can be applied to both Euclidean and non-Euclidean d measures, although in the second case, some distance values may be complex numbers after centering, which requires discarding the imaginary part (and hence will imply small distortions).The resulting centered dissimilarity matrix D cent can be used to calculate distances between trajectories, as before, using D SDSP .Note that centering of trajectories can be done with respect to the initial composition, rather than trajectory centroids, if the focus of CTA is on community divergence.

Spatial variation in community dynamics
Let {T 1 , T 2 , . .., T r } be a set of r community trajectories and let n 1 , n 2 , . .., n r (n i > 1) be the size (i.e., number of observations) of each trajectory.Then, let N P be the total number of observations (i.e., N P ¼ P r i¼1 n i ) and D be a N P 9 N P symmetric matrix calculated by applying d to all pairs of observations.Matrix D defines the resemblance between all observations (community states), which may belong to the same or different trajectories.Let now D T be an r 9 r symmetric matrix of dissimilarity values between pairs of trajectories obtained using D SDSP , or another suitable measure (Clarke et al. 2006).In the following, we describe how variation in community dynamics can be assessed and summarized with well-known multivariate methods, but using D T as input, instead of D or the usual community data tables (e.g., site-by-species compositional data).Although we require D to be a metric, this does not imply that D T will have this property (in fact, D SDSP is not a metric) and this fact must be taken into account when studying the variation in D T .
Overall variation in community dynamics.-Betadiversity is generally defined as the variation in community composition among sites within a geographical area (Whittaker 1960).Several frameworks have been proposed to quantify beta diversity (Anderson et al. 2011, Legendre andDe C aceres 2013), as well as to display and summarize spatial patterns of community variation (Dray et al. 2012).More recently, temporal beta diversity or temporal turnover have been defined as the change in community composition through time (Dornelas et al. 2014, Legendre and Gauthier 2014, Shimadzu et al. 2015).Being able to measure the resemblance between pairs of community trajectories allows variation in community dynamics to be quantified.We will use the term "dynamic beta diversity" (dBD) to refer to the overall variation of community dynamics observed in a given area during a given period of time.While in the previous sections, we allowed trajectories to be compared without restrictions, assessing dBD requires that trajectories refer to a definite set of r sampling units and to a common sampling period.In what follows, we briefly repeat part of the framework of Legendre and De C aceres (2013), adapting it to variation in community dynamics (other beta diversity frameworks could be adapted).First, the overall variation in D T is found by summing the squared dissimilarities in the upper or lower half of matrix where the notation SS is chosen by analogy with the sum-of-squares obtained from compositional data.SS Total is not a comparable measure of variation because it depends on the number of trajectories, but forms the basis for defining an index of dBD Eq. 14 converts the sum-of-squares into an unbiased estimator of variance, whose values can be compared between different data sets (Legendre and De C aceres 2013).In addition to quantifying variation in community dynamics, it is possible to split SS Total into the set of r contributions of individual community trajectories.This requires the calculation of the Gower-centered matrix G (Gower 1966) where I is an identity matrix of size r, 1 is a vector of ones and A ¼ ½a hi ¼ ½À0:5D 2 T;hi .The set of diagonal elements of G, that is diag(G) = {SS 1 , SS 2 , . .., SS r }, are the squared dissimilarities of the trajectories to the multivariate centroid of the principal coordinate space of D T .From here, the vector of local contributions to dynamic beta diversity (LCdBD) can be obtained as LCdBD indices measure the degree of uniqueness of the sampling unit i in terms of community dynamics, in comparison with the dynamics exhibited by the other sampling units.Eqs.14-17 can be applied directly to D T regardless of Euclidean or metric properties.The interpretation of LCdBD values depends on the choices of d and the distance between trajectories.For example, if d measures differences in composition and one uses D SDSP without centering trajectories, high LCdBD values may occur in sites whose initial states are similar to others but undergo unique dynamics, or in sites whose initial states were already unique.
Displaying and summarizing variation in community dynamics.-Unconstrainedordination methods, such as PCoA (Gower 1966), metric MDS, or nonmetric MDS (Borg and Groenen 2005), are commonly used to display trajectories in the space of community resemblance (i.e., matrix D).As shown by Clarke et al. (2006), these can be complemented by ordination analyses of matrix D T to display the resemblance in community dynamics.Because the matrices D and D T are related, their corresponding ordination diagrams can be displayed side by side; the former displays the similarity among community states (with trajectories joining points), and the latter displays the overall resemblance among community trajectories.Furthermore, groups of similar trajectories in the space of D can be identified by applying clustering procedures (Everitt et al. 2011) on the space of D T (or to an Euclidean representation approximating its dissimilarity values if required by the clustering algorithm), and the clustering results can be represented by colors, symbols or ellipses in the two ordination diagrams derived from D and D T .Representative trajectories may be defined using centroids or medoids in the space of D T and then interpreted using the original spatiotemporal community data (see section Examples of Application).

Software availability
To facilitate conducting CTA, functions to calculate length, angles, directionality, projection, and distances between segments/trajectories (Eqs.1-13) have been included in a new version of the R package vegclust (De C aceres et al. 2010), available in CRAN and GitHub repositories.Eqs.14-17 can be calculated using the R package adespatial (Dray et al. 2017), and multivariate analyses can be conducted using several packages.A vignette explaining how to conduct CTA using vegclust has been added to the package.

SIMULATION STUDY
Community trajectories are meant to represent community dynamics, but their geometric properties are affected by the spatial and temporal resolution of community sampling.We conducted a simulation study to illustrate the behavior and usefulness of CTA when applied to different kinds of compositional dynamics, while assessing the variability of results derived from differences in the size of sampling units and the frequency of surveys.To model species-environment relationships, we considered a theoretical environmental space with two dimensions (gradients).Species responses were simulated using Gaussian curves on each dimension independently, and initial composition was taken as that corresponding to the center of the environmental space.All simulations included random mortality of individuals, whereas species identity of recruits depended on the simulated scenario.Three kinds of community dynamics were simulated: (1) stabilizing selection, (2) post-disturbance recovery, and (3) directional selection.By "stabilizing selection," we mean that species recruitment probabilities were held constant and equal to the initial community composition, so that variability in transient compositional states (i.e., ecological drift) could be observed as a result of stochastic death and recruitment, but no long-term community change occurred.Under this scenario, CTA results depended on the interaction between ecological drift and sampling decisions, hence constituting a benchmark against which to compare the outcome of directional dynamics.The "disturbance recovery" scenario was similar to the previous one, but where 80% of the individuals had been randomly removed from the initial community.Under this scenario, the initial (post-disturbance) composition of small sampling units could differ due to the stochastic mortality caused by disturbance, but recruitment probabilities were the same as before, so that composition was directed toward the pre-disturbance dynamic equilibrium.Finally, "directional selection" included a progressive change toward a different dynamic equilibrium.This was obtained by setting recruitment probabilities to those corresponding to a position in the environmental space different from the one used to define the initial composition.Each simulated data set focused on one scenario of community dynamics and included the community dynamics of 16 different simulations (i.e., 16 trajectories).For stabilizing selection and disturbance recovery scenarios, all 16 simulations were conducted with the same parameters.For the directional selection, scenario-simulated community dynamics were paired.Here, each of eight pairs of initial communities was assigned a different target position of the environmental space (i.e., a different set of recruitment probabilities), located five units away from the center of the environmental space in the direction of one of the eight cardinal/intercardinal directions.Therefore, community dynamics within the same pair were independent stochastic realizations of the same directional pattern, but different pairs exhibited dynamics in different directions.
Compositional differences were measured using the percentage difference coefficient, and we calculated for each trajectory the total path length (Eq.1), average angle between consecutive segments, and overall directionality (DIR; Eq. 3).We also calculated the dissimilarity (D SDSP ; Eq. 13) between the 16 trajectories in each data set and the resulting dynamic beta diversity (dBD Total ; Eq. 15).Additionally, the design of the "directional selection" scenario allowed us to evaluate (1) the power of the proposed convergence/divergence test for pairs of parallel trajectories and for diverging trajectories, (2) whether the pairs of communities departing in different directions of the environmental space could be recovered by applying k-means clustering on the space of trajectory resemblance.The details of the simulation study are given in Appendix S3, and computer code to replicate simulations is provided in Data S1.We summarize the main results here: 1) Regardless of the type of dynamics simulated, the total path length of compositional trajectories generally decreases when the size of sampling units increases (due to a reduced effect of drift) or when the frequency of surveys increases (Fig. 5a).2) For stabilizing selection, the average angle between consecutive segments generally increases, and DIR decreases, when the size of sampling units increases or when the frequency of surveys decreases (Fig. 5b,  c).DIR is higher for disturbance recovery and directional selection than for stabilizing selection, the difference increasing with the size of sampling units.
3) The proposed convergence/divergence test allows trajectories departing in different directions to be distinguished from parallel trajectories.Factors affecting the power of this test are trajectory size (n) and the size of sampling units.4) The proposed dissimilarity between trajectories (D SDSP ; Eq. 13) coupled with cluster analysis allows trajectories moving in similar directions to be identified and grouped together.However, in small sampling units, the effect of ecological drift lowered the degree of agreement between clustering results and simulated dynamics.5) Regardless of the type of dynamics simulated, dBD Total tends to decrease when the size of sampling units or the frequency of surveys increases (Fig. 5d).

Compositional dynamics in the Barro Colorado Island permanent forest plot
The 50-ha (1,000 9 500 m) permanent forest plot of Barro Colorado Island (BCI; Panama; 9°9ʹ N, 79°51ʹ W) is mostly an old-growth, evergreen, tropical forest.Rainfall reaches 2,500 mm/yr but includes a four-month dry season.The plot is on a rather flat terrain, but six habitats have been distinguished by Harms et al. (2001;Fig. 6a; Table 1).Whereas a small portion of BCI consists of young forests (F), the remaining habitats correspond to old-growth forest occurring under different hydrologic regimes.The high plateau (H) is the driest area of BCI, and the andesitic cap beneath it accumulates water and creates springs along the slopes (S), which therefore have shorter duration of drought.At the other moisture extreme, the seasonal swamp (W) is inundated during the wet season, which does not occur in the low plateau (L).Finally, streamside (R) habitats occurring adjacent to seasonal streams tend to contain water into the dry season.During 1983 "El Niño," the dry season was unusually severe, lasting five months and leading to elevated mortality rates, particularly of drought-sensitive canopy trees (Condit et al. 1996).Indeed, mortality rates have been estimated to be around 3% per year for the 1982-1985 period but only 2% per year afterward (Condit et al. 1995), corresponding to 10% in a 5-yr period.Despite being influenced by spatially, and temporally, variable selection processes, community dynamics in BCI are strongly driven by ecological drift and immigration (Condit et al. 2012a).Eight plot surveys (1982, 1985, 1990, 1995, 2000, 2005, 2010, and 2015) have been conducted in BCI, spanning 33 years of compositional dynamics.
We applied our CTA framework on BCI tree data to (1) describe the overall tree compositional dynamics of the forest, (2) determine whether the different BCI habitats exhibit similar direction in compositional dynamics, (3) test for compositional convergence/divergence among habitats, and (4) test for differences in trajectory speed and directionality among habitats.
To study compositional dynamics at the whole forest level, we first calculated the basal area (m 2 /ha) of all living stems larger than or equal to 1 cm in diameter by species and survey, obtaining a data table with 8 rows (surveys) and 328 columns (tree species).Then, we built a compositional dissimilarity matrix D BCI using the percentage difference (Bray-Curtis) coefficient.Fig. 6b displays the trajectory of the BCI forest in this compositional space.The two longest segments correspond to changes that occurred between the first and second surveys (1982)(1983)(1984)(1985) and between the second and the third (1985)(1986)(1987)(1988)(1989)(1990), reflecting the above-normal mortality rates associated with the 1983 drought and subsequent recruitment.Angles between consecutive segments were all between 60°and 70°.These values can be compared with our simulations for 960 individuals and surveys every two steps (i.e., largest sample unit size and 10% individual turnover between surveys), which yielded means of 89°and 59°under stabilizing selection and directional selection, respectively.Regarding overall directionality, DIR(T) = 0.63 for BCI whereas our corresponding simulations yielded means of DIR(T) = 0.41 and DIR(T) = 0.79 under stabilizing selection and directional selection, respectively.While caution should always be taken when comparing simulations with real data sets, our comparisons are consistent with directional selection operating on the long-term compositional dynamics of the BCI forest, in addition to ecological drift.We interpret this temporal pattern as the result of species replacement derived from mortality/ recruitment during the 1983 drought and the subsequent growth of the newly established individuals.
To describe and compare compositional dynamics among BCI habitats, we first calculated the basal area  by species for 20 9 20 m quadrats, obtaining a data table with 1,250 rows (quadrats) and 328 columns (tree species).We then averaged basal area values within each habitat, obtaining a table of 48 rows (eight surveys times six habitats).The percentage difference coefficient calculated on pairs of rows of this table led to D HAB , the matrix of compositional dissimilarities between habitats.
Like the overall BCI trajectory, habitat trajectories were fastest during the first surveys (Appendix S2: Fig. S1).
Afterward, the swamp (W), the young forest (F), and stream sides (R) appeared to have faster compositional dynamics than the other habitats, but this was an artifact of their smaller area.Angles between consecutive surveys were more variable than those of the overall BCI trajectory (Appendix S2: Fig. S1).An ordination plot of the habitat trajectories in D HAB was not very informative because of large compositional distances between some of them (Fig. 7a).When we assessed the resemblance between habitat trajectories using D SDSP (Eq.13), the resulting dissimilarity matrix D T very strongly correlated with the submatrix of D HAB corresponding to the first BCI survey (r = 0.975; P = 0.0014 in a Mantel test).All these results indicated that the compositional variability among habitats obscured compositional dynamics and that D T represented mostly time-invariant differences in composition between habitats (compare Fig. 7a and b).We thus centered trajectories to focus on compositional changes (i.e., differences in trajectory shape, size, and direction), obtaining matrix D Cent HAB .Compared to the (non-centered) D HAB , matrix D Cent HAB allowed the shape and direction of habitat trajectories to be more easily displayed (Fig. 7a, c).Using D Cent HAB as input for D SDSP , we recalculated the resemblance between habitat trajectories (Fig. 7d).The resulting matrix D cent T still correlated with the submatrix of D HAB corresponding to the first survey to some extent (r = 0.795; P = 0.026 in a Mantel test), indicating that habitats similar in composition also tend to follow similar dynamics.The representation of D cent T (Fig. 7d) indicates that composition in the two plateaus (L and H) and slope (S) habitats is changing in a similar way (i.e., that the same drought-deciduous species are becoming more frequent in these habitats), whereas the swamp (W), stream sides (R), and young forest (F) are following more idiosyncratic compositional pathways (W appears to follow the same direction as L, H, and S in Fig. 7c, but the ordination plot only displays 21% of variance).The proposed asymmetric test of convergence/ divergence between pairs of habitat trajectories indicated that BCI habitats had different trends (Table 2): F is compositionally approaching the old-growth BCI forest; H is becoming compositionally more similar to S; the composition in L, while moving in a similar direction to H, is at the same time slowly diverging from it; S does not show significant trends, but seems to be approaching the composition that the plateaus had in the first surveys; finally, the composition in W is approaching that of L and converging with R. Overall, convergence in composition occurred more frequently than divergence.
Lastly, we addressed compositional dynamics at the level of 20 9 20 m quadrats.For that we assembled a data table with 10,000 rows (1,250 quadrats times 8 surveys) and 328 columns (species) and built the corresponding dissimilarity matrix D 20920 .We estimated beta diversity (BD) following Legendre and De C aceres (2013), for the eight submatrices of D 20920 corresponding to the forest surveys.The resulting BD values (0.269, 0.257, 0.251, 0.249, 0.249, 0.248, 0.245, 0.244) exhibited a negative trend (s = À0.929;P = 0.002; Mann-Kendall test), indicating an overall compositional convergence within BCI.We used one-way ANOVA to test for differences in lengths and directionalities of quadrat trajectories among habitats.Shorter trajectory lengths (hence slower dynamics) occurred in quadrats of the high plateau and young forest (F 5, 1178 = 6.605;P < 0.0001), which may be interpreted as drier conditions affecting tree growth rates.In contrast, trajectories were significantly more erratic (i.e., lower DIR values) in the swamp (F 5, 1178 = 6.712;P < 0.0001), which could be an artifact of lower tree density in this area (Table 1).

Structural dynamics on permanent plots in a mountain beech forest
In this second example, we illustrate how the CTA framework can be used to summarize and characterize dynamics in forest structure, that is, changes in the distribution of tree sizes.This contrasts to the previous example that focused on dynamics in composition.The study area is mountainous and centered on the Craigieburn Range (Southern Alps), South Island, New Zealand (43°10ʹ S, 171°35ʹ E; Fig. 8).Fuscospora cliffortioides (mountain beech) is the dominant tree species, usually forming monospecific stands.Previously, the forests consisted largely of mature stands, but a 1973 snowstorm damaged trees in 30% of stands (Harcombe et al. 1998).Profiting from woody debris created by snowfall damage and windstorms, an outbreak of the native pinhole beetle (Platypus spp., Coleoptera) and an associated fungal pathogen resulted in dispersed mortality of large trees and a decline in tree biomass over the following decade (Harcombe et al. 1998).This dispersed mortality meant that individual tree growth was strongly influenced by neighborhood competition as well as variation in site conditions (Coomes and Allen 2007).While forests were relatively unaffected by major disturbance events between 1983 and 1993, in 1994, an earthquake (M w 6.7 in magnitude), with an epicenter 10 km northwest of the study area, caused substantial size-indiscriminant mortality of trees as a result of landslides (Allen et al. 1999, Hurst et al. 2011).The Craigieburn data set consists of 250 permanent plots sampling 9,000 ha of forest, established systematically along 98 compass lines between 1970 and 1972, with line origins located randomly along streams (Fig. 8).On each 20 9 20 m plot all tree stems with diameter at breast height >3 cm had their diameter measured and recorded; tree data from nine different surveys are used in this study : 1970-1972, 1974, 1978, 1983, 1987, 1993, 1999, 2004, and 2009.These were chosen from the full set of surveys to obtain a more or less homogeneous spacing between consecutive surveys (4-6 yr).The data set includes records of 12 tree species but 99.6% of them are from F. cliffortioides.Considering the 250 plots and nine surveys gives 250 9 9 = 2,250 tree community states, but three plots were completely destroyed by the 1994 earthquake, reducing the number of states to 2,241.
We used the CTA framework to summarize structural dynamics in this forested area, which result from the interplay between stand development, natural disturbances and subsequent recruitment.We calculated resemblance between plot data in terms of size structure, using the cumulative abundance profile (hereafter CAP) approach of De C aceres et al. (2013).To account for differences in tree diameter, while emphasizing regeneration, we defined 19 quadratic diameter bins (in cm): {(2.25, 4], (4, 6.25], (6.25, 9] . . .(110.25, 121]}; species composition was also taken into account (albeit not important for this example, given the dominance of F. cliffortioides) and we took the number of stems to represent abundance.To reduce the influence of large numbers of individuals, CAPs were log-transformed before community resemblance calculations.Distance values between pairs of tree community states were calculated using a generalization of the Manhattan metric (De C aceres et al. 2013) resulting in a 2,241 9 2,241 matrix D.
As before, we used PCoA on matrix D to display community trajectories in the 250 forest plots (Fig. 9a).We then assessed the resemblance between the 250 plot trajectories using D SDSP , and PCoA was again employed to display the resemblance matrix D T (Fig. 9b).The dBD Total was 4.17, and LCdBD values are displayed as point sizes in Fig. 7b.We conducted k-means clustering on the PCoA space of D T to summarize trajectories of forest plots.For this, we used the R function cascadeKM in package vegan, where the SSI criterion suggested a clustering structure into 5, 13, or 19 groups.We chose the smallest value (k = 5) to facilitate interpretability.To characterize forest structural dynamics in each group we inspected the temporal changes in the median diametersize class distribution of F. cliffortioides trees (Fig. 10); we also calculated the median basal area for each group by survey, the median segment length/speed for each group by segment, and the median angle for each group by consecutive segments (see Appendix S2: Fig. S2; Rydgren et al. 2004).Groups 1, 2, and 4 were characterized by tree growth (i.e., shifts over time toward larger diameter classes) and self-thinning (i.e., progressive decrease in FIG. 8. Map of the study area for the analysis of structural dynamics on mountain beech forests.The location of transects (containing plots) in the study area is indicated using black lines.
the number of smaller, presumably shaded, individuals) as part of stand development.However, group 1 had higher basal area and exhibited a much narrower diameter distribution than group 4 (Fig. 10 and Appendix S2: Fig. S2), whereas group 2 was characterized by low numbers of trees spread across size classes.Like the former three, group 5 dynamics reflected growth and self-thinning, but the last two surveys showed the incorporation of new recruits.Finally, group 3 exhibited the fastest structural dynamics, including decreasing numbers of large individuals between 1970 and 1978, and a high level of recruitment since 1987 (Fig. 10 and Appendix S2: Fig. S2).Although we focused on structural dynamics in this example, we believe CTA is a useful framework to summarize the dynamic patterns of a set of sites by grouping them according to similarity in dynamics in any kind of community property, provided the user makes an informed choice of d.

Advantages, applications, and limitations of the proposed framework
We presented here a framework to describe and analyze community dynamics based on geometric analysis of trajectories in a chosen space of community dissimilarity.Unlike approaches that focus on testing for particular kinds of community dynamics (Bagchi et al. 2017, Kalyuzhny andShnerb 2017), or approaches that allow the response of communities to specific treatments to be examined with respect to dynamics under control conditions (Van Den Brink and Ter Braak 1999), CTA provides a more general way to describe and compare trajectories.Among other advantages, CTA allows the relationships in community dynamics to be summarized in dissimilarity values.As a result, any dissimilaritybased statistical tool designed for the (static) analysis of community variation can be also used to address questions of variation in community dynamics.Therefore, the goals of CTA are very different to those of recent statistical frameworks that pursue an explicit modelling and prediction of species interactions and community dynamics (Hampton et al. 2013, Thorson et al. 2016).
A key element of our contribution is the ability to evaluate resemblance in community dynamics between sampling units through the geometric comparison of their representation as trajectories, which leads to the definition of a space of trajectory resemblance.Clarke et al. (2006) developed a similar idea in their "secondstage" community analysis, a term coined by Somerfield and Clarke (1995).In their approach, the resemblance between a pair of trajectories T 1 and T 2 is measured by calculating the nonparametric Spearman correlation between the two dissimilarity matrices that describe each trajectory.This requires synchronous surveys of the sampling units, because one needs to match the elements of the two dissimilarity matrices (following our notation, one needs to define fdðx 1i ; x 1j Þ; dðx 2i ; x 2j Þg pairs, for all i and j surveys).As it discards the information contained in dissimilarities between states of one trajectory and states of the other (i.e., d(x 1i , x 2j ) values), the Mantel correlation is appropriate to compare the shape of the two trajectories, but fails include several other geometric aspects.Moreover, the correlation will be À1 if the pathway of T 1 goes from x to y and the pathway of T 2 follows an exact inverse pathway from y to x, but it will be +1 if both trajectories have the same exact overall shape but T 1 goes from x to y and T 2 goes from x to z.Therefore, the sense is the only aspect of trajectory direction to which this approach is sensitive.Our framework extends the approach of Clarke et al. (2006)  match.Second, it allows differences in geometric aspects beyond shape (i.e., position, size, and direction) to be taken into account.Finally, in CTA, we showed how to compare compositional trajectories after centering them (i.e., after removing differences in trajectory position), which allows static spatial variation to be discarded to focus on the spatiotemporal component of community resemblance.Clarke et al. (2006) claimed that their approach allowed this, but in our opinion, it suffers from the limitation of not appropriately accounting for differences in direction, and neglecting differences in size, aspects of geometry that may be relevant to study spatiotemporal interaction.
There are multiple applications of CTA.(1) Basic geometric properties of community trajectories (length, speed, or overall directionality) can be compared across sites or regressed against explanatory factors.For example, it is increasingly accepted that the order and timing of species arrival during community assembly (i.e., priority effects) alters ecosystem structure and functioning (Tan et al. 2012), but it is much less clear how variation in community assembly history influences the rate at which trajectories diverge or the speed of community changes.(2) Projection of community states onto trajectories, asymmetric convergence/divergence tests, and trajectory asymmetric resemblance assessments are three tools that can be used to study the dynamics observed at particular sampling units with respect to pathways defined a priori from known models of community dynamics, such as models of primary/secondary succession or seasonal cyclical dynamics.(3) Classification and ordination analyses on (symmetric) matrix D T can be used to display and summarize patterns of community dynamics in a set of monitored sites (or the output of simulation models), as illustrated here for forest structural dynamics.(4) Instead of comparing community dynamics among sampling units, dBD Total values (and local contributions) can be used to compare the overall variation in community dynamics among sets of sampling units surveyed with comparable sampling methods, analogously to the comparison of beta diversity across forests (De C aceres et al. 2012).(5) In the same way that multivariate regression models are routinely used in ecology to hypothesize drivers of observed static patterns (Legendre and Legendre 2012), the application of distance-based regression frameworks, such as PERMA-NOVA (Anderson 2001(Anderson , 2017)), on matrix D T could be used to test hypotheses about the causes of variation in community dynamics.This use of CTA was already illustrated by Clarke et al. (2006), who employed one-way ANOSIM (Clarke and Green 1988) to test for the consistency of marine benthic community dynamics across sites.The numerous papers applying their method (Caro et al. 2010, Kroeker et al. 2013, Schulz et al. 2013) demonstrate the usefulness of analyzing the variation in D T against all kinds of factors.Our contribution provides a more robust basis for defining D T and a broader framework within which to naturally embed such analyses.When analyzing the variation in D T , it is important to remember that the beginning and end of a trajectory is dictated by sampling and is often arbitrary (Dornelas et al. 2012).This means the influence of the starting point must be considered when interpreting community dynamics in terms of potential predictors.
Regardless of the application, users of CTA should also be aware of the multiple factors that may affect the geometry of community trajectories, most importantly the choice of d but also the different artifacts derived from spatiotemporal sampling decisions and their interaction with ecological processes.To be able to compare CTA results among different data sets, trajectory analyses should be conducted after homogenizing as much as possible the size of sampling units (e.g., by aggregation of data in data sets whose sampling unit size is smallest) and the frequency of surveys (e.g., by selecting surveys of data sets with the highest temporal resolution), so that the effects of temporal and spatial resolution affect all data sets in the same way.Null models of community dynamics coupled with of sampling decisions may also be helpful as benchmarks to provide context for CTA results in real data sets, as we did in our simulation study.

Improvements and extensions of the current framework
When developing our framework, we evaluated the degree to which procedures from spatial movement analysis or spatial trajectory data mining should be modified to be used in the multidimensional spaces typical of community ecology.This led us to propose new definitions of trajectory angles and overall directionality and to define indices to compare trajectories while incorporating the direction of segments D DS and D SDSP .The behavior of DIR and D SDSP was satisfactory in our study of spatiotemporal dynamics using simulations and in the two examples above, but better alternatives may be defined.Hence, additional work will be needed to decide which properties are desirable for the appropriate comparison of the dynamics of ecological communities.In the same way, we calculated D SDSP after centering of trajectories to discard differences in trajectory position; methods related to Procrustes analysis could be tested and added to our framework to compare trajectories solely in terms of their shape (Adams and Collyer 2009).
A potential extension of this framework concerns the analysis of matrix D S , containing dissimilarities between segments.Instead of clustering trajectories, we could have chosen to group directed segments, from the same or different trajectories, to define what could be called "consensus community trajectories."In trajectory data mining, spatial consensus trajectories have been obtained using density-based cluster analysis applied to D S (Lee et al. 2007, Tripathi et al. 2016).An example of consensus community trajectory in forest structural dynamics is the one representing stand development.
As in other interdisciplinary tools such as cluster analysis, literature on methods of trajectory analysis lies dispersed across different disciplines (mathematics, engineering, movement ecology, evolution, etc.).Therefore, trajectory analysis methods of which we are unaware may be adapted and added to the CTA toolbox.Searching the literature of these related fields may be a fruitful endeavor to enlarge and improve the present framework.Further, although we presented our framework as focused on community ecology, CTA can readily be used to study trajectories in other fields where temporal series are multivariate or where spatial transects result in trajectories in multivariate spaces (Adams andCollyer 2009, Lohman et al. 2017).

FIG. 3 .
FIG.3.Illustration of different ways to study the convergence/divergence between two trajectories T 1 (dashed segments) and T 2 (continuous segments) of sizes n and m, respectively: (a) symmetric convergence analysis by analyzing the sequence of n = m distances between synchronous surveys; (b) asymmetric analysis of the n distances from states of T 1 to T 2 ; (c) asymmetric analysis of the m distances from states of T 2 to T 1 .
FIG. 4. (a, b)  The calculation of resemblance between a pair of directed segments s i and s j .In panel a, D HS (Eq.8) and D DS (Eq.9) are equal and the distance value results from finding the maximum (black dashed line) among the four endpoint-to-segment distances (dotted lines).When the direction of segment s j is reversed (b), to calculate D DS (Eq.9), the distance between the final endpoint of one segment and the other segment is modified (Eq.10).Here, the red dashed paths indicate the modified distances D 0 ps ðx iþ1 ; s j Þ and D 0 ps ðx jþ1 ; s i Þ. (c-h) Examples of D HS and D DS values for different pairs of directed segments (assuming Cartesian coordinates and the Euclidean distance for d).In panels c-e, dashed lines indicate the (equal) value of D HS and D DS ; in panels f-h, black dashed lines have lengths equal to D HS values, and red dashed paths have lengths equal to D DS values.Note the progressive increase in D DS when the second segment is tilted until its direction has opposite sense with respect to the first.
FIG. 5. Results of the simulation study for (a) mean trajectory length, (b) mean angle between consecutive segments, (c) mean trajectory directionality, and (d) dynamic beta diversity of the data set.Each panel shows the mean value of these statistics across 20 simulated data sets for combinations of the size of sampling units (in our case, the maximum number of individuals in the simulated community) and the number of steps between surveys.Columns indicate the kind of simulated dynamics.
FIG. 6.(a) Barro Colorado Island (BCI; Panama) forest topography and habitat classification according to Harms et al. (2001).Habitat types are W, swamp; S, slopes; R, streamside; L, low plateau; H, high plateau; F, young forest.(b) Ordination of the community states of the overall BCI compositional trajectory (PCoA).Survey years are indicated, as well as turning angles between consecutive segments.

FIG. 7 .
FIG. 7. Barro Colorado Island forest habitat trajectories.(a) Ordination of the habitat compositional states before centering (D HAB ), with trajectories indicated using arrows; (b) ordination of dissimilarities between trajectories (D T ) derived from D HAB ; (c) ordination of habitat compositional states after centering (D Cent HAB ), with trajectories indicated using arrows; (d) ordination of dissimilarities between trajectories (D cent T ) derived from D Cent HAB .The four ordinations were obtained by PCoA with correction factors for negative eigenvalues indicated in the top right corner of each panel.
FIG. 9. (a) Ordination of tree community states (PCoA on D, after applying Cailliez's [1983] correction for negative eigenvalues; c = 88.51), with trajectories indicated using arrows; (b) ordination of plot trajectories in the space of structural dynamics (PCoA on D T , after applying correction for negative eigenvalues; c = 42.86).Arrow colors in panel a and colored symbols in panel b indicate the five groups obtained using k means.Point sizes in panel b are proportional to LCdBD (local contributions to dynamic beta diversity) values.

FIG. 10 .
FIG. 10.Temporal changes in diameter class distribution for different forest plot groups.For each diameter class and group, the median across plots belonging to that group is shown.Different surveys are indicated using gray tones.Note log scale of y-axes.

TABLE 1 .
Number of 20 9 20 quadrats, number of trees, mean tree density per quadrat, quadrat mean trajectory length (L 20920 ) and quadrat mean directionality (DIR 20920 ) by Barro Colorado Island (BCI; Panama) forest habitats and for the whole forest.
Note: Letters indicate homogeneous groups of habitats according to Tukey's honest significant differences.

TABLE 2 .
(Mann 1945)the asymmetric convergence/divergence test for combinations of habitats in Barro Colorado Island.Values in the table correspond to the statistic (tau) and significance level of the Mann-Kendall test(Mann 1945)on the sequence of distances between states of the row habitat trajectory with respect to the column habitat trajectory.Negative values indicate that the row trajectory is approaching the column trajectory, and positive values indicate that it is moving away from it.*P ≤ 0.05; **P ≤ 0.01; ***P ≤ 0.001; †P ≤ 0.1; NS, P > 0.1. Notes: