This tutorial introduces how to perform both univariate andmultivariate early warning signal (EWS) assessments usingEWSmethods
. It will give examples of rolling and expandingwindow approaches for univariate data, introduce trait-based compositeEWSs and then conclude with an example of multivariate EWSs.
Greater detail on each function can be found at the Referencepage.
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.
simTransComms | CODrecovery |
---|---|
191.5 | 2055 |
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.
Rolling (A) vs expanding (C) window concept diagram. Panels B and Dintroduce the circ*mstance when a ‘warning’ is signalled in bothapproaches
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.
EWS indicator | Description | uniEWS() metric code |
---|---|---|
Standard deviation | Increasing variance/standard deviation is observed approaching atransition | "SD" |
Coefficient of variation | Equivalent to SD as is simply SD at time t divided by the mean SD ofthe time series | "cv" |
Autocorrelation at lag1 | Autocorrelation (similarity between successive observations)increases approaching a transition. The value of this indicator can beestimated as either the autocorrelation coefficient estimated from afirst order autoregressive model or the estimated autocorrelationfunction at lag1 | "ar1" - autoregressive model,"acf" - autocorrelation function |
Skewness | At a transition, the distribution of values in the time series canbecome asymmetric | "skew" |
Kurtosis | Kurtosis represents the system reaching more extreme values in thepresence of a transition. Due to the increased presence of rare valuesin the time series, the tails of the observation distribution widen | "kurt" |
Return rate | The inverse of the first-order term of a fitted autoregressive AR(1)model. Return rate is the primary quantity impacted by CSD – return ratedecreases as a tipping point is approached | "rr" |
Density ratio | Spectral reddening (high variance at low frequencies) occurs neartransition. The density ratio quantifies the degree of reddening as theratio of the spectral density at low frequency to the spectral densityat high frequency | "dr" |
Rolling windows
The rolling window approach is the most commonly used form of EWScomputation due to the work of Dakos et al 2012 and theearlywarnings
package. uniEWS()
accepts amethod =
and a winsize =
argument which callsthe expanding window method and creates a rolling windowwinsize
% of the total time series length.
Lets use an example where we are interested in the autocorrelation,variance and skewness of one of the five species inpre_simTransComms
. First we supply a dataframe/matrix of nx 2 dimensions (first column is an equally time sequence and the secondis the abundance/biomass time series) and the EWS indicator metrics. Theremaining arguments specify the form of computation, window size andplotting characteristics.
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")
Note how all EWS indicators begin to trend upwards attime ~= 170
which results in the positive Kendall Taucorrelation coefficient indicative of an oncoming transition/tippingpoint.
Expanding windows
Let’s explore the alternative expanding window approach. All we needto change in uniEWS()
is the method =
argument, and replace winsize =
withburn_in =
. Instead of specifying the size of the rollingwindow, burn_in =
dictates the number of datapointsuniEWS()
is to use to ‘train’ the algorithm. This mitigatesthe high number of false-positive signals resulting from the short timeseries length and high variability when few data points are supplied atthe beginning of assessment (O’Brien & Clements, 2021).
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")
Similar to the rolling window approach, EWS indicators calculated byexpanding windows have exceeded the 2σ threshold for more than twoconsecutive time points and thus identified warnings fromtime ~= 170
. However we are more confident in thisconclusion as the composite metrics also display this warning(ar1 + SD + skew
). These composite metrics simple sum thestandardised individual indicator strengths together and are known toprovide a more reliable signal than lone indicators (Clements &Ozgul, 2016).
Trait information
The final contribution by uniEWS()
is the ability tointegrate multiple information sources in the assessment. For example,including body size estimates improves assessment reliability byreducing false positive rate whilst increasing the number of truepositives (Clements and Ozgul 2016, Baruah et al.2020).uniEWS()
consequently accepts a trait =
argument where an additional trait time series can be combined with theother abundance-based EWSs as a composite metric. This capability isonly available if method = "expanding"
andmetrics =
contains "trait"
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.
Multivariate EWS indicator | Description | multiEWS() metric code | Average or dimension reduction technique |
---|---|---|---|
Mean standard deviation | Average variance across all time series representing the system | "meanSD" | Average |
Max standard deviation | The variance of the time series with the highest variance of allassessed time series | "maxSD" | Average |
Mean autocorrelation at lag1 | Average autocorrelation across all time series representing thesystem | "meanAR" | Average |
Max autocorrelation at lag1 | The autocorrelation of the time series with the highestautocorrelation of all assessed time series | "maxAR" | Average |
Dominant MAF (maximum autocorrelation factor)eigenvalue | The minimum eigenvalue of the system following MAF dimensionreduction | "eigenMAF" | Dimension reduction |
MAF autocorrelation | The autocorrelation of the data projected on to the first MAF –i.e.the autocorrelation of the first MAF. | "mafAR" | Dimension reduction |
MAF standard deviation | The variance of the data projected on to the first MAF – i.e.thevariance of the first MAF | "mafSD" | Dimension reduction |
First PC (principal component) autocorrelation | The autocorrelation of the data projected on to the first PC –i.e.the autocorrelation of the first PC | "pcaAR" | Dimension reduction |
First PC standard deviation/ Explained variance | The variance of the data projected on to the first PC – i.e.thevariance of the first PC | "pcaSD" | Dimension reduction |
First PC standard deviation/ Explained variance | The variance of the data projected on to the first PC – i.e.thevariance of the first PC | "pcaSD" | Dimension reduction |
Dominant eigenvalue of the covariance matrix | The maximum eigenvalue of the covariance matrix between allrepresentative time series | "eigenCOV" | Neither |
Maximum covariance | The maximum value of the covariance matrix between allrepresentative time series. | "maxCOV" | Neither |
Mutual information | The ‘amount of information’ one time series contains onanother. | "mutINFO" | Neither |
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?
Rolling windows
The simplicity of the rolling window approach also limits itsusefulness. A ‘warning’ is indicated when an EWS displays a strongpositive Kendall Tau correlation with time. However, it is unclear whatconstitutes a ‘strong’ correlation in this context with publishedwarnings ranging from 0.5 through to 0.9 (Dakos et al. 2012,Dablander et al. 2022, Southall et al. 2022). Thestrength of correlation therefore appears to be context dependent andsystem specific.
An alternative approach suggested by Dakos et al. (2012) isto generate random permutations of the assessed time series and thencompare the estimated Kendall Tau coefficients to that of the ‘true’time series. If the true coefficient is stronger than 95% of thepermuted coefficients then that represents a warning.
Expanding windows
Expanding windows have a stronger evidence base for what constitutesa ‘warning’. Clements et al. 2019 show that two consecutivetransgressions of the 2σ threshold reduce false-postive rates andimprove the number of true-postives. This has been validated by Southallet al. 2022 although they suggest that more than twoconsecutive signals should be aimed for.
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:
Detrending method code | Description |
---|---|
linear | Linear detrending - returns the residuals of a fitted linear modelbetween the time index and time series |
loess | Local polynomial regression smoothing - subtracts a smooth curvefitted by local polynomial regression from the observed time series |
gaussian | Gaussian kernel smoothing - subtracts a smooth curve estimated by akernel based weighted moving average from the observed time series |
first.difference | First differencing - subtracts the lagged time series from theunlagged 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.
Clements, C.F. & Ozgul, A. (2016) Including trait-based earlywarning signals helps predict population collapse. NatureCommunications, 7, 10984. doi:10.1038/ncomms10984
Clements, C.F., McCarthy, M.A.& Blanchard, J.L. (2019) Earlywarning signals of recovery in complex systems. NatureCommunications, 10, 1681. doi:10.1038/s41467-019-09684-y
Cleveland, R.B., Cleveland, W.S., McRae, J.E., & Terpenning, I.J.(1990) STL: A seasonal-trend decomposition procedure based on loess.Journal of Official Statistics, 6, 1, 3–33. http://bit.ly/stl1990
Dablander, F., Heesterbeek H., Borsboom D. & Drake J.M. (2022)Overlapping timescales obscure early warning signals of the secondCOVID-19 wave. Proceedings of the Royal Society B , 289,20211809. doi:10.1098/rspb.2021.1809
Dakos, V. (2018) Identifying best-indicator species for abrupttransitions in multispecies communities. Ecological Indicators,94, 494–502. doi:10.1016/j.ecolind.2017.10.024
Dakos, V., Carpenter, S.R., Brock, W.A., Ellison, A.M., Guttal, V.,Ives, A.R., et al.(2012) Methods for detecting early warnings ofcritical transitions in time series illustrated using simulatedecological data. PLoS One, 7, e41010. doi:10.1371/journal.pone.0041010
Gama Dessavre, A., Southall, E., Tildesley, M.J. & Dyson, L.(2019) The problem of detrending when analysing potential indicators ofdisease elimination, Journal of Theoretical Biology, 481,183-193. doi:10.1016/j.jtbi.2019.04.011
Ladiray, D., & Quenneville, B. (2001) Seasonal Adjustment withthe X-11 Method. Springer. doi:10.1007/978-1-4613-0175-2
Lenton T M., Livina V.N., Dakos V., van Nes E.H. & Scheffer M.(2012) Early warning of climate tipping points from critical slowingdown: comparing methods to improve robustness PhilosophicalTransactions of the Royal Society A 370, 1185–1204. doi:10.1098/rsta.2011.03047
O’Brien, D.A. & Clements, C.F. (2021) Early warning signalreliability varies with COVID-19 waves. Biology Letters, 17,20210487. doi:10.1098/rsbl.2021.0487
Persons, W.M. (1919) General considerations and assumptions. TheReview of Economics and Statistics, 1, 1,, 5–107. doi:10.2307/1928754
Southall E., Tildesley, M.J. & Dyson, L. (2022) How early can anupcoming critical transition be detected? medRxiv2022.05.27.22275693. doi:10.1101/2022.05.27.22275693
Weinans, E., Quax, R., van Nes, E.H. & van de Leemput, I.A.(2021) Evaluating the performance of multivariate indicators ofresilience loss. Scientific Reports, 11, 9148. doi:10.1038/s41598-021-87839-y