Run cuspra on empirical data

1. Load package and dataset

Load the cusp and cuspra R packages. To make sure that you have the last version of the package, you can run devtools::install_github(repo = "rfrelat/cuspra", dependencies = TRUE, build_vignettes = TRUE).

library(cusp)
library(cuspra)

Let’s now load the dataset:

data("ecosystem_ns")
# or 
# data("fishtraits_med")
# data("cod_nea") 

If you use another dataset, make sure to change the name of the object in the cusp() model with the correct name of columns (e.g. PC1 should be replaced by Biomass if you are using the cod_nea dataset).

2. Fit the cusp model

fitNS <- cusp(y ~ PC1, alpha ~ Fishing,
             beta ~ Temperature,  data=ecosystem_ns)
# Make sure to check the model results
# with summary(), evalcusp(), or visually with plot()
evalcusp(fitNS)
#>    rsquared delta.aicc   percin   pval.state
#> 1 0.8847192   61.91278 41.37931 1.321463e-33

Make sure to check the model results with summary(fitNS), plot(fitNS) and/or evalcusp(fitNS). The output of evalcusp() shows: - rsquared: the pseudo R2 of the cusp model (should be >0.3)} - delta.aicc: the difference between AICc of the linear and cusp model (should be >0)} - percin: percentage of points within the cusp area (should be >10)} - pval.state: p-value of the state variable (should be <0.05)}

3. Estimate the resilience

raNS <- cuspra(fitNS)

The output of the cusp resilience assessment provide 4 columns: - cuspRA: the resilience assessment, with values between 0 (low resilience) and 1 (high resilience) - dalpha: the horizontal component as distance to the cusp area (if negative, the point is inside the cusp area) - dbeta: the vertical component as distance to linearity (if negative the system is discontinuous) - sumd: the weighted average of the horizontal and vertical components (by default the horizontal component has double weight)

4. Visualize the results

layout(matrix(c(2,1,1), ncol=1))
pal <- plotra(raNS$cuspRA, fitNS)
par(mar=c(2,4,1,1))
plot(ecosystem_ns$Year, ecosystem_ns$PC1, pch=16,
     ylab="State variable", type="b",
     col=as.character(cut(raNS$cuspRA, breaks = pal$bk, labels = pal$pal)))