Ecological Dynamic Regime framework

1.1. Ecological Dynamic Regimes and the EDR framework

An Ecological Dynamic Regime (EDR) is defined as the “fluctuations of ecosystem states around some trend or average resulting from an intricate mix of internal processes and external forces that, in the absence of perturbations, keep the system within specific domains of attraction” (Sánchez-Pinillos et al., 2023).

EDRs are composed of multiple ecological trajectories showing similar processes in the development, interaction, and reorganization of some state variables over time. Thus, ecological trajectories are the essential units of EDRs and can be defined as sequences of observations capturing the temporal changes of some state variables. Usually, ecological trajectories are defined based on multiple state variables (e.g., species, functional traits, land uses) that make challenging the characterization and comparison of EDRs.

The EDR framework is a set of algorithms and metrics aiming to characterize and compare EDRs composed of ecological trajectories in multidimensional state spaces. You can find more information in the following publication, including the formal definition of the EDR concept, the methodological details of the framework, and some illustrative examples with artificial and real data sets:

• Sánchez-Pinillos M., Kéfi, S., De Cáceres, M., Dakos, V. 2023. Ecological Dynamic Regimes: Identification, characterization, and comparison. Ecological Monographs. https://doi.org/10.1002/ecm.1589

1.2. About this vignette

This vignette aims to give you an overview of the EDR framework and its implementation in ecoregime to characterize and compare EDRs.

You can install ecoregime directly from CRAN or from my GitHub account (development version):

# install.packages("ecoregime")
# devtools::install_github(repo = "MSPinillos/ecoregime", dependencies = T, build_vignettes = T)

Once you have installed ecoregime you will have to load it:

library(ecoregime)

citation("ecoregime")
#> To cite 'ecoregime' in publications use:
#> 
#>   Sánchez-Pinillos M, Kéfi S, De Cáceres M, Dakos V (2023). "Ecological
#>   dynamic regimes: Identification, characterization, and comparison."
#>   _Ecological Monographs_, e1589. <https://doi.org/10.1002/ecm.1589>.
#> 
#>   Sánchez-Pinillos M, Dakos V, Kéfi S (2024). "Ecological dynamic
#>   regimes: A key concept for assessing ecological resilience."
#>   _Biological Conservation_, 110409.
#>   <https://doi.org/10.1016/j.biocon.2023.110409>.
#> 
#>   Sánchez-Pinillos M (2023). _ecoregime: Analysis of Ecological Dynamic
#>   Regimes_. <https://doi.org/10.5281/zenodo.7584943>.
#> 
#> To see these entries in BibTeX format, use 'print(<citation>,
#> bibtex=TRUE)', 'toBibtex(.)', or set
#> 'options(citation.bibtex.max=999)'.

3.1. Identifying EDRs from clustering analyses

Identifying EDRs involves finding groups of ecological trajectories showing dynamic patterns more similar between each other than with any other ecological trajectory in the trajectory space. In this sense, it is possible to use the trajectory dissimilarity matrix (\(D_T\)) to identify EDRs through clustering. However, applying clustering analyses to trajectory data is not always straightforward. Depending on the characteristics of your data, your goals, and the system that you are evaluating, the choice of the clustering algorithm may have important consequences.

Currently, this is out of the scope of ecoregime but you can find more details in the paper introducing the EDR framework (Sánchez-Pinillos et al., 2023)

3.2. Defining EDRs based on ecological properties

Clustering analyses are not the only way to define EDRs. As I stated before, the ecological trajectories of our simulated data can be considered part of the same EDR for having similar dynamic patterns. Thus, you could define groups of trajectories that you want to evaluate for sharing similar aspects in the state variables. For example, you could consider groups of trajectories with the same number of observations in which the state variables change following similar patterns.

4.1. Identifying representative trajectories with RETRA-EDR

RETRA-EDR (REpresentative TRAjectories in Ecological Dynamic Regimes) is an algorithm that aims to identify representative trajectories in EDRs based on the distribution and density of their ecological trajectories (Sánchez-Pinillos et al., 2023). Whereas many algorithms with similar goals where developed for moving objects in low-dimensional spaces, RETRA-EDR can be applied to high-dimensional spaces, capturing the complexity of empirical data.

