rEDM的用法（5）——Scenario exploration（待续）

（代码有点小问题，还没有全部解决）

Load package and time series

In this demonstration, we use the same time series generated from the 5-species model (Resource, Consumer1, Consumer2, Predator1, Predator2) following Deyle et al. (2016).

## Load package and data
library(rEDM) 
d <- read.csv("ESM5_Data_5spModel.csv")

Set parameters

In this demonstration, we focus on the dynamics of Consumer1 (C1). We investigate how changes in Resource (R) induce changes in C1 using scenario exploration (Deyle et al. 2013). First, we determine the best embedding dimension of C1 using simplex projection

# Normalize data
std_C1 <- as.numeric(scale(d[,"C1"]))

# Determine the best embedding dimension by simplex projection
E_tested <- simplex(std_C1, E=1:10, silent=T)

 # The best E is determined based on MAE
E_C1 <- E_tested[which.min(E_tested$mae),"E"]

The, we generate data frame for scenario exploration by adding one dimension (i.e., R(t)) to the univariate embedding. The multivariate embedding follows:

# Normalize data
std_R <- as.numeric(scale(d[,"R"]))

# Make time-delayed embedding
embed_C1 <- embed(std_C1, dimension = E_C1)

# Add R as an external force
block0 <- cbind(embed_C1, std_R[E_C1:length(std_R)])

We prepare a new "R" column that includes simulated changes in Resource (0.5*sd(R) increase or decrease in Resource). The magnitude of change is arbitrarily defined here. The users can decide the magnitude of change, depending on their own research questions.

# Make new R (increased/decreased situation)
std_R_increase <- std_R + 0.5
std_R_decrease <- std_R - 0.5

# Make new embeddings (increased/decreased situation)
block_inc0 <- cbind(embed_C1, std_R_increase[E_C1:length(std_R)])
block_dec0 <- cbind(embed_C1, std_R_decrease[E_C1:length(std_R)])

# Combine the original and simulated time series
block_inc <- rbind(block0, rep(NaN,4), block_inc0) # NaNs were added to separate two data
block_dec <- rbind(block0, rep(NaN,4), block_dec0) # NaNs were added to separate two data
colnames(block_inc) <- colnames(block_dec) <- c("C1_t","C1_t_1","C1_t_2","R_t")

Do scenario exploration using simplex projection

In this demonstration, we use simplex projection to forecast the effects of changing resource (R).

# Normal multivariate simplex projection
scenario_res <- block_lnlp(block0,
                           method = "simplex",
                           columns = 1:4, # columns used to reconstruct the state space
                           target_column = 1,
                           stats_only = F,
                           silent = T)

# "Resource increased" situation
# Attractor is reconstructed using the original time series, but the dynamics is forecasted based on the resource-increased situation
scenario_inc <- block_lnlp(block_inc,
                           lib = c(1,nrow(block0)),
                           pred = c((nrow(block0)+2),nrow(block_inc)),
                           method = "simplex",
                           columns = 1:4, # columns used to reconstruct the state space.
                           target_column = 1,
                           stats_only = F,
                           silent = T)

## "Resource decreased" situation
# Attractor is reconstructed using the original time series, but the dynamics is forecasted based on the resource-decreased situation
scenario_dec <- block_lnlp(block_dec,
                           lib = c(1,nrow(block0)),
                           pred = c((nrow(block0)+2),nrow(block_dec)),
                           method = "simplex",
                           columns = 1:4, # columns used to reconstruct the state space.
                           target_column = 1,
                           stats_only = F,
                           silent = T)

# Extract model outputs
pred_nochange <- scenario_res[[1]]$model_output
pred_increase <- scenario_inc[[1]]$model_output[2000:nrow(scenario_inc[[1]]$model_output),]
pred_decrease <- scenario_dec[[1]]$model_output[2000:nrow(scenario_dec[[1]]$model_output),]


#### Visualize results
# Extract time window (for figure)
trange <- 1:50
plot(pred_nochange[trange,"obs"],
     type="l", col="black", lwd=2, xlab="time",
     ylab="Normalized C1 values", main="Scenario exploration")

# Add predicted values using noraml multivariate simplex projection
points(pred_nochange[trange,"pred"],bg="black",pch=21)
# Add predicted values under R increased situation
points(pred_increase[trange,"pred"],bg="red",pch=24)
# Add predicted values under R increased situation
points(pred_decrease[trange,"pred"],bg="blue",pch=25)
# Add legend
legend("topleft",c("predicted","R increased","R decreased"),
       pt.bg=c(1,2,4), pch=c(21,24,25), bty="n", cex = 0.8)

As can be seen in Figure 8 in the main text, increase in R does not always result in increased C1. This is because the effect of R on C1 depends additionally on the condition of other species, a phenomenon known as state-dependent behavior in dynamical systems.

References

Deyle ER, Fogarty M, Hsieh CH, Kaufman L, MacCall AD, Munch SB, Perretti CT, Ye H, Sugihara G (2013) Predicting climate effects on Pacific sardine. Proc Natl Acad Sci USA 110: 6430-6435. DOI: 10.1073/pnas.1215506110(PDF)

Deyle ER, May RM, Munch SB, Sugihara G (2016) Tracking and forecasting ecosystem interactions in real time. Proc R Soc Lond B Biol Sci 283: 20152258(PDF)