Introduction:

Hsieh et al. 2008 propose a method to detect non-linearity in ecological time series using the simplex and s-map projections developed by Sugihara and co-workers. Generally, non-linear time series analysis requires long-time series to detect signals. The approach by Hsieh et al. combines individual time series of ecologically similar species into one long time series to which simplex and s-maps can be applied (termed Dewdrop regression by the authors). To illustrate their approach, a five species competition model was used to generate time series and the ability of the method to get back the original parameter values of the simulation tested under various levels of noise and lengths of the time series.

The goals for this reproduction are: * set up a simulation model of multi-species competition * understanding simplex and s-maps as forecasting tools using the recently released rEDM package. * understand the approach taken to combine several short time series into one long time series to detect non-linearity and make useful forecasts

Step 1: Simulate 5 species dynamics (as described in detail in Appendix B):

First we reproduced the dynamics of the 5 species community as follows:

The following R code implements the simulation model:

# This code implements a multispecies competition model based on the logistic map (https://en.wikipedia.org/wiki/Logistic_map)
# the equilibrium density depends on the growth rate (steady state at values between 1 and 2, extinction below 1, chaos > 3.7, fluctuations > 3)

# set seed
set.seed(190981)

## Specify the number of species
S <- 5

#growth rate in the eight-mode oscillation regime below the threshold to chaos
r <- 3.57

# specify alphas
a <- 0.05

# simulate 1000 time steps
t=seq(1,1000, by=1)

# specify matrix to hold results
N <- matrix(nrow=S, ncol=length(t))

# set initial conditions randomly
N[,1] <- runif(S)

# simulate dynamics using the discrete logistic growth model; add process noise at the end
for (i in 1:(length(t)-1)){
  
  for (j in 1:5){

    N[j,i+1] = (r * N[j,i] * (1- N[j,i]) - a*sum(N[-j,i]) ) * runif(1,0.85,1.15)
    
  }
}

# add some observation noise to time series
N <- N * runif(length(N), 0.95, 1.05)

output <- as.data.frame(t(N[1:S,]))
names(output) <- paste0("species ", 1:S)

# check for correlation in abundances
#plot(output$`species 1`, output$`species 5`)

#write.csv(output, "/Users/Frank/Documents/Git projects/RREEBES/Hsieh_et_al_2008/data/5sp_comp_sim_data.csv", row.names=F)

five_sp_comp <- read.csv("/Users/Frank/Documents/Git projects/RREEBES/Hsieh_et_al_2008/data/5sp_comp_sim_data.csv")
#names(five_sp_comp) <- c("time", paste0("species.", 1:5))
five_sp_comp$time <- seq(1,1000)

five_sp_comp <- five_sp_comp[,c(6,1:5)]

Plot shows time series for the 5 species competition model.

# plot results
matplot(t,t(N[1:S,]), "l",ylim=c(0,1), xlim=c(0,100),col=2:5)