I will not explain the details of RETRA-EDR in this vignette (you can find them in Sánchez-Pinillos et al., 2023), but it is important that you know its main steps to understand the output returned by retra_edr():

1. The ecological trajectories forming the EDR are split into segments. Then, a segment space analogous to the trajectory space is generated from a matrix containing segment dissimilarities.
2. The segment space is divided into regions with a minimum number of trajectory segments (minSegs). For that, RETRA-EDR recursively divides the space into halves following the axes of the segment space and generating a kd-tree structure. The greater the density of segments in a particular region, the larger the number of partitions required until obtaining regions with minSegs segments. The number of partitions is known as the depth of the kd-tree (kdTree_depth) and the number of segments in the final region is the Density.
3. For each region with at least minSegs segments, RETRA-EDR extracts the medoid as the representative segment.
4. All representative segments are joined forming a network of representative trajectories.
5. RETRA-EDR returns a set of attributes characterizing the representative trajectories identified in the EDR, including the value of minSegs; the sequence of representative segments (Segments) and the number of states (Size) forming each representative trajectory; the total length of the trajectory (Length); the links (Link) generated to join representative segments and the dissimilarity between the connected states (Distance); the number of segments in each region represented by each representative segment (Density), and the number of partitions of the segment space until obtaining a region with minSegs segments or less (kdTree_depth).

Let’s apply retra_edr() to our data. As the EDR only contains 10 trajectories, we will assign a value of minSegs relatively low.

# Use set.seed to obtain reproducible results of the segment space in RETRA-EDR
set.seed(123)

# Apply RETRA-EDR
repr_traj <- retra_edr(d = state_dissim, trajectories = rep(ID_sampling, each = 3),
                   states = rep(ID_obs, 10), minSegs = 2)

RETRA-EDR has returned two representative trajectories. We can summarize their attributes using summary():

summary(repr_traj)
#>    ID Size    Length   Avg_link   Sum_link Avg_density Max_density Avg_depth
#> T1 T1    6 0.3577606 0.02996245 0.05992491           3           3  3.666667
#> T2 T2    7 0.6048257 0.03925918 0.07851836           3           3  3.500000
#>    Max_depth
#> T1         4
#> T2         4

The output of retra_edr() is an object of class RETRA. We can extract specific attributes as we do with lists.

Segments is probably the most important attribute since it informs about the states that form the representative trajectories.

• T1 is composed of six states of three different trajectories (B, I, and E).
• T2 is composed of the three states of the sampling unit C and converges towards T1 through a link between the last state of C and the second state of I.

lapply(repr_traj, "[[", "Segments")
#> $T1
#> [1] "B[1-2]" "I[2-3]" "E[2-3]"
#> 
#> $T2
#> [1] "C[1-2]" "C[2-3]" "I[2-3]" "E[2-3]"

4.2. Visualizing representative trajectories in the state space

We can also visualize the representative trajectories in the state space using the function plot().

# Plot the representative trajectories of an EDR
plot(repr_traj, # <-- This is a RETRA object returned by retra_edr()
     # data to generate the state space
     d = state_mds, trajectories = rep(ID_sampling, each = 3), states = rep(ID_obs, 10),
     # use the colors previously used for individual trajectories. 
     # (make them more transparent to highlight representative trajectories)
     traj.colors = alpha(traj.colors, 0.3),
     # display representative trajectories in blue    
     RT.colors = "blue",
     # select T2 to be displayed with a different color (black)
     select_RT = "T2", sel.color = "black",
     # Identify artificial links using dashed lines (default) and a different color (red)
     link.lty = 2, link.color = "red",
     # We can use other arguments in plot()
     main = "Representative trajectories")

# Add a legend
legend("bottomright", c("T1", "T2", "Link"), 
       col = c("blue", "black", "red"), lty = c(1, 1, 2))

4.3. Defining representative trajectories based on trajectory features

RETRA-EDR returns the longest possible trajectory according to the procedure defined before. However, you could be interested in only a fraction of a representative trajectory, split representative trajectories based on some criteria (e.g. long artificial links), or maybe, define your own representative trajectory.

Let’s say that we want to split T1 and T2 from the last segment ("E[2-3]"). For that, we can use the function define_retra().

There are two ways of defining our “new” trajectories. The first one consists in providing all the details in a data frame:

