1. The data

EWSmethods comes bundled with two data objects which allow you to practice using the uniEWS() and multiEWS() functions in both transitioning and non-transitioning data before applying it to your own use case.

"simTransComms" contains three replicate datasets of a simulated five species community that has been driven to transition by the introduction of an invasive species (following Dakos 2018). This will be our multivariate dataset when using multiEWS() although we may also use each time series in isolation in uniEWS().

"CODrecovery" contains three replicate datasets of a simulated cod ( Gadus morhua ) population that transitions from an overfished to a recovered state following the relaxation of fishing pressure. This data was first published by Clements et al. 2019. While univariate, "CODrecovery" provides extra information on the body size of cod individuals which will improve composite EWSs estimated by uniEWS().

#Load the two datasets in to the session
data("simTransComms")

data("CODrecovery")

We can visualise a community from each of these datasets using the code below:

matplot(simTransComms$community1[,3:7], type = "l", xlab = "Time", ylab = "Density",  main = "Transitioning five species community")

plot(x = CODrecovery$scenario2$time, y = CODrecovery$scenario2$biomass, type = "l", xlab = "Year", ylab = "Abundance", main = "Recovering cod population")

These plots show that a transition takes place at time ~= 180 in "simTransComms$community1" and year ~= 2050 in "CODrecovery$scenario2". EWSmethods helpfully provides this information in each dataset under the inflection_pt column.

However, EWS assessments are only meaningful if performed on data prior to a transition. As EWsmethods provides the time point of transition for both datasets, we can truncate our time series to pre-transition data only.

pre_simTransComms <- subset(simTransComms$community1,time < inflection_pt)

pre_CODrecovery <- subset(CODrecovery$scenario2,time < inflection_pt)

In reality, EWSs will be assessed in real-time with the presence of past/present tipping points often unknown. If past transitions are known to have occurred, it may be prudent to follow the suggestions of O’Brien & Clements (2021) who show that the occurrence of a historic transition can mask an oncoming event and that only using data post the historic transition improves EWS reliability.

Now the data has been loaded and truncated, it can now be passed to uniEWS() and multiEWS() to perform EWS assessments.

2. Univariate EWSs

EWSmethods provides two computational approaches to calculate univariate EWSs via the uniEWS() function - rolling vs expanding windows. The difference between the two is evident in the figure below but simply, rolling windows estimate EWSs in subsets of the overall time series before ‘rolling’ on one data point and reassessing, Conversely, expanding windows add data points sequentially in to an ‘expanding’ assessment and then standardises against the running mean and standard deviation of the previous window.

Rolling (A) vs expanding (C) window concept diagram. Panels B and D introduce the circumstance when a ‘warning’ is signalled in both approaches

Both computational approaches are able to calculate the same EWS indicators. A brief outline of each can be found in the following table as well as their reference code in uniEWS() for the metrics = argument.

rolling_ews_eg <- uniEWS(data = pre_simTransComms[,c(2,5)],
                         metrics = c("ar1","SD","skew"),
                         method = "rolling",
                         winsize = 50)

plot(rolling_ews_eg,  y_lab = "Density")

expanding_ews_eg <- uniEWS(data = pre_simTransComms[,c(2,5)],
                         metrics = c("ar1","SD","skew"),
                         method = "expanding",
                         burn_in = 50,
                         threshold = 2)

plot(expanding_ews_eg, y_lab = "Density")

trait_ews_eg <- uniEWS(data = pre_CODrecovery[,c(2,3)],
                         metrics = c("ar1","SD","trait"), #note "trait" is provided here
                         method = "expanding",
                         trait = pre_CODrecovery$mean.size, #and here
                         burn_in = 15, #small burn_in due to shorter time series
                         threshold = 2)

plot(trait_ews_eg, y_lab = "Density", trait_lab = "Mean size (g)")

3. Multivariate EWSs

A more powerful and informative form of EWS are multivariate EWSs. These indicators combine multiple time series/measurements of the focal system to provide a community level assessment of transition risk.

There are two primary forms of multivariate EWS, those which are averaged univariate EWS across all time series and those which are assessments made on a dimension reduction of all representative time series. A brief outline of each can be found in the following table as well as their reference code in multEWS() for the metrics = argument. See Weinans et al. (2021) for a rigorous testing of these multivariate EWSs in a simulated system.

Using multiEWS() we can estimate each of these multivariate indicators in the same way as uniEWS() - specifying the method =, winsize = and/or burn_in = - but must provide a n x m dataframe/matrix of all representative time series. The first column must again be an equally spaced time sequence.

A rolling window assessment would therefore be coded as such:

