History of ETA

• In 2019, De Cáceres et al, presented Community Trajectory Analysis (CTA) as a framework for trajectory analysis in community ecology.

• The initial framework was extended with new metrics and visualisation modes in Sturbois et al. (2021a).

• The same framework was applied to stable isotope data in Sturbois et al. (2021b), referring to it as Stable Isotope Trajectory Analysis (SITA).

• Recently, the framework has been extended to cyclical data (Djeghri et al. in prep.), under the name of Cyclical Ecological Trajectory Analysis (CETA).

Flexibility

Since it can be applied to multiple target entities and multivariate spaces, we collectively refer to the framework as Ecological Trajectory Analysis (ETA).

Installation

Package ecotraj (ver. 1.1.0) is distributed from CRAN and can be installed using:

install.packages("ecotraj")

More frequent updates can be obtained if installing from GitHub (compilation tools are required):

remotes::install_github("emf-creaf/ecotraj")

Documentation

Although the package comes with function documentation, the complete function reference and articles explaining how to use the package can be found at https://emf-creaf.github.io/ecotraj/.

Loading

We begin by loading ecotraj (ver. 1.1.0) to the workspace:

library(ecotraj)

Metadata

In this session we will primarily use a small data set where three entities have been observed four times:

entities <- c("1","1","1","1","2","2","2","2","3","3","3","3")
times <- c(1.0,2.0,3.0,4.0,1.0,1.75,2.5,3.25,1.0,1.5,2.0,2.5)

Distance matrix

We will assume that ecological states can be compared using the Euclidean distance and only two variables have been measured:

xy<-matrix(0, nrow=12, ncol=2)
xy[2,2]<-1
xy[3,2]<-2
xy[4,2]<-3
xy[5:6,2] <- xy[1:2,2]
xy[7,2]<-1.5
xy[8,2]<-2.0
xy[5:6,1] <- 0.25
xy[7,1]<-0.5
xy[8,1]<-1.0
xy[9:10,1] <- xy[5:6,1]+0.25
xy[11,1] <- 1.0
xy[12,1] <-1.5
xy[9:10,2] <- xy[5:6,2]
xy[11:12,2]<-c(1.25,1.0)

The distance matrix is obtained using:

d <- dist(xy)

To build them, we need to combine the distance matrix (\(\mathbf{D}\)) and the entity/survey information in a single object using function defineTrajectories():

x <- defineTrajectories(d, entities)

The function returns an object (a list) of class trajectories that contains all the information for analysis:

class(x)

[1] "trajectories" "list"        

This object contains two elements:

names(x)

[1] "d"        "metadata"

Element d contains the input distance matrix (\(\mathbf{D}\)), whereas metadata is a data frame including information of observations:

x$metadata

   sites surveys times
1      1       1     1
2      1       2     2
3      1       3     3
4      1       4     4
5      2       1     1
6      2       2     2
7      2       3     3
8      2       4     4
9      3       1     1
10     3       2     2
11     3       3     3
12     3       4     4

In our examples we will also use another object where we include observation times:

xt <- defineTrajectories(d, entities, times = times)
xt$metadata

   sites surveys times
1      1       1  1.00
2      1       2  2.00
3      1       3  3.00
4      1       4  4.00
5      2       1  1.00
6      2       2  1.75
7      2       3  2.50
8      2       4  3.25
9      3       1  1.00
10     3       2  1.50
11     3       3  2.00
12     3       4  2.50

One particularity of this second trajectory object xt is that trajectories are not synchronous (observation times are not the same across trajectories):

is.synchronous(x)

[1] TRUE

is.synchronous(xt)

[1] FALSE

Since \(\Omega_0\) has only two dimensions in this example, the Principal Coordinates Analysis (PCoA) on matrix \(\mathbf{D}\) displays the complete space.

One can obtain the length of trajectory segments and the total path length:

trajectoryLengths(x)

  S1       S2        S3     Path
1  1 1.000000 1.0000000 3.000000
2  1 0.559017 0.7071068 2.266124
3  1 0.559017 0.5590170 2.118034

In addition to segment lengths, one can also calculate the distance (i.e. the Net change) between all states and the initial, which is taken as reference:

trajectoryLengths(x, relativeToInitial = TRUE)

  Lt1_t2   Lt1_t3   Lt1_t4     Path
1      1 2.000000 3.000000 6.000000
2      1 1.520691 2.136001 4.656692
3      1 1.346291 1.414214 3.760505

When observation times are available, it may be of interest to calculate segment or trajectory speeds using:

trajectorySpeeds(xt)

        S1       S2       S3     Path
1 1.000000 1.000000 1.000000 1.000000
2 1.333333 0.745356 0.942809 1.007166
3 2.000000 1.118034 1.118034 1.412023

Finally, one may calculate the internal variation of states within each trajectory using:

trajectoryInternalVariation(x)

       ss_1      ss_2      ss_3      ss_4 internal_ss internal_variance
1 2.2500000 0.2500000 0.2500000 2.2500000    5.000000         1.6666667
2 1.3281250 0.0781250 0.1406250 1.0156250    2.562500         0.8541667
3 0.8007812 0.1757812 0.2070312 0.4257813    1.609375         0.5364583

Zero angles indicate that the three points (e.g. the two consecutive segments) are in a straight line. The larger the angle value, the more is trajectory changing in direction.

Function trajectoryAngles() allows calculating the angles between consecutive segments:

