---
title: "Predicting ecological trajectories from ecological dynamic regimes"
author: "Martina Sánchez-Pinillos"
output: 
  rmarkdown::html_vignette:
    toc: TRUE
vignette: >
  %\VignetteIndexEntry{Predicting ecological trajectories from ecological dynamic regimes}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = requireNamespace("vegan", quietly = TRUE)
)
```

```{r vegan, eval = !requireNamespace("vegan", quietly = TRUE), echo = FALSE, comment = NA}
message('This vignette requires the "vegan" package. Please, install it: install.packages("vegan").')
```

## 1. Introduction

### 1.1. Predicting ecological trajectories within the EDR framework

**Ecological dynamic regimes (EDRs)** are defined as a group of trajectories reflecting the temporal changes in a set of state variables under certain external conditions and in the absence of perturbations. As such, EDRs provide valuable information to infer the dynamics of analogous systems (i.e., defined by the same state variables) for which we lack long temporal series.

**PETRA-EDR (_Predicted Ecological TRAjectories in Ecological Dynamic Regimes_**; [Sánchez-Pinillos et al., 2026](https://doi.org/10.1111/2041-210x.70372)) is an algorithm implemented in `ecoregime` that leverages ecological dynamic regimes defined in a multidimensional state space to forecast the unknown trajectory of any ecosystem state. PETRA-EDR is based on the rationale that the successional changes in the state variables of ecological systems under comparable and relatively stable external conditions follow similar patterns. PETRA-EDR does not aim to outperform other data-driven and process-based models to forecast ecological dynamics, but to provide the EDR framework with predictive capacity, improving ecological resilience analyses, and enabling the detection of dynamic regime shifts when long time series are unavailable.

In particular, PETRA-EDR can be used to predict the dynamics of a disturbed system if that disturbance had not occurred. That expected undisturbed trajectory can then be used as the reference to quantify the deviation of the disturbed trajectory through amplitude, recovery, and net change indicators and assess the system's ecological resilience (see `vignette("Resilience")`).

Moreover, forecasting the post-disturbance dynamics of a system within a stability landscape composed of alternative EDRs can be useful to identify potential regime shifts.

PETRA-EDR is complemented by the metric **MPD (*Mean Predicted Deviation*)**, which quantifies the prediction accuracy of PETRA-EDR outputs assuming that, in the absence of stochasticity and observational noise, there is one and only one ecological trajectory passing by the target.

To learn more about the technical details of PETRA-EDR and MPD, see this publication:

- Sánchez-Pinillos M., Fortin, M-J., Messier, C., Kneeshaw, D. 2026. Forecasting ecological trajectories from ecological dynamic regimes to improve resilience analysis. *Methods in Ecology and Evolution*. <https://doi.org/10.1111/2041-210x.70372>

Additional information about the analysis of ecological dynamic regimes can be found in `vignette("EDR_framework")` and this publication:

- 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>

To assess ecological resilience using EDRs, you can check `vignette("Resilience")` and this publication:

- Sánchez-Pinillos M., Dakos, V., Kéfi, S. 2024. Ecological dynamic regimes: A   key concept for assessing ecological resilience. *Biological Conservation*. <https://doi.org/10.1016/j.biocon.2023.110409>

### 1.2. About this vignette

This vignette focuses on the PETRA-EDR algorithm proposed in [Sánchez-Pinillos et al., (2026)](https://doi.org/10.1111/2041-210x.70372) to predict the ecological trajectories from ecological dynamic regimes. This algorithm was implemented in the function `petra_edr()` of the R package `ecoregime`. In particular, this vignette introduces the arguments and outputs of the function `petra_edr()`, summarizes the workflow to predict trajectories, quantifies the prediction accuracy through the metric MPD implemented in `MPD()`, and represents predicted trajectories in a state space using the function `plot()`.

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

```{r install, eval=FALSE}
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:

```{r setup}
library(ecoregime)
```

```{r citation}
citation("ecoregime")
```

## 2. Artificial data

First, we will generate three hypothetical targets from the EDR data included in `ecoregime`. We will use them along the vignette to predict their trajectories using EDR1 as the reference.