# Generate a data frame indicating the states forming the new trajectories
new_traj_df <- data.frame(
  # name of the new trajectories (as many times as the number of states)
  RT = c(rep("T1.1", 4), rep("T2.1", 5), rep("T1.2", 2)),
  # name of the trajectories (sampling units)
  RT_traj = c(rep("B", 2), rep("I", 2), # for the first trajectory (T1.1)
              rep("C", 3), rep("I", 2), # for the second trajectory (T2.1)
              rep("E", 2)), # for the third trajectory (T1.2)
  # states in each sampling unit 
  RT_states = c(1:2, 2:3, # for the first trajectory (T1.1)
                1:3, 2:3, # for the second trajectory (T2.1)
                2:3), # for the third trajectory (T1.2)
  # representative trajectories obtained in retra_edr()
  RT_retra = c(rep("T1", 4), rep("T2", 5), 
               rep("T1", 2)) # The last segment belong to both (T1, T2), choose one
)

new_traj_df
#>      RT RT_traj RT_states RT_retra
#> 1  T1.1       B         1       T1
#> 2  T1.1       B         2       T1
#> 3  T1.1       I         2       T1
#> 4  T1.1       I         3       T1
#> 5  T2.1       C         1       T2
#> 6  T2.1       C         2       T2
#> 7  T2.1       C         3       T2
#> 8  T2.1       I         2       T2
#> 9  T2.1       I         3       T2
#> 10 T1.2       E         2       T1
#> 11 T1.2       E         3       T1

The second one consists in providing a list with the sequences of segments defining each trajectory. This is particularly useful when we want to select a subset in the output of RETRA-EDR:

# List including the sequence of segments for each new trajectory
new_traj_ls <- list(
  # First part of T1 excluding the last segment
  repr_traj$T1$Segments[1:(length(repr_traj$T1$Segments) - 1)],
  # First part of T2 excluding the last segment
  repr_traj$T2$Segments[1:(length(repr_traj$T2$Segments) - 1)],
  # Last segment of T1 and T2: segment composed of states 2 and 3 of the sampling unit E ("E[2-3]")
  "E[2-3]"
)

new_traj_ls
#> [[1]]
#> [1] "B[1-2]" "I[2-3]"
#> 
#> [[2]]
#> [1] "C[1-2]" "C[2-3]" "I[2-3]"
#> 
#> [[3]]
#> [1] "E[2-3]"

In any case, define_retra() will return the same output

# Define representative trajectories using a data frame
new_repr_traj <- define_retra(data = new_traj_df, 
             # Information of the state space
             d = state_dissim, trajectories = rep(ID_sampling, each = 3), 
             states = rep(ID_obs, 10), 
             # RETRA object returned by retra_edr()
             retra = repr_traj)

# Define representative trajectories using a list with sequences of segments
new_repr_traj_ls <- define_retra(data = new_traj_ls, 
             # Information of the state space
             d = state_dissim, trajectories = rep(ID_sampling, each = 3), 
             states = rep(ID_obs, 10), 
             # RETRA object returned by retra_edr()
             retra = repr_traj)

if (all.equal(new_repr_traj, new_repr_traj_ls)) {
  print("Yes, both are equal!")
}
#> [1] "Yes, both are equal!"

define_retra() returns an object of class RETRA that you can use to plot the new trajectories:

plot(new_repr_traj, # <-- This is the RETRA object returned by define_retra()
     # data to generate the state space
     d = state_mds, trajectories = rep(ID_sampling, each = 3), states = rep(ID_obs, 10),
     # display individual trajectories in light blue
     traj.colors = "lightblue",
     # display representative trajectories in dark blue    
     RT.colors = "darkblue",
     # select T1.2 to be displayed in a different color (red)
     select_RT = "T1.2", sel.color = "red",
     # Identify artificial links using dashed lines (default), but use the same 
     # color than the representative trajectories (default)
     link.lty = 2, link.color = NULL,
     # We can use other arguments in plot()
     main = "Defined representative trajectories")

# Add a legend 
legend("bottomright", c("T1.1, T2.1", "T1.2", "Link"), 
       col = c("darkblue", "red", "darkblue"), lty = c(1, 1, 2))

4.4. Variation of ecological properties along representative trajectories

The definition of representative trajectories is based on the real states of the sampling units. Thus, it is possible to assess the variation of some ecological properties along the representative trajectories, and therefore, across the EDR.

