Performing Early Warning Signal Assessments (2024)

set.seed(123) #to ensure reproducible datalibrary(EWSmethods)

1. The data

EWSmethods comes bundled with two data objects whichallow you to practice using the uniEWS() andmultiEWS() functions in both transitioning andnon-transitioning data before applying it to your own use case.

"simTransComms" contains three replicate datasets of asimulated five species community that has been driven to transition bythe introduction of an invasive species (following Dakos 2018). Thiswill be our multivariate dataset when using multiEWS()although we may also use each time series in isolation inuniEWS().

"CODrecovery" contains three replicate datasets of asimulated cod ( Gadus morhua ) population that transitions froman overfished to a recovered state following the relaxation of fishingpressure. This data was first published by Clements et al.2019. While univariate, "CODrecovery" provides extrainformation on the body size of cod individuals which will improvecomposite EWSs estimated by uniEWS().

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

We can visualise a community from each of these datasets using thecode 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 attime ~= 180 in "simTransComms$community1" andyear ~= 2050 in "CODrecovery$scenario2".EWSmethods helpfully provides this information in eachdataset under the inflection_pt column.

However, EWS assessments are only meaningful if performed on dataprior to a transition. As EWsmethods provides the timepoint of transition for both datasets, we can truncate our time seriesto 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 ofpast/present tipping points often unknown. If past transitions are knownto have occurred, it may be prudent to follow the suggestions of O’Brien& Clements (2021) who show that the occurrence of a historictransition can mask an oncoming event and that only using data post thehistoric transition improves EWS reliability.

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

2. Univariate EWSs

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

Both computational approaches are able to calculate the same EWSindicators. A brief outline of each can be found in the following tableas well as their reference code in uniEWS() for themetrics = 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 focalsystem to provide a community level assessment of transition risk.

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

Using multiEWS() we can estimate each of thesemultivariate indicators in the same way as uniEWS() -specifying the method =, winsize = and/orburn_in = - but must provide a n x m dataframe/matrix ofall representative time series. The first column must again be anequally 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 andtherefore are ‘warnings’.

We could also use expanding windows to achieve a similar result.Note - no composite metric is computed inmultiEWS() as it is currently unknown how combiningmultivariate 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 circ*mstance, many of the indicators are warning at differenttimes (e.g."eigenMAF" at time ~= 65 or"meanAR" at time ~= 100) but that the vastmajority are warning in the last 20 time points. This highlights theusefulness of expanding windows over rolling as the exact time point ofwarning can be determined, and supports Weinans et al.’s (2021)suggestion that there is no superior multivariate EWS indicator; thebest 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 toimprove the reliability of early warning signal indicators (Dakos etal. 2012, Gama Dessavre et al. 2019, Lenton etal. 2012). EWSmethods therefore provides a simpledetrending function detrend_ts() which provides fourmethods 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 tounivariate EWS analysis. However for non-average multivariateindicators, we suggest assessments should be made on non-detrended dataas the trend is informative. Mutual information in particular suffersfollowing detrending.

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

Early warning signal indicators are particularly sensitive tocyclical data as the repeated non-linearity throughout the year/cycleperiod will be interpreted as a transition initially, before maskingfuture cycles and tipping points. It is therefore sensible to deseasonseasonal data prior to assessment.

EWSmethods provides a suite of average and time seriesdecomposition deseasoning techniques via the deseason_ts()function. This function takes a n x m dataframe of time series to bedeaseasoned. The first column must be a vector of dates with theincrement = and order = arguments indicatingthe data resolution (year, month, day) and the order of the date vector(ymd/dmy/ydm). The form of deseasoning can then be selected using themethod = 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 pronouncedcycles 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 stepsmatplot(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 todeseason using LOESS (locally weighted smoothing) by estimating thecyclical component for each time series and then subtracting it.method = "decompose"/method = "x11" perform a similarprocess but use classical decomposition and x11 ARIMA modelling (Ladiray& Quenneville, 2001) respectively. method = "average"is the simplest method where the average increment value is estimatedfor each unique increment and that value subtracted from each data pointthat shares that increment key. deseason_ts() simplyprovides the default procedure for each of these methods andconsequently we visualise the results of deseason_ts()before using it in downstream analyses.

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