Although the data included in `ecoregime` refer to ecological communities and species abundances, it is important to note that other ecological units (e.g., individuals, populations) and state variables (e.g., number of individuals, functional traits) can be used. Some examples are detailed in the Supporting Information (Appendix S1) of the associated paper ([Sánchez-Pinillos et al., 2026](https://doi.org/10.1111/2041-210x.70372))

```{r data}
# Matrix including the state variables (sp1-sp12) of the EDR trajectories
edr <- EDR_data$EDR1$abundance

# The first target is composed of one state resulting from averaging the state
# variables of two states in the reference EDR
target1 <- data.frame(matrix(colMeans(edr[traj == 3 & state %in% 1:2, paste0('sp', 1:12)]),
                             ncol = 12, 
                             dimnames = list(1, paste0('sp', 1:12))))
target1$traj <- 'target1'
target1$state <- 1

# The second target is composed of three states resulting from averaging the 
# state variables of four states in the reference EDR
target2 <- data.frame(t(sapply(1:3, function(istate){
  matrix(colMeans(edr[traj == 6 & state %in% istate:(istate+1), 
                      paste0('sp', 1:12)]))
})))
names(target2) <- paste0('sp', 1:12)
target2$traj <- 'target2'
target2$state <- 1:3

# For the third target, we will consider a trajectory of a different EDR
target3 <- EDR_data$EDR2$abundance[1:5, 3:ncol(EDR_data$EDR2$abundance)]
target3$traj <- 'target3'

```

## 3. Understanding `petra_edr()`

### 3.1. The arguments of `petra_edr()`

We can classify the arguments of `petra_edr()` into four main groups:

- Arguments to define the state variables and identify trajectories and states: *state_var, trajectories, states, targets*

- Arguments to define the state space: *d_function, d_args, d*

- Parameters of PETRA-EDR: *k, eps, minPts, w_function, alpha, w, method, direction*

- Arguments to define the contents of the output: *return_args*

Let's use the first target defined in the previous section as an example to better understand the arguments of `petra_edr()`.

- **Arguments to define the state variables and identify EDR trajectories, states, and targets:**

  - `state_var`: the state variables of each EDR and target state. In our example, `state_var` is a `data.frame` including the species abundances for all states in the EDR and the target:

```{r state_var}
# Select the columns containing the state variables (sp1, ..., sp12) and 
# include the information of EDR and target states in the same data.frame 
state_var1 <- data.frame(rbind(edr[, paste0("sp", 1:12)], 
                               target1[, paste0("sp", 1:12)]))
head(state_var1)
```

- `trajectories`: the ID of all trajectories or sites, including both the EDR trajectories and the targets.

```{r trajectories1}
trajectories1 <- c(edr$traj, target1$traj)
head(trajectories1)
```

- `states`: the ID of all states or surveys, including both the EDR states and the target states.

```{r states1}
states1 <- as.integer(c(edr$state, target1$state))
head(states1)
```

- `targets`: the ID of the targets for which we want to predict their dynamics. In our example, `"target1"`.

- **Arguments to define the state space:** In the EDR framework, the state space is defined by a dissimilarity matrix. `petra_edr()` needs to recalculate the state space as it advances predicting new states. As a consequence, we need to indicate the function used to define the state space and a list of arguments to apply the dissimilarity function:

  - `d_function`: name of the function used to compute the state dissimilarity matrix in the form `"package::function"` (e.g., `"vegan::vegdist"`). It is important to consider a function returning an object of class `dist`.

  - `d_args`: a list of arguments needed to apply `d_function`. This list depends on `d_function`, for example, in the function `vegdist` in `vegan`, the minimum argument that we need to specify is `x` (i.e., the community data matrix, which in our case is `state_var`). If we want to modify the other arguments of the function, we will have to include it in `d_args`. For example, to apply the Bray-Curtis dissimilarity, we need to include de argument `method` in `d_args`. Thus, in this example, we would define `d_args = list(x = state_var, method = 'bray')`.

  - `d`: we can include a matrix containing the dissimilarities between each pair of states in the EDR and targets in the same order than the specified in `trajectories` and `states`. In our example, `d` can be calculated by applying `d_function` to `state_var`. If `d` is not specified, `petra_edr()` will compute it using the arguments `d_function` and `d_args`.

```{r dStates1}
# Compute state dissimilarities from state_var
dStates1 <- vegan::vegdist(x = state_var1, method = "bray")
```

- **Parameters of PETRA-EDR**: These parameters are specific of the PETRA-EDR algorithm. You can find useful information in [Sánchez-Pinillos et al., (2026)](https://doi.org/10.1111/2041-210x.70372) and its Supporting Information (Appendix S1).

  - `k`: number of the nearest states to the target. Small values of `k` returns forecasts highly dependent on the nearest trajectories, whereas large values consider the trends of many trajectories in the EDR. Let's see the differences of applying `petra_edr()` in our example when `k = 2L` and `k = 20L`:
    
```{r petra k, fig.width=8, fig.height=4, warning=FALSE}
# Compute petra_edr using a small k
petra_k1 <- petra_edr(state_var = state_var1, 
                      trajectories = trajectories1,
                      states = states1, 
                      targets = "target1",
                      d_function = "vegan::vegdist", 
                      d_args = list(x = state_var1, method = "bray"),
                      k = 2L,
                      minPts = 2L,
                      return_args = T)

# Compute petra_edr using a large k
petra_k2 <- petra_edr(state_var = state_var1, 
                      trajectories = trajectories1,
                      states = states1, 
                      targets = "target1",
                      d_function = "vegan::vegdist", 
                      d_args = list(x = state_var1, method = "bray"),
                      k = 20L,
                      minPts = 2L,
                      return_args = T)

# Use plot to see the results
par(mfrow = c(1, 2))
plot(x = petra_k1, 
     target.colors = "red", 
     petra.colors = "blue", 
     xlab = "MDS D1", 
     ylab = "MDS D2", 
     main = "k = 2")
legend("bottomleft", 
       c("Target", "Predicted trajectory", "EDR trajectories"), 
       lwd = 2,
       col = c("red", "blue", "grey"), 
       cex = 0.8, 
       bty = "n")

plot(x = petra_k2, 
     target.colors = "red", 
     petra.colors = "blue", 
     xlab = "MDS D1", 
     ylab = "MDS D2", 
     main = "k = 20")
legend("bottomleft", 
       c("Target", "Predicted trajectory", "EDR trajectories"), 
       lwd = 2, 
       col = c("red", "blue", "grey"), 
       cex = 0.8, 
       bty = "n")
```
    

When `k = 2L`, `petra_edr()` selects two states that are very close to the target:

```{r petra k1}
petra_k1$k_dist
```

In contrast, when `k = 20L`, `petra_edr()` selects 20 states being the furthest state at a dissimilarity larger than 0.07:

```{r petra k2}
tail(petra_k2$k_dist)
```

While setting `k=20L` returns a longer trajectory, the predicted trajectory has a greater uncertainty, since it can be biased by the furthest states.

- `eps`: dissimilarity threshold beyond which EDR states are disregarded in the computation of the predicted trajectories. This argument can be used to avoid biased outputs due to far states.

```{r petra eps, fig.width=8, fig.height=4, warning=FALSE}
# Compute petra_edr using a small eps
petra_eps1 <- petra_edr(state_var = state_var1, 
                        trajectories = trajectories1, 
                        states = states1,
                        targets = "target1",
                        d_function = "vegan::vegdist", 
                        d_args = list(x = state_var1, method = "bray"),
                        k = 20L, 
                        minPts = 2L, 
                        eps = 0.03, 
                        return_args = T)

# Compute petra_edr using a large eps
petra_eps2 <- petra_edr(state_var = state_var1, 
                        trajectories = trajectories1, 
                        states = states1, 
                        targets = "target1",
                        d_function = "vegan::vegdist", 
                        d_args = list(x = state_var1, method = "bray"),
                        k = 20L, 
                        minPts = 2L, 
                        eps = 0.1, 
                        return_args = T)

# Plot PETRA outputs
par(mfrow = c(1, 2))
plot(x = petra_eps1, 
     target.colors = "red", 
     petra.colors = "blue", 
     xlab = "MDS D1", 
     ylab = "MDS D2", 
     main = "eps = 0.03")
legend("bottomleft", 
       c("Target", "Predicted trajectory", "EDR trajectories"), 
       lwd = 2,
       col = c("red", "blue", "grey"), 
       cex = 0.8, 
       bty = "n")

plot(x = petra_eps2, 
     target.colors = "red", 
     petra.colors = "blue", 
     xlab = "MDS D1", 
     ylab = "MDS D2", 
     main = "eps = 0.1")
legend("bottomleft", 
       c("Target", "Predicted trajectory", "EDR trajectories"), 
       lwd = 2,
       col = c("red", "blue", "grey"), 
       cex = 0.8, 
       bty = "n")
         
```

Even if we specify `k = 20L`, we can restrict the number of selected states to those within a radius `eps` from the target states. For example, when we set `eps = 0.03`, only three out of the 20 states are used in the analyses:

```{r petra eps1}
petra_eps1$k_dist
```

- `minPts`: minimum number of states required to calculate the predicted states. In general, smaller values of `minPts` lead to longer predicted trajectories. However, the states calculated from a small number of states could be associated with a higher uncertainty.

```{r petra minPts, fig.width=8, fig.height=4, warning=FALSE}
# Compute petra_edr using a small minPts
petra_minPts1 <- petra_edr(state_var = state_var1, 
                           trajectories = trajectories1, 
                           states = states1,
                           targets = "target1",
                           d_function = "vegan::vegdist", 
                           d_args = list(x = state_var1, method = "bray"),
                           k = 6L, 
                           minPts = 2L, 
                           return_args = T)

# Compute petra_edr using a large minPts
petra_minPts2 <- petra_edr(state_var = state_var1, 
                           trajectories = trajectories1, 
                           states = states1,
                           targets = "target1",
                           d_function = "vegan::vegdist", 
                           d_args = list(x = state_var1, method = "bray"),
                           k = 6L, 
                           minPts = 6L, 
                           return_args = T)

# Plot PETRA outputs
par(mfrow = c(1, 2))
plot(x = petra_minPts1, 
     target.colors = "red", 
     petra.colors = "blue", 
     xlab = "MDS D1", 
     ylab = "MDS D2", 
     main = "minPts = 2")
legend("bottomleft", 
       c("Target", "Predicted trajectory", "EDR trajectories"), 
       lwd = 2,
       col = c("red", "blue", "grey"), 
       cex = 0.8, 
       bty = "n")

plot(x = petra_minPts2, 
     target.colors = "red", 
     petra.colors = "blue", 
     xlab = "MDS D1", 
     ylab = "MDS D2", 
     main = "minPts = 10")
legend("bottomleft", 
       c("Target", "Predicted trajectory", "EDR trajectories"), 
       lwd = c(NA, 2, 2), 
       pch = c(20, NA, NA),
       col = c("red", "blue", "grey"), 
       cex = 0.8, 
       bty = "n")
```

When we set `minPts = 2L`, `petra_edr()` returns five predicted states, two of
them calculated from four and two EDR states, respectively. In contrast, by
setting `minPts = 6L`, we ensured that all predicted states are calculated from
six or more EDR states, indicating a higher probability of finding those state
variables in nature:

```{r petra minPts1 minPts2}
petra_minPts1$predicted_dist[, c("target", "predicted_state", "N")]
petra_minPts2$predicted_dist[, c("target", "predicted_state", "N")]
```

- `w_function`: instead of using the argument `eps`, we can assign different weights as a function of the dissimilarity of the target to the EDR states to improve the results. `petra_edr()` permits using six different weighting functions: `"linear"`, `"power"`, `"exponential"`, `"Gaussian"`, `"hyperbolic"`, or `"spherical"`. You can see the relevance of this argument in [Sánchez-Pinillos et al., (2026)](https://doi.org/10.1111/2041-210x.70372).
- `alpha`: when `w_function` is one of `"power"`, `"exponential"`, or `"hyperbolic"`, `alpha` is the parameter modulating the shape of these functions.
- `w`: it is possible to assign different weights using some values pre-defined by the user. For example, one could be interested in defining `w` depending on some external variables.
- `method`: it indicates how the predicted states are calculated.

```{r method}
# Compute petra_edr using method = "mean"
petra_method1 <- petra_edr(state_var = state_var1, 
                           trajectories = trajectories1, 
                           states = states1,
                           targets = "target1",
                           d_function = "vegan::vegdist", 
                           d_args = list(x = state_var1, method = "bray"),
                           k = 20L, 
                           minPts = 2L, 
                           method = "mean",
                           return_args = T)

# Compute petra_edr using method = "medoid"
petra_method2 <- petra_edr(state_var = state_var1, 
                           trajectories = trajectories1, 
                           states = states1, 
                           targets = "target1",
                           d_function = "vegan::vegdist", 
                           d_args = list(x = state_var1, method = "bray"),
                           k = 20L, 
                           minPts = 2L, 
                           method = "medoid",
                           return_args = T)
```

If `method = "mean"`, the predicted states are calculated averaging the state variables of the EDR states. If `method = "medoid"`, the predicted states coincide with the medoids of the EDR states. We can zoom in one of the predicted states to see the differences between both methods:

```{r method1, fig.width=8, fig.height=4, warning=FALSE}
par(mfrow = c(1, 2))
plot(petra_method1, 
     target.colors = "red", 
     petra.colors = "blue", 
     xlim = c(-0.1, 0.3), 
     ylim = c(0.1, 0.4),
     xlab = "MDS D1", 
     ylab = "MDS D2", 
     main = "method = mean")

plot(petra_method2, 
     target.colors = "red", 
     petra.colors = "blue", 
     xlim = c(-0.1, 0.3), 
     ylim = c(0.1, 0.4),
     xlab = "MDS D1", 
     ylab = "MDS D2", 
     main = "method = medoid")
```

- `direction`: we can specify whether we want to predict the dynamics before and/or after the target.

```{r direction, fig.width=9, fig.height=3, warning=FALSE}
# Compute petra_edr before the target
petra_direction1 <- petra_edr(state_var = state_var1, 
                              trajectories = trajectories1, 
                              states = states1, 
                              targets = "target1",
                              d_function = "vegan::vegdist", 
                              d_args = list(x = state_var1, method="bray"),
                              k = 20L, 
                              minPts = 2L, 
                              direction = -1, 
                              return_args = T)

  # Compute petra_edr after the target
petra_direction2 <- petra_edr(state_var = state_var1, 
                              trajectories = trajectories1, 
                              states = states1, 
                              targets = "target1",
                              d_function = "vegan::vegdist", 
                              d_args = list(x = state_var1, method="bray"),
                              k = 20L, 
                              minPts = 2L, 
                              direction = 1, 
                              return_args = T)

  # Compute petra_edr before and after the target
petra_direction3 <- petra_edr(state_var = state_var1, 
                              trajectories = trajectories1, 
                              states = states1, 
                              targets = "target1",
                              d_function = "vegan::vegdist", 
                              d_args = list(x = state_var1, method="bray"),
                              k = 20L, 
                              minPts = 2L, 
                              direction = 2,
                              return_args = T)

  # Plot PETRA outputs
par(mfrow = c(1, 3))
plot(x = petra_direction1, 
     target.colors = "red", 
     petra.colors = "blue", 
     xlab = "MDS D1", 
     ylab = "MDS D2", 
     main = "direction = -1")
legend("bottomleft", 
       c("Target", "Predicted trajectory", "EDR trajectories"), 
       lwd = 2,
       col = c("red", "blue", "grey"), 
       cex = 0.8, 
       bty = "n")

plot(x = petra_direction2, 
     target.colors = "red", 
     petra.colors = "blue", 
     xlab = "MDS D1", 
     ylab = "MDS D2", 
     main = "direction = 1")
legend("bottomleft", 
       c("Target", "Predicted trajectory", "EDR trajectories"), 
       lwd = c(NA, 2, 2), 
       pch = c(20, NA, NA),
       col = c("red", "blue", "grey"), 
       cex = 0.8, 
       bty = "n")

plot(x = petra_direction3, 
     target.colors = "red", 
     petra.colors = "blue", 
     xlab = "MDS D1", 
     ylab = "MDS D2", 
     main = "direction = 2")
legend("bottomleft", 
       c("Target", "Predicted trajectory", "EDR trajectories"), 
       lwd = 2,
       col = c("red", "blue", "grey"), 
       cex = 0.8, 
       bty = "n")
```

- **Argument to define the contents of the output**

  - `return_args`: if `TRUE`, the output contains the information included in the arguments of `petra_edr()`. This is essential for computing other functions, like `plot()` or `MPD()`.

### 3.2. The outputs of `petra_edr()`

`petra_edr()` returns a list of elements including the state variables of the predicted trajectory and some statistics associated with the targets and the predicted states.

- `state_var`: element of the same class than the argument `state_var` including the state variables for all states in the predicted trajectory.

```{r state_var2}
petra_k2$state_var
```

- `trajectories` and `states`: these elements are useful to identify the trajectory and state ID of the predicted trajectory in `state_var`.

In our example, as we only used one target, `trajectories` indicate that all rows in `state_var` refer to target1. This output is, however, useful when we include several targets in `petra_edr()` as we will see later.

```{r trajectories2}
petra_k2$trajectories
```

`states` shows the order of the predicted states. In this case, target1 contains only one state, which is identified by `1`. The states `-3:0` and `2:5` correspond to the predicted states preceding and following the target, respectively.

```{r states2}
petra_k2$states
```

`trajectories` and `states` serve to identify the states in `state_var`:

```{r state_var trajectories states}
state_var1_k2 <- petra_k2$state_var
state_var1_k2$traj <- petra_k2$trajectories
state_var1_k2$state <- petra_k2$states
state_var1_k2[, c("traj", "state", paste0("sp", 1:12))]
```

- `k_dist`: data.table containing the ID and dissimilarity of the k-nearest states in the EDR to each extreme of the target.

For example, when we set `k = 20L` and `eps = 0.03`, `petra_edr` selected three EDR states: the state 2 of trajectory 9, the state 2 of trajectory 16, and the state 2 of trajectory 10. `k_dist` also includes the dissimilarity between the target states and each of the k-nearest EDR states.

```{r k_dist}
petra_eps1$k_dist
```

- `predicted_dist`: data table including some statistics associated with each of the predicted states.

Following the previous example (`k = 2L` and `eps = 0.03`), the predicted trajectory for target1 contains four predicted states (0, 2, 3, 4). `N` includes the number of EDR states used to calculate the state variables of each predicted state (3 in all predicted states); `mean_dist`, `sd_dist`, `min_dist`, and `max_dist` respectively include the mean, standard deviation, minimum, and maximum dissimilarities between the predicted state and the EDR states used to calculate their state variables.

The states used to calculate the state variables of the predicted state 0 are, in average further than those used for the other predicted states, showing the highest values in all statistics and, therefore, indicating a greater uncertainty in its prediction. Note that these statistics depend on the dissimilarity metric used to define the state space (i.e., `d_function`) and, therefore, the uncertainty of the predicted state must be interpreted accordingly.

```{r predicted_dist}
petra_eps1$predicted_dist
```

## 4. Forecasting ecological trajectories

### 4.1. Preliminary analyses

Before applying `petra_edr()` to predict the ecological trajectories of our three targets, it is useful to carry out some preliminary analyses that will inform us about the adequacy of our data to apply PETRA-EDR.

The function `plot_edr()` can be useful for evaluating the position of the targets relative to the EDR. For that, we first need to define the state space using a dissimilarity matrix including all EDR trajectories and targets. In our example, we will use the Bray-Curtis dissimilarity, but other metric can be more adequate depending on your data.

```{r dStates}
# State variables and dissimilarity metric of the EDR and the first target
state_var1 <- rbind(edr[, paste0("sp", 1:12)], 
                    target1[, paste0("sp", 1:12)])
dStates1 <- vegan::vegdist(state_var1, method = "bray")

# State variables and dissimilarity metric of the EDR and the second target
state_var2 <- rbind(edr[, paste0("sp", 1:12)], 
                    target2[, paste0("sp", 1:12)])
dStates2 <- vegan::vegdist(state_var2, method = "bray")

# State variables and dissimilarity metric of the EDR and the third target
state_var3 <- rbind(edr[, paste0("sp", 1:12)], 
                    target3[, paste0("sp", 1:12)])
dStates3 <- vegan::vegdist(state_var3, method = "bray")

```

Now we can use the function `plot_edr()` to visualize the position of each target in the EDR:

```{r plot_edr, fig.width=9, fig.height=3, warning=FALSE}
par(mfrow = c(1, 3))

# Number of trajectories in the EDR
Ntraj <- length(unique(edr$traj))

# Location of target1. 
# As target1 is composed of one state, we need to specify type = 'states'
plot_edr(x = dStates1, 
         trajectories = c(edr$traj, target1$traj),
         states = as.integer(c(edr$state, target1$state)),
         type = "states",
         state.colors = c(rep("grey", length(edr$traj)), "red"),
         xlab = "MDS D1", 
         ylab = "MDS D2", 
         main = "Location of the first target")
legend("bottomleft", 
       c("EDR trajectories", "Target 1"), 
       lwd = c(2, NA), 
       pch = c(NA, 20), 
       col = c("grey", "red"), 
       cex = 0.8, 
       bty = "n")


# Location of target2
plot_edr(x = dStates2, 
         trajectories = c(edr$traj, target2$traj),
         states = as.integer(c(edr$state, target2$state)),
         traj.colors = c(rep("grey", Ntraj), "red"),
         xlab = "MDS D1", 
         ylab = "MDS D2", 
         main = "Location of the second target")
legend("bottomleft", 
       c("EDR trajectories", "Target 2"), 
       lwd = 2, 
       col = c("grey", "red"), 
       cex = 0.8, 
       bty = "n")

# Location of target3
plot_edr(x = dStates3, 
         trajectories = c(edr$traj, target3$traj),
         states = as.integer(c(edr$state, target3$state)),
         traj.colors = c(rep("grey", Ntraj), "red"),
         xlab = "MDS D1", 
         ylab = "MDS D2", 
         main = "Location of the third target")
legend("bottomleft", 
       c("EDR trajectories", "Target 3"), 
       lwd = 2, 
       col = c("grey", "red"), 
       cex = 0.8, 
       bty = "n")
```

Whereas the first and second targets are within the boundaries of the EDR, the third target is far from the EDR trajectories. This indicates that the forecast of the third target will not be accurate.

We can get similar conclusions by computing the dynamic dispersion (*dDis*) of each target relative to the EDR trajectories:

```{r dDis}
# We cannnot compute dDis for the first target because it is composed of a single state. We will calculate an equivalent metric by calculating the average dissimilarity between the target and the EDR states.
dDis1 <- mean(as.matrix(dStates1)[-nrow(state_var1), nrow(state_var1)])
names(dDis1) <- "dDis (ref. target1)"

# For targets 2 and 3, we use the function dDis
dDis2 <- dDis(d = dStates2, 
              d.type = 'dStates', 
              trajectories = c(edr$traj, target2$traj),
              states = c(edr$state, target2$state), 
              reference = 'target2')
dDis3 <- dDis(d = dStates3, 
              d.type = 'dStates', 
              trajectories = c(edr$traj, target3$traj),
              states = c(edr$state, target3$state), 
              reference = 'target3')

dDis1; dDis2; dDis3

```

**DISCLAIMER:** to predict the trajectories of target3, we should either use a different EDR as the reference or discard the analyses for the inadequacy of the data. However, we will perform the analyses so we can compare the results.

### 4.2. Forecasting ecological trajectories of multiple targets using different conditions

We can predict the ecological trajectories of multiple targets using a common EDR.

For that, `state_var` must include the state variables of the reference EDR and the targets:

```{r state_var3}
# data.table including the state variables, trajectories, and states of the EDR and the targets
data <- rbind(edr[, -1], target1, target2, target3)

# state_var needs to be a data.frame with only the state variables
state_var <- data.frame(data[, paste0("sp", 1:12)])
```

We can define different conditions for each target:

```{r target_pars, echo=FALSE}
target_pars <- data.frame(Target = paste0("target", 1:3),
                          k = c(10L, 50L, 50L),
                          eps = c(NA, 0.05, 0.5), 
                          minPts = c(3L, 2L, 2L),
                          w_function = c(NA, "exponential", "linear"), 
                          alpha = c(NA, 3, NA))

knitr::kable(target_pars, row.names = F)
```

Some arguments (`d_function`, `method`, `direction`) need to be equal for all the targets. Thus, we will use `method = "mean"` and Bray-Curtis dissimilarities. We need to set `return_args = TRUE` to be able to assess the prediction quality and represent the predicted trajectories later.

```{r petra}
petra <- petra_edr(state_var = state_var,
                   trajectories = data$traj,
                   states = as.integer(data$state),
                   targets = c("target1", "target2", "target3"),
                   d_function = "vegan::vegdist", 
                   d_args = list(x = state_var, method = "bray"),
                   k = c(10L, 50L, 50L), 
                   eps = c(NA, 0.05, 0.5), 
                   minPts = c(3L, 2L, 2L),
                   w_function = c(NA, "exponential", "linear"), 
                   alpha = c(NA, 3, NA),
                   direction = 2, 
                   method = "mean", 
                   return_args = T)
```

### 4.3. Evaluating the predicted trajectories

#### The k-nearest EDR states

The predicted trajectory of **target1** was calculated from the 10 nearest states in the EDR trajectories, which were located at a dissimilarity that ranged between 0.025 and 0.054.

```{r petra k_dist 1}
petra$k_dist[target == "target1"]
```

For **target2**, although we considered the 50 nearest states, the algorithm restricted the analyses to the 17 states that were closer than 0.05

```{r petra k_dist 2}
petra$k_dist[target == "target2"]
```

In the case of **target3**, despite being less restrictive (`eps = 0.5`), the analyses were conducted from only nine states, with dissimilarities ranging between 0.467 and 0.495. The large dissimilarities of the k-nearest states indicates that the predicted trajectory of target3 is not very accurate.

```{r petra k_dist 3}
petra$k_dist[target == "target3"]
```

#### The predicted states

The predicted trajectory of **target1** contains six predicted states calculated from a minimum number of four states in the EDR trajectories, with mean dissimilarities ranging between 0.025 and 0.048 and standard deviations between 0.007 and 0.019.

```{r predicted_dist 1}
petra$predicted_dist[target == "target1"]
```

The predicted trajectory of **target2** has eight predicted states calculated from a minimum number of two states in the EDR trajectories, with mean dissimilarities ranging between 0.007 and 0.059 and standard deviations between 0.002 and 0.017. Although the mean dissimilarity of the EDR states to the predicted state -2 is relatively high (0.059), using the exponential weighting function can significantly reduce potential biases due to the furthest states.

```{r predicted_dist 2}
petra$predicted_dist[target == "target2"]
```

Despite the large dissimilarities of the k-nearest EDR states to the **target3**, the mean dissimilarities of the predicted states to the states of the EDR trajectories are much smaller. Although the results in `predicted_dist` are useful to compare the uncertainty of the predicted states in the surroundings of the EDR trajectories, additional analyses are required to assess the accuracy of the predicted trajectory regarding the target states.

```{r predicted_dist 3}
petra$predicted_dist[target == "target3"]
```

## 5. Prediction accuracy

To assess the accuracy of the predicted trajectories, we can use the MPD metric, which quantifies the average dissimilarity between the target states and new forecasts generated from the predicted states. Thus, MPD is expressed in the same units than the dissimilarity metric set in `d_function`.

In contrast to the predicted trajectories for target1 and target2, the large MPD value for target3 confirms the inaccuracy of the predicted trajectory.

```{r MPD}
MPD(x = petra)
```

Besides informing about the prediction accuracy, it is important to note that MPD can be used in optimization techniques to improve our results by adjusting the values of `k`, `eps`, `minPt`, and `alpha` (if necessary). You can see more details in [Sánchez-Pinillos et al., (2026)](https://doi.org/10.1111/2041-210x.70372).

## 6. Representing predicted trajectories

Finally, we can represent the targets, their predicted trajectories, and the EDR trajectories used to compute PETRA-EDR in a common multidimensional state space.

```{r plot predicted, fig.width=5, fig.height=5, warning=FALSE}
plot(x = petra,
     xlab = "MDS D1", 
     ylab = "MDS D2")
legend("bottomleft", 
       c("Predicted trajectories", "EDR trajectories"),
       lwd = 2, 
       col = c("red", "grey"), 
       cex = 0.8, 
       bty = "n")
```

The function `plot()` includes multiple arguments that can be set depending on the elements to be highlighted.

We can represent each predicted trajectory with a different color:

```{r plot predicted col, fig.width=5, fig.height=5, warning=FALSE}
plot(x = petra, 
     petra.colors = grDevices::palette.colors(6, "Paired")[c(2, 4, 6)],
     xlab = "MDS D1", 
     ylab = "MDS D2")
legend("topleft", 
       c("Predicted trajectory 1", "Predicted trajectory 2",
         "Predicted trajectory 3", "EDR trajectories"),
       lwd = 2, 
       col = c(grDevices::palette.colors(6, "Paired")[c(2, 4, 6)], "grey"), 
       cex = 0.8, 
       bty = "n")
```

We can also define different colors for the observed states (i.e., target states) and the predicted states. For example, we could use darker colors for the observed states:

```{r plot targets, fig.width=5, fig.height=5, warning=FALSE}
plot(x = petra, 
     traj.colors = grDevices::palette.colors(9, "Set 3")[9],
     petra.colors = grDevices::palette.colors(6, "Paired")[c(1, 3, 5)],
     target.colors = grDevices::palette.colors(6, "Paired")[c(2, 4, 6)],
     xlab = "MDS D1", 
     ylab = "MDS D2")

legend("bottomleft", 
       c("Target 1", "Predicted trajectory 1", 
         "Target 2", "Predicted trajectory 2",
         "Target 3", "Predicted trajectory 3", 
         "EDR trajectories"),
       col = c(grDevices::palette.colors(6, "Paired")[c(2, 1, 4, 3, 6, 5)],
               grDevices::palette.colors(9, "Paired")[9]),
       lwd = 2, 
       ncol = 2, 
       cex = 0.8, 
       bty = "n")
```

Finally, we can use the results in `predicted_dist` and represent the uncertainty associated with each predicted state. For example, based on the average dissimilarity between the predicted states and the EDR trajectories:

```{r plot mean_dist, fig.width=5, fig.height=5, warning=FALSE}
plot(x = petra, 
     petra.colors = grDevices::hcl.colors(5, "Viridis")[3],
     target.colors = grDevices::hcl.colors(5, "Viridis")[1], 
     uncert.metric = "mean_dist", 
     uncert.colors = grDevices::hcl.colors(5, "Viridis"),
     xlab = "MDS D1", 
     ylab = "MDS D2")

legend("topleft", 
       legend = c(paste0("mean_dist = ", 
                         round(min(petra$predicted_dist$mean_dist), 2)),
                  rep(NA, 18), 
                  paste0("mean_dist = ",
                         round(max(petra$predicted_dist$mean_dist), 2))),
       fill = grDevices::hcl.colors(20, "Viridis"), 
       border = NA, 
       y.intersp = 0.2,
       cex = 0.8, 
       bty = "n")
```