Let’s see how species diversity varies along the representative trajectories returned by retra_edr(). First, we need to calculate an index of species diversity (e.g., Shannon index) for the states in each representative trajectory and set the order of the states in the trajectory:

# Set an ID in the abundance matrix
abundance$ID <- paste0(abundance$sampling_unit, abundance$time)

# Calculate the Shannon index in all trajectory states
abundance$Shannon <- vegan::diversity(abundance[, c("sp1", "sp2")], index = "shannon")

# Identify the states forming both representative trajectories
traj_states <- lapply(repr_traj, function(iRT){
  segments <- iRT$Segments
  seg_components <- strsplit(gsub("\\]", "", gsub("\\[", "-", segments)), "-")
  traj_states <- vapply(seg_components, function(iseg){
    c(paste0(iseg[1], iseg[2]), paste0(iseg[1], iseg[3]))
  }, character(2))
  traj_states <- unique(as.character(traj_states))
  traj_states <- data.frame(ID = traj_states, RT_states = 1:length(traj_states))
})

# Associate the states of the representative trajectories with their corresponding
# value of the Shannnon index
RT_diversity <- lapply(traj_states, function(iRT){
  data <- merge(iRT, abundance, by = "ID", all.x = T)
  data <- data[order(data$RT_states), ]
})

RT_diversity$T2
#>   ID RT_states sampling_unit time       sp1        sp2    Shannon
#> 1 C1         1             C    1 0.3761468 0.62385321 0.66214624
#> 2 C2         2             C    2 0.5930319 0.40696810 0.67573602
#> 3 C3         3             C    3 0.7754031 0.22459689 0.53266482
#> 6 I2         4             I    2 0.8009390 0.19906100 0.49909794
#> 7 I3         5             I    3 0.9044378 0.09556219 0.31522122
#> 4 E2         6             E    2 0.9574203 0.04257971 0.17605768
#> 5 E3         7             E    3 0.9809725 0.01902749 0.09422978

Then, we can evaluate the variation of species diversity along the representative trajectories graphically:

# Plot the variation of species diversity in T1
plot(x = RT_diversity$T1$RT_states + 1, y = RT_diversity$T1$Shannon,
     type = "l", col = "blue", xlim = c(1, 7), ylim = c(0, 0.7), 
     xlab = "RT state", ylab = "Shannon index",
     main = "Variation of species diversity")
# Add the variation of species diversity in T2
lines(x = RT_diversity$T2$RT_states, y = RT_diversity$T2$Shannon,
      col = "black")
# Add a legend
legend("topright", c("T1", "T2"), col = c("blue", "black"), lty = 1)

5.1. Dynamic dispersion (dDis)

Dynamic dispersion (dDis) is calculated as the average dissimilarity between the trajectories in an EDR and another trajectory taken as a reference (Sánchez-Pinillos et al., 2023).

If the trajectory taken as the reference is one of the representative trajectories (e.g. T2), dDis can be used as an indicator of its fitness to the data and as an overall metric of the dispersion of the trajectories in the EDR.

First, you need to prepare the data. Even if the states of the representative trajectories are already in your abundance and dissimilarity matrices, you will need to duplicate them to indicate that this is a different trajectory. We can directly use the data frame that we produced to calculate species diversity.

# Change the trajectory identifier and the number of the state
abundance_T2 <- RT_diversity$T2
abundance_T2$sampling_unit <- "T2"
abundance_T2$time <- abundance_T2$RT_states

# Add representative trajectories' data to the abundance matrix
abundance_T2 <- rbind(abundance, abundance_T2[, names(abundance)])

# Calculate state dissimilarities including the representative trajectory
state_dissim_T2 <- vegan::vegdist(abundance_T2[, c("sp1", "sp2")], method = "bray")

Now, we can compute dDis using dDis():

# Compute dDis taking T2 as reference
dDis_T2 <- dDis(d = state_dissim_T2, d.type = "dStates", 
         trajectories = abundance_T2$sampling_unit, states = abundance_T2$time, 
         reference = "T2")
dDis_T2
#> dDis (ref. T2) 
#>      0.1463935

Instead of using the representative trajectory as the reference, we could use any other trajectory of the EDR to quantify how “immersed” it is in the EDR. Let’s compare dDis for trajectories of the sampling units I, which is part of both representative trajectories, and F, which was displayed relatively far from the rest of trajectories in the state and trajectory spaces.