trajectoryAngles(x)

     S1-S2    S2-S3     mean         sd       rho
1  0.00000  0.00000  0.00000 0.00000000 1.0000000
2 26.56505 18.43495 22.50000 0.07097832 0.9974842
3 63.43495 53.13010 58.28253 0.08998746 0.9959593

Mean and standard deviation statistics of angles are calculated according to circular statistics.

We can use the same function to calculate angles relative to the initial state:

trajectoryAngles(x, relativeToInitial = TRUE)

     t1-S2    t1-S3     mean        sd       rho
1  0.00000  0.00000  0.00000 0.0000000 1.0000000
2 26.56505 35.53768 31.05136 0.0783410 0.9969360
3 63.43495 94.76364 79.09930 0.2751266 0.9628599

To measure the overall directionality of a trajectory (i.e. if the path consistently follows the same direction in \(\Omega\) ), we recommend using function trajectoryDirectionality():

trajectoryDirectionality(x)

        1         2         3 
1.0000000 0.8274026 0.5620859 

Another function, called trajectoryWindowMetrics() calculates trajectory metrics on moving windows over trajectories, but will not be illustrated here.

Function trajectoryProjection() allows obtaining the relative position of each ecological state of a trajectory:

trajectoryProjection(x$d, 1:4, 1:4)

  distanceToTrajectory segment relativeSegmentPosition
1                    0       1                       0
2                    0       1                       1
3                    0       2                       1
4                    0       3                       1
  relativeTrajectoryPosition
1                  0.0000000
2                  0.3333333
3                  0.6666667
4                  1.0000000

The same function can also be used to perform an orthogonal projection of arbitrary ecological states onto a given reference trajectory.

For example we can study the projection of third state of the trajectory of entity ‘2’ (i.e. state 7) onto the trajectory of entity ‘1’ (i.e. states 1 to 4), which happens to be in the half of the trajectory:

trajectoryProjection(x$d, 7, 1:4)

  distanceToTrajectory segment relativeSegmentPosition
7                  0.5       2                     0.5
  relativeTrajectoryPosition
7                        0.5

We can quantify those differences using function trajectoryShifts(), which internally uses orthogonal projection.

To illustrate this function, we will use a small data set of two parallel trajectories, but where the second is shifted:

entities3 <- c("1","1","1","1","2","2","2","2")
times3 <- c(1,2,3,4,1,2,3,4)
xy3<-matrix(0, nrow=8, ncol=2)
xy3[2,2]<-1
xy3[3,2]<-2
xy3[4,2]<-3
xy3[5:8,1] <- 0.25
xy3[5:8,2] <- xy3[1:4,2] + 0.5 # States are all shifted with respect to entity "1"
x_shift <- defineTrajectories(dist(xy3), entities3, times = times3)

If the trajectories are synchronous (i.e., their observations where done at the same times), one can analyze the sequence of distances between consecutive observations and address the question of whether the states of the two trajectories become closer/farther over time using a statistical test:

trajectoryConvergence(x, type = "pairwise.symmetric")

$tau
          1         2         3
1        NA 0.9128709 0.9128709
2 0.9128709        NA 0.9128709
3 0.9128709 0.9128709        NA

$p.value
          1         2         3
1        NA 0.1485617 0.1485617
2 0.1485617        NA 0.1485617
3 0.1485617 0.1485617        NA

Note that this test can be straightforwardly generalized to multiple trajectories:

trajectoryConvergence(x, type = "multiple")

$tau
[1] 0.9128709

$p.value
[1] 0.1485617

Finally, trajectory convergence/divergence can also be studied in non-synchronous trajectories. In this case, the question addressed is whether the states of the target trajectory become closer/farther to the reference trajectory over time.

trajectoryConvergence(xt, type = "pairwise.asymmetric")

There are multiple ways of assessing the dissimilarity in dynamics, e.g.:

Function trajectoryDistances() allows calculating several of them, and returns a distance matrix containing the dissimilarity between pairs of trajectories:

trajectoryDistances(xt, distance.type = "TSPD")

          1         2
2 0.6276322          
3 1.2065718 0.6822515

Analogously to trajectoryInternalVariation(), function dynamicVariation() returns the sum of squares of dynamic variation, an unbiased dynamic variance estimator and the relative contribution of individual trajectories to the overall sum of squares.

Function dynamicVariation(), makes internal calls to trajectoryDistances(), which means that we may get slightly different results if we change the trajectory dissimilarity coefficient:

dynamicVariation(xt, distance.type = "TSPD")

$dynamic_ss
[1] 0.7717349

$dynamic_variance
[1] 0.3858674

$relative_contributions
         1          2          3 
0.46561876 0.03786032 0.49652092 

Trajectory interpolation allows recalculating positions along trajectory pathways so that observation times are the same across all trajectories, hence obtaining a synchronous data set.

For example, here we interpolate trajectories in xt to times c(1, 1.5, 2, 2.5) (the observation times of entity ‘3’) to obtain a synchronous data set:

xt_inter <- interpolateTrajectories(xt, times = c(1, 1.5, 2.0, 2.5))

The following trajectory plots show the effect of interpolation visually:

Trajectory centering is done using function centerTrajectories():

x_cent <- centerTrajectories(x)

The following trajectory plots show the effect of centering visually:

Temporal smoothing of trajectories done using function smoothTrajectories(), which applies a multivariate moving average over each trajectory and uses a Gaussian kernel to specify average weights.

The following plots illustrate the effect of smoothing: