EWSmethods <a href="https://duncanobrien.github.io/EWSmethods/"><img src="man/figures/logo.svg" align="right" height="150"/></a>

EWSmethods is a user friendly interface to various methods of performing Early Warning Signal (EWS) assessments. This R package allows the user to input univariate or multivariate data and perform either traditional rolling window (e.g. Dakos et al. 2012) or expanding window (Drake and Griffen, 2010) EWS approaches. Publication standard and ggplot inspired figures can also be generated during this process. EWSmethods also provides an R interface to EWSNet, a deep learning modelling framework for predicting critical transitions (Deb et al. 2022).

This README is a quick-fire introduction to the package, but indepth tutorials are available for the following topics:

early warning signal assessments: tutorial
machine learning: tutorial
resilience indicators: tutorial

Installation

You can install the stable version of EWSmethods from CRAN with:

install.packages("EWSmethods")

Alternatively, you can install the development version from GitHub with:

install.packages("devtools")
devtools::install_github("duncanobrien/EWSmethods", dependencies = TRUE)

⚠️ Warning! - Due to the large file size of the EWSNet model weights (~220MB), EWSmethods does not come bundled with them. Users must download weights using ewsnet_reset() before interfacing with EWSNet.

Getting Started

library(EWSmethods)

The remainder of this page will introduce the two primary ways of interacting with EWSmethods for your critical transition forecasting needs. For specific function help, please refer to the Reference page.

1. Early Warning Signals

Early warning signals are a collection of summary statistics that attempt to characterise the phenomenon of critical slowing down (CSD). As a system approaches a tipping point (or bifurcation), it takes longer and longer for it to recover when it is pushed away from stability (Dakos et al. 2012). This increased return rate is a manifestation of CSD and can be detectable in data. EWSmethods provides a collection of these summary statistics which can be calculated from either univariate or multivariate time series using uniEWS() and multiEWS() respectively.

The univariate approach

Imagine we have 50 years of monitoring data for a local population of skylarks (Alauda arvensis) and have measured mean body mass data throughout this period as well.

We could calculate either rolling or expanding window EWSs using the abundance data via uniEWS() but decide to initially focus on rolling windows. We would therefore parameterise the function to do so as below:

set.seed(125) #seed to ensure reproducible results

skylark_data <- data.frame(time = seq(1:50), abundance = rnorm(50,mean = 100,sd=20), trait = rnorm(50,mean=40,sd=5)) #dummy skylark dataset

ews_metrics <- c("SD","ar1","skew") #the early warning signal metrics we wish to compute

roll_ews <- uniEWS(data = skylark_data[,1:2], metrics =  ews_metrics, method = "rolling", winsize = 50) #lets use a rolling window approach

roll_ews$EWS$cor #return the Kendall Tau correlations for each EWS metric
#>            SD        ar1       skew
#> tau 0.5446154 -0.3907692 -0.7046154

We can then use the resulting figures to identify oncoming transitions. In this case, we expect no transition as the data is randomly sampled from a normal distribution and this is evident in the Kendall Tau values, with no strong positive correlation with time:

plot(roll_ews,  y_lab = "Skylark abundance")

Alternatively, we may be more interested in expanding windows as that approach standardises the changing EWS metrics over time and therefore allows the strength of multiple signals to be combined. We could achieve this using the following code:

exp_ews <- uniEWS(data = skylark_data[,1:2], metrics =  ews_metrics, method = "expanding", burn_in = 10, threshold = 2,  tail.direction = "one.tailed") #lets use a rolling window approach

head(exp_ews$EWS) #return the head of the EWS dataframe
#>   time metric.score metric.code rolling.mean rolling.sd threshold.crossed
#> 1   10    0.0000000         ar1    0.0000000         NA                 0
#> 2   11   -0.7071068         ar1   -0.3535534  0.5000000                 0
#> 3   12   -0.9223669         ar1   -0.5431579  0.4825449                 0
#> 4   13   -0.1500572         ar1   -0.4448827  0.4403012                 0
#> 5   14    0.8428688         ar1   -0.1873324  0.6906950                 0
#> 6   15   -1.8668306         ar1   -0.4672488  0.9229121                 0
#>   count.used        str
#> 1  108.91249         NA
#> 2   93.56813 -0.7071068
#> 3  109.56960 -0.7858522
#> 4  103.92341  0.6695997
#> 5  114.29655  1.4915428
#> 6   80.79726 -1.5164845

And again, we can then use the resulting figures to identify oncoming transitions. Whilst we have some trangressions of the 2σ threshold, we only consider these signals “warnings” if two or more consecutive signals are identified (Clements et al. 2019).

plot(exp_ews, y_lab = "Skylark abundance")

A second benefit of the expanding window approach is that additional information can be used to improve the reliability of the assessment. Including trait information has been shown to decrease the likelihood of both false positive and false negative signals (Clements and Ozgul, 2016; Baruah et al. 2019) and therefore should be considered if possible.

For example, in our hypothetical skylark dataset, we have measured average population body mass. This data can then be delivered to the uniEWS() function in EWSmethods, using the trait argument.

trait_metrics <- c("SD", "ar1", "trait")

exp_ews_trait <- uniEWS(data = skylark_data[,1:2], metrics =  trait_metrics, trait = skylark_data$trait, method = "expanding", burn_in = 10, threshold = 2, tail.direction = "one.tailed")

plot(exp_ews_trait, y_lab = "Skylark abundance", trait_lab = "Body mass (g)", trait_scale = 5)

The multivariate approach

If we had data from multiple timeseries/measurements of the same system, we might be interested in multivariate early warning signals. These indicators exploit either dimension reduction techniques (such as Principal Component Analysis) or community average estimates to give an overall measure of system resilience (see Weinans et al. 2021 for an overview of each indicator).

Here we’ve constructed another hypothetical dataset representing five related populations of Caribbean reef octopus (Octopus briareus) in Bahamian salt water lakes (O’Brien et al. 2020) and are interested in assessing the resilience of this metapopulation. The following code shows how we would achieve this using the EWSmethods function multiEWS().

set.seed(123)

octopus_spp_data <- matrix(nrow = 50, ncol = 5)
octopus_spp_data <- as.data.frame(cbind("time"=seq(1:50),sapply(1:dim(octopus_spp_data)[2], function(x){octopus_spp_data[,x] <- rnorm(50,mean=500,sd=200)}))) #create our hypothetical, uncollapsing ecosystem

oct_exp_ews <- multiEWS(data = octopus_spp_data, method = "expanding", threshold = 2, tail.direction = "one.tailed") #lets use an expanding window approach

The figure again shows that one multivariate EWS indicator has expressed a warning, but that overall, no transition is anticipated.

plot(oct_exp_ews)

2. EWSNet

The other half of EWSmethods allows you to query the Python-based EWSNet via an easy to use R workflow. Here is a simple example that details how to first prepare your R session to communicate with Python (using the excellent reticulate R package) and then calls EWSNet to assess the probability of a transition occurring in the skylark time series. This is a two step process where we must a) call ewsnet_init() before b) using ewsnet_predict().

However, because this is the first time using EWSNet via EWSmethods, we must first download the pretrained model weights from https://ewsnet.github.io.

options(timeout = max(600, getOption("timeout"))) #due to possible internet issues, increase the timeout options from 60 seconds to 600

ewsnet_reset(remove_weights = FALSE, auto = TRUE)
#> Model weights downloaded

Now we can setup our R session and interface with EWSNet.

ewsnet_init(envname = "EWSNET_env", auto = TRUE) #prepares your workspace using 'reticulate' and asks to install Anaconda (if no appropriate Python found) and/or a Python environment before activating that environment with the necessary Python packages
#> Attention: may take up to 10 minutes to complete
#> EWSNET_env successfully found and activated. Necessary python packages installed

library(reticulate)

reticulate::py_config() #confirm that "EWSN

EWSmethods

Install / Use

README