In this case, we will use the trajectory dissimilarity matrix generated for the original data:

# dDis: reference C
dDis_I <- dDis(d = traj_dissim, d.type = "dTraj", 
         trajectories = ID_sampling, 
         reference = "I")
# dDis: reference F
dDis_F <- dDis(d = traj_dissim, d.type = "dTraj", 
         trajectories = ID_sampling, 
         reference = "F")

dDis_I
#> dDis (ref. I) 
#>     0.1847822
dDis_F
#> dDis (ref. F) 
#>     0.5809146

As expected, the dispersion when trajectory F is taken as the reference is larger than for trajectory I, indicating that F is more isolated than I within the EDR.

Alternatively, you can weight each trajectory depending on some criteria. Currently, dDis() is able to compute dDis weighting ecological trajectories by their size (w.type = "size") and length (w.type = "length"). However, if you want to use different criteria, you can use a set of weights of your choice (w.type = "precomputed") and indicate them through the w.values argument. For example, imagine that you want to weight each trajectory by the initial abundance of sp2:

# Define w values
initial_sp2 <- abundance[which(abundance$time == 1), ]$sp2

# Identify the index of the reference trajectories
ind_I <- which(ID_sampling == "I")
ind_F <- which(ID_sampling == "F")

# Compute dDis with weighted trajectories:
# Considering I as reference
dDis_I_w <- dDis(d = traj_dissim, d.type = "dTraj", 
         trajectories = ID_sampling, reference = "I", 
         w.type = "precomputed", w.values = initial_sp2[-ind_I])
# Considering F as reference
dDis_F_w <- dDis(d = traj_dissim, d.type = "dTraj", 
         trajectories = ID_sampling, reference = "F", 
         w.type = "precomputed", w.values = initial_sp2[-ind_F])
dDis_I_w
#> dDis (ref. I) 
#>     0.2517555
dDis_F_w
#> dDis (ref. F) 
#>     0.4742493

In comparison with the previous results, dDis has increased and decreased for I and F, respectively, as a result of the greater abundance of sp2 in F than in I. This result could be interpreted as a lower isolation of F in relation to the trajectories initially dominated by sp2.

Weighting ecological trajectories is equivalent to modifying the location of the trajectory taken as the reference. Therefore, I do not recommend you to weight trajectories if the representative trajectory is the reference of dDis because it could make difficult the interpretation of the results.

5.2. Dynamic beta diversity (dBD)

Dynamic beta diversity (dBD) was proposed by De Cáceres et al. (2019; Ecol. Monogr., https://doi.org/10.1002/ecm.1350) to quantify the overall variability of a group of trajectories. Conceptually, it measures the average distance to the centroid of the trajectories. Therefore, it is quite similar to dDis when it is calculated taking a representative trajectory as the reference.

dBD is implemented in the function dBD(). As you can imagine, dBD in our EDR is very small because the dynamical patterns of all trajectories are very similar.

# Calculate dBD
dBD(d = state_dissim, d.type = "dStates", 
    trajectories = rep(ID_sampling, each = 3), states = rep(ID_obs, 10))
#>        dBD 
#> 0.06252172

5.3. Dynamic evenness (dEve)

Dynamic evenness (dEve) quantifies the regularity with which an EDR is filled by ecological trajectories (Sánchez-Pinillos et al., 2023). This metric informs about the existence of groups of trajectories within the EDR.

You can compute dEve with the function dEve():

# Calculate dEve
dEve(d = traj_dissim, d.type = "dTraj", trajectories = ID_sampling)
#>      dEve 
#> 0.5827095

Like dDis, dEve can also be calculated using different weights for the individual trajectories of the EDR. In that case, dEve is equivalent to the functional evenness index proposed by Villéger et al. (2003; Ecology, https://doi.org/10.1890/07-1206.1).

For example, weighting ecological trajectories by the initial abundance of sp2 is equivalent to reducing the dissimilarity of the trajectories dominated by sp2 and dEve increases.

# Calculate dEve weighting trajectories by the initial abundance of sp2
dEve(d = traj_dissim, d.type = "dTraj", trajectories = ID_sampling,
         w.type = "precomputed", w.values = initial_sp2)
#>      dEve 
#> 0.6296034