multi_ews_eg <- multiEWS(data = pre_simTransComms[,2:7],
                         metrics = c("meanAR","maxAR","meanSD","maxSD","eigenMAF","mafAR","mafSD","pcaAR","pcaSD","eigenCOV","maxCOV","mutINFO"),
                         method = "rolling",
                         winsize = 50)

plot(multi_ews_eg)

Many of these indicators are postively correlated with time and therefore are ‘warnings’.

We could also use expanding windows to achieve a similar result. Note - no composite metric is computed in multiEWS() as it is currently unknown how combining multivariate EWS indicators influences prediction reliability.

multi_ews_eg2 <- multiEWS(data = pre_simTransComms[,2:7],
                         method = "expanding",
                         burn_in = 50,
                         threshold = 2)

plot(multi_ews_eg2)

In this circumstance, many of the indicators are warning at different times (e.g. "eigenMAF" at time ~= 65 or "meanAR" at time ~= 100) but that the vast majority are warning in the last 20 time points. This highlights the usefulness of expanding windows over rolling as the exact time point of warning can be determined, and supports Weinans et al.’s (2021) suggestion that there is no superior multivariate EWS indicator; the best fit depends on the scenario the system is subject t0.

1. How do I interpret EWSs?

2. Do I need to detrend my data?

The is a body of work that suggests detrending may be necessary to improve the reliability of early warning signal indicators (Dakos et al. 2012, Gama Dessavre et al. 2019, Lenton et al. 2012). EWSmethods therefore provides a simple detrending function detrend_ts() which provides four methods of detrending a time series:

Example:

detrend_dat <- detrend_ts(data =  pre_simTransComms[,2:7],  method = "loess", span = 0.75, degree = 2)

matplot(x = detrend_dat$time, y = detrend_dat[,2:6], type = "l", xlab = "Date", ylab = "Density",  main = "LOESS detrended five species community")

The current consensus is that detrending is required prior to univariate EWS analysis. However for non-average multivariate indicators, we suggest assessments should be made on non-detrended data as the trend is informative. Mutual information in particular suffers following detrending.

3. My data is seasonal/contains cycles. Can I still use Early Warning Signals?

Early warning signal indicators are particularly sensitive to cyclical data as the repeated non-linearity throughout the year/cycle period will be interpreted as a transition initially, before masking future cycles and tipping points. It is therefore sensible to deseason seasonal data prior to assessment.

EWSmethods provides a suite of average and time series decomposition deseasoning techniques via the deseason_ts() function. This function takes a n x m dataframe of time series to be deaseasoned. The first column must be a vector of dates with the increment = and order = arguments indicating the data resolution (year, month, day) and the order of the date vector (ymd/dmy/ydm). The form of deseasoning can then be selected using the method = argument.

Lets create some dummy monthly data that we can deseason.

spp_data <- matrix(nrow = 5*12, ncol = 5)

seasonal_cycle <-  20*sin(2*pi*(1:5*12)/12)

spp_data <- sapply(1:dim(spp_data)[2], function(x){

  spp_data[,x] <- ts(rnorm(5*12,mean = 20, sd = 3) + seasonal_cycle, freq = 12, start = c(2000, 1)) #add seasonal cycle to random noise
  
  })

multi_spp_data <- cbind("time" = base::seq(base::as.Date('2000/01/01'), base::as.Date('2004/12/01'), by = "month"), as.data.frame(spp_data))
  
matplot(x = multi_spp_data$time, y = multi_spp_data[,2:6], type = "l", xlab = "Date", ylab = "Density",  main = "Seasonal five species community")

As multi_spp_data is random, there are few pronounced cycles but for the sake of this tutorial, deseason_ts() would be applied to it as such:

deseas_dat <- deseason_ts(data = multi_spp_data, increment = "month", method = "average", order = "ymd")
#> data successfully aggregated into monthly time steps

matplot(x = deseas_dat$date, y = deseas_dat[,2:6], type = "l", xlab = "Date", ylab = "Density",  main = "Deseasoned five species community")

The method = "stl" argument shows that we have chosen to deseason using LOESS (locally weighted smoothing) by estimating the cyclical component for each time series and then subtracting it. method = "decompose"/method = "x11" perform a similar process but use classical decomposition and x11 ARIMA modelling (Ladiray & Quenneville, 2001) respectively. method = "average" is the simplest method where the average increment value is estimated for each unique increment and that value subtracted from each data point that shares that increment key. deseason_ts() simply provides the default procedure for each of these methods and consequently we visualise the results of deseason_ts() before using it in downstream analyses.

In our example, we can see that large monthly values shared across multiple years have been shrunk - e.g. the black dashed species at ~2001 - whereas anomalous values have been maintained - e.g. the green dotted species at ~2002.75.