ecosyst.Rmd

---
title: "Threats to ecosystems"
subtitle: "Quantification of threats on the Norwegian Red List of ecosystems"
author: "Hanno Sandvik"
date: "24 May 2023"
output:
  md_document:
    toc: yes
---

This **R** code can be used to run the analyses of the Norwegian Red List for 
ecosystems and habitat types described in the paper "Metrics for quantifying 
how much different threats contribute to red lists of species and ecosystems" 
([Sandvik & Pedersen 2023](https://doi.org/10.1111/cobi.14105)).



## Variables

The following variables can be used to adjust the output.

**(1) Name of the data file.** The default downloads the Norwegian Red List data 
for ecosystems from [doi:10.5281/zenodo.7893216](https://doi.org/10.5281/zenodo.7893216). 
To analyse other Red Lists, use `url = ""` and provide the file name of the 
dataset as `file` (including file path, if needed).
```{r}
url  <- "https://zenodo.org/record/7893216/files/ecosyst.csv"
file <- "ecosyst.csv"
```

**(2) Handling of DD systems.** 
Decides whether Data Deficient ecosystems are excluded (if FALSE)
or randomly assigned to other Red List Categories (if TRUE).
```{r}
includeDD <- TRUE
```

**(3) Handling of unknown threats.**
Decides whether (if TRUE) or not (if FALSE) unknown threat factors
should be inferred from the distribution of the known threat factors.
```{r}
inferThreats <- FALSE
```

**(4) Disaggregation.**
Decides whether ecosystem major types should be disaggregated into their minor types.
```{r}
disaggregate <- FALSE
```

**(5) Weighting schemes.** 
Defines the weighting scheme for the Red List Index and the Expected Loss of Species.
(Defaults to "equal-steps" for RLI and the thresholds of the IUCN Red List 
Criterion E for ELS; other options are the IUCN Red List Criteria
"A1", "A2", "B1", "B2", "C1", "C2", "D1" and "D2" as well as "Ev2", "Ev3".)
```{r}
weightingRLI <- "equal-steps"
weightingELS <- "E"
```

**(6) Column names.** 
Column names in the dataset which contain Red List Categories, threat factors,
reasons for category change, and generation time, respectively. 
The three former ones need to be followed by the year of assessment
(for change, the year of the _second_ of the two relevant assessments).
So if the column name containing Red List Categories is _not_ named
something like "Categ18" or "Categ2018", this needs to be adjusted here!
```{r}
Categ  <- "Categ"
Threat <- "Threat"
Change <- "Change"
GTime  <- "GenTime"
```
Note the following formatting requirements of these columns:

* Red List Categories in the "Categ" column must match with the constants specified (see next section).
* Threat columns must contain text strings specifying threat. Each threat must be described as a sequence of (abbreviations for) (i) threat factor, (ii) timing, (iii) scope and (iv) severity, which are separated by _colons_; different threats to the same ecosystem are separated by _commas_.
* Change columns are needed only if the dataset contains results from more than one Red List. It must contain no more than one reason for change in Red List Category per ecosystem.
* Generation time is not applicable to ecosystems. However, `GTime` needs to be specified for the functions to work correctly.

**(7) Abbreviations used.**
What are the abbreviations used for unknown threats and for real status change? 
The four "unknowns" can occur in the `Threat` column(s), see previous item. 
(Defaults to the abbreviations used in the dataset analysed in the paper.
May need to be adjusted for other datasets.) The latter is only needed if 
Red List Categories are to be "back-cast" to earlier Red List assessments. 
It must occur in the `Change` column(s), see previous item.
```{r}
unknownThreat   <- "unknownf"
unknownTiming   <- "unknownt"
unknownScope    <- "unknownp"
unknownSeverity <- "unknownd"
realChange      <- "realchng"
```

**(8) Timings to include.** What is (are) the abbreviation(s) of the timing 
categories that should be considered (defaults to "ongoing").
```{r}
inclTiming <- "ongoingt"
```
If _all_ threats are to be included, irrespective of timing, this would need to be 
replaced (in terms of the abbreviations used in this dataset) by
`inclTiming <- c("onlypast", "suspendd", "ongoing", "onlyfutu", "unknownt")`.

**(9) Calculation of threat scores.** 
Decides whether threat scores are based on the product of scope and severity 
(if TRUE) or on severity alone (if FALSE).
```{r}
useScope <- FALSE
```
IUCN now [states](https://www.iucnredlist.org/resources/threat-classification-scheme) 
that severities should describe the decline _within the scope_ 
of a particular threat. If this definition has been followed, `useScope` 
should be `TRUE`. Previously, however, the definition of severity was ambiguous. 
In the Norwegian Red List analysed here, severity was used to describe the decline 
of the _entire area_. The default value (FALSE) assumes the latter situation, 
in which severity should _not_ be multiplied with score.
(See [here](scopesev.md) for some more detail.)

**(10) Number of simulations.** 
NB: the default takes a while. 
For exploration purposes, `nsim <- 1000` will suffice.
```{r}
nsim <- 100000
```

**(11) Re-create published estimates?** 
Decides whether the estimation of confidence intervals
should be re-created exactly as published (if TRUE) or be based on 
novel random numbers (if FALSE).
```{r}
re.create <- TRUE
```

**(12) File names of figures.** If you want to display the figures on screen,
keep the default. If you want to create PNG files, specify the file names
(including paths).
```{r}
fig4  <- ""
figS4 <- ""
```



## Constants

Constants should not normally need to be changed.
Changing them entails modifying some underlying assumptions.

**(1) Red List Categories** and their weights, extinction probabilities etc.
This data frame needs to contain all Red List Categories used in the 
Red List analysed of ecosystems that have been evaluated:
```{r}
RLcateg <- data.frame(
  name  = c(  "LC",   "NT",   "VU",   "EN",   "CR",   "CO"),
  LC    = c(  TRUE,  FALSE,  FALSE,  FALSE,  FALSE,  FALSE),
  EX    = c( FALSE,  FALSE,  FALSE,  FALSE,  FALSE,   TRUE),
  wt    = c(     0,      1,      2,      3,      4,      5),
  lowP  = c(  0.00,   0.05,   0.10,   0.20,   0.50,   1.00),
  uppP  = c(  0.00,   0.10,   0.20,   0.50,   1.00,   1.00),
  lowT  = c(   100,    100,    100,     50,     50,     50),
  uppT  = c(   100,    100,     50,     50,     50,     50),
  lowG  = c(     0,      0,      0,      0,      0,      0),
  uppG  = c(     0,      0,      0,      0,      0,      0),
  lowA1 = c(  0.00,   0.20,   0.30,   0.50,   0.80,   1.00),
  uppA1 = c(  0.00,   0.30,   0.50,   0.80,   1.00,   1.00),
  lowA2 = c(  0.00,   0.20,   0.30,   0.50,   0.80,   1.00),
  uppA2 = c(  0.00,   0.30,   0.50,   0.80,   1.00,   1.00),
  lowB1 = c( 55000,  55000,  50000,  20000,   2000,      0),
  uppB1 = c( 55000,  50000,  20000,   2000,      0,      0),
  lowB2 = c(  5500,   5500,   5000,   2000,    200,      0),
  uppB2 = c(  5500,   5000,   2000,    200,      0,      0),
  lowC  = c(  0.00,   0.15,   0.24,   0.40,   0.64,   1.00),
  uppC  = c(  0.00,   0.24,   0.40,   0.64,   1.00,   1.00),
  lowD  = c(  0.00,   0.15,   0.24,   0.40,   0.64,   1.00),
  uppD  = c(  0.00,   0.24,   0.40,   0.64,   1.00,   1.00),
  distr = c("unif", "unif", "unif", "unif", "decr", "unif"),
  beta  = c(    NA,     NA,     NA,     NA,     NA,     NA),
  stringsAsFactors = FALSE
)
```
The values in the dataframe are based on [IUCN (2016)](https://doi.org/10.2305/IUCN.CH.2016.RLE.2.en), 
[Bland et al. (2017)](https://doi.org/10.2305/IUCN.CH.2016.RLE.3.en) and the 
Norwegian guidance document ([Artsdatabanken 2018](https://artsdatabanken.no/Files/27210/)). 
The columns have the following meanings:

* The column "LC" identifies the Red List Category "Least Concern" (defaults to IUCN's abbreviation).
* The column "EX" identified the Red List Categories for ecosystem collapse (defaults to IUCN's abbreviations).
* The column "wt" provides the Red List Weight of the Category (defaults to equal-steps weighting).
* The columns "lowP" and "uppP" provide the lower and upper threshold values for extinction probability according to IUCN Red List Criterion E.
* The columns "lowT" and "uppT" provide the lower and upper threshold values for extinction time frames in _years_ according to IUCN Red List Criterion E.
* The columns "lowG" and "uppG" are not used for ecosystems but should contain zeros.
* The columns "lowA1" and "uppA1" provide the lower and upper threshold values for reduction in geographic distribution according to IUCN Red List Criterion A1.
* The columns "lowA2" and "uppA2" provide the lower and upper threshold values for reduction in geographic distribution according to IUCN Red List Criterion A2. (Note that A3 is not implemented.)
* The columns "lowB1" and "uppB1" provide the lower and upper threshold values for extents of occurrence (EOO) according to IUCN Red List Criterion B1.
* The columns "lowB2" and "uppB2" provide the lower and upper threshold values for areas of occupancy (AOO) according to IUCN Red List Criterion B2.
* The columns "lowC" and "uppC" provide the lower and upper threshold values for environmental degradation according to IUCN Red List Criteria C1 and C2 (estimated as the _product_ of _extent_ and _relative severity_ over a 50-year period; note that C3 is not implemented).
* The columns "lowD" and "uppD" provide the lower and upper threshold values for disruption of biotic processes or interactions according to IUCN Red List Criteria D1 and D2 (estimated as the _product_ of _extent_ and _relative severity_ over a 50-year period; note that D3 is not implemented).
* The column "distr" provides the distribution of extinction probabilities within the interval.
* The column "beta" is not currently needed (but may be needed if "distr" is changed).

**(2) Special Red List Categories.** What are the abbreviations used for 
data-deficient ecosystems and for ecosustems that have _not_ been evaluated? 
(Defaults to IUCN's abbreviations.)
```{r}
DD <- "DD"
notEval <- c("NA", "NE")
```

**(3) Downlisting.** What is added to a Red List Category to indicate downlisting? 
(Defaults to the degree symbol.)
If a Red List Category if followed by this symbol, it is assumed to have been 
_downlisted_ by _one_ Red List Category.
(This is not relevant for ecosystems but included for compatibility.)
```{r}
downlistSymbol <- "\u00b0"  # degree symbol in unicode
downlistSymbol <- iconv(downlistSymbol, Encoding(downlistSymbol), "latin1")
```

**(4) Scopes** and their threshold values. 
This data frame needs to contain all scope classes of threats used in the 
Red List analysed. Two versions of the data frame are provided, one for analysis 
of the Norwegian Red List data (the default) and one with the scope classes defined 
by IUCN ([2023](https://www.iucnredlist.org/resources/threat-classification-scheme)). 
The default is to use the Norwegian scope classes.
```{r}
ScopeNorway <- data.frame(
  name  = c("neglarea", "minority", "majority", "wholarea", "unknownp"),
  lower = c(      0.00,       0.05,       0.50,       0.90,       0.00),
  upper = c(      0.05,       0.50,       0.90,       1.00,       1.00),
  distr = c(    "unif",     "unif",     "unif",     "unif",     "beta"),
  beta  = c(        NA,         NA,         NA,         NA,          2),
  stringsAsFactors = FALSE
)

ScopeIUCN <- data.frame(
  name  = c("minority", "majority", "wholarea", "unknownp"),
  lower = c(      0.00,       0.50,       0.90,       0.00),
  upper = c(      0.50,       0.90,       1.00,       1.00),
  distr = c(    "unif",     "unif",     "unif",     "beta"),
  beta  = c(        NA,         NA,         NA,          2),
  stringsAsFactors = FALSE
)

Scope <- ScopeNorway
```
Values are the proportions of the total area affected by a threat 
([Artsdatabanken 2018](https://artsdatabanken.no/Files/27210/); cf. 
[IUCN 2023](https://www.iucnredlist.org/resources/threat-classification-scheme)).
The columns have the following meanings:

* The column "name" contains the abbreviations used for the scope classes.
* The column "lower" contains the lower limit of the respective interval.
* The column "upper" contains the upper limit of the respective interval.
* The column "distr" contains the distribution of values within the respective
interval (possible values: "unif", "incr", "decr", "beta").
* The column "beta" contains the beta parameter of a Beta distribution
(a numeric values if `distr == "beta"`, and `NA` otherwise).

**(5) Severities** and their threshold values.
This data frame needs to contain all severity classes of threats used in the Red 
List analysed. Two versions of the data frame are provided, one for analysis of the 
Norwegian Red List data (the default) and one with the severity classes defined 
by IUCN ([2023](https://www.iucnredlist.org/resources/threat-classification-scheme)). 
The default is to use the Norwegian severity classes.
```{r}
SeverityNorway <- data.frame(
  name  = c("negldecl", "slowdecl", "rapidecl", "unknownd"),
  lower = c(      0.00,       0.02,       0.20,       0.00),
  upper = c(      0.02,       0.20,       1.00,       1.00),
  distr = c(    "incr",     "unif",     "decr",     "beta"),
  beta  = c(        NA,         NA,         NA,         20),
  stringsAsFactors = FALSE
)

SeverityIUCN <- data.frame(
  name  = c("nodeclin", "negldecl", "slowdecl", "rapidecl", "veryrapd", "fluctuat", "unknownd"),
  lower = c(      0.00,       0.00,       0.02,       0.20,       0.30,       0.02,       0.00),
  upper = c(      0.00,       0.02,       0.20,       0.30,       1.00,       0.20,       1.00),
  distr = c(    "unif",     "unif",     "unif",     "unif",     "decr",     "unif",     "beta"),
  beta  = c(        NA,         NA,         NA,         NA,         NA,         NA,         20),
  stringsAsFactors = FALSE
)

Severity <- SeverityNorway
```
Values correspond to the declines in area of occupancy over 50 years caused by 
a threat ([Artsdatabanken 2020](https://artsdatabanken.no/Files/41216/); cf. 
[IUCN 2023](https://www.iucnredlist.org/resources/threat-classification-scheme)).
The columns have the following meanings:

* The column "name" contains the abbreviations used for the severity classes.
* The column "lower" contains the lower limit of the respective interval.
* The column "upper" contains the upper limit of the respective interval.
* The column "distr" contains the distribution of values within the respective 
interval (possible values: "unif", "incr", "decr", "beta";
see [here](scopesev.md) for the rationale).
* The column "beta" contains the beta parameter of a Beta distribution
(a numeric values if `distr == "beta"`, and `NA` otherwise).

**(6) Time frame** for the Expected Loss of Systems, in years
(defaults to 50 years).
```{r}
TimeFrame <- 50
```



## Preliminaries

Load the set of functions belonging to this repository:
```{r}
eval(parse(text = readLines("function.R")))
```

Define further required variables, 
based on the variables and constants specified above:
```{r}
LC      <- RLcateg$name[RLcateg$LC]  # abbreviation(s) for systems of Least Concern
extinct <- RLcateg$name[RLcateg$EX]  # abbreviation(s) for collapsed systems
LC.EX   <- RLcateg$name              # Red List Categories of evaluated systems
RedListCat <- c(LC.EX, DD, notEval)  # all Red List Categories
```



## Read and check the data

Read the dataset "Norwegian Red List for ecosystems and habitat types":
```{r}
{
  foundFile <- FALSE
  if (file.exists(file)) {
    foundFile <- TRUE
  } else {
    if (nchar(url)) {
      downl <- try(download.file(url, file))
      foundFile <- !inherits(downl, "try-error")
    }
  }
  if (foundFile) {
    RL <- read.csv2(file, as.is=TRUE, dec=".", na.strings="n/a", encoding="latin1")
  } else {
    cat("The datafile was not found.\n")
  }
}
```

Check whether the data are as expected:
```{r}
{
  usedCategories <- checkRL(RL)
  years <- identifyYears(RL)
  cat("\nYears included in this dataset:\n")
  print(years)
  threats <- identifyThreats(RL)
  cat("\nThreat factors reported in this dataset:\n")
  print(threats)
}
```
We can ignore these two warnings in this case.

Ensure that `RedListCat` and `LC.EX` only contain categories 
that are actually used:
```{r}
RedListCat <- RedListCat %A% usedCategories
LC.EX      <-      LC.EX %A% usedCategories
```



## Prepare the data frame for analyses

**(1) Reverse downlisting**:
```{r}
RL <- uplist(RL) 
```
We wouldn't have needed this step, since the dataset does not contain 
downlisted Red List Categories.

**(2) Back-cast** knowledge from the most recent Red List to earlier ones:
```{r}
RL <- backCast(RL)
```
That's correct, so we wouldn't have needed this step either.

**(3) Calculate collapse probabilites** for all ecosystems:
```{r}
RL <- calcLoss(RL)
```
This warning does not have to bother us when we are analysing Red Lists 
of ecosystems.

**(4) Add columns for all threat factors**:
```{r}
RL <- addThreats(RL)
```



## Summarise the data

Summarise the Red List:
```{r}
tab <- summariseRL(RL, exclude = NULL)
```



## Analysis of threat factors

Estimate ΔRLI:
```{r}
DRLI <- DeltaRLI(RL)
```

Estimate δRLI and ELS~50~:
```{r}
drli <- dRLI(RL)
print(drli)
```

Confidence intervals on RLI:
```{r, eval=FALSE, echo=TRUE}
print(confidenceRLI(RL, nsim, "Categ18"))
```
```{r, eval=TRUE, echo=FALSE}
# With the default `nsim`, the above line would take quite a while... 
# Instead, cached results are loaded:
load("cache.rdata")
print(confE)
```

Confidence intervals on ΔRLI, δRLI and ELS~50~:
```{r, eval=FALSE, echo=TRUE}
results <- simulateDRLI(RL, nsim)
```
```{r, eval=TRUE, echo=FALSE}
results <- simulateDRLI(RL, nsim, showProgress = FALSE)
```



## Figure 4

The following script recreates Figure 4.

Simplify the data by collapsing minor threats:
```{r}
drli. <- drli$dRLI[, 1]
ELS.  <- drli$ELS50[,1]
drli.["otherthr"] <- sum(drli.[c("alienspe", "huntgath", "natcatas",
                                 "nativesp", "otherthr", "unknownf")])
ELS. ["otherthr"] <- sum(ELS. [c("alienspe", "huntgath", "natcatas",
                                 "nativesp", "otherthr", "unknownf")])
drli. <-     drli.[-which(names(drli.) %in% c("alienspe", "huntgath", "natcatas",
                                              "nativesp", "unknownf"))]
ELS.  <-     ELS. [-which(names(ELS.)  %in% c("alienspe", "huntgath", "natcatas",
                                              "nativesp", "unknownf"))]
ELS.  <- ELS.[order(drli., decreasing=T)]
drli. <- drli.[order(drli., decreasing=T)]
```

Plot a graph for δRLI:
```{r}
{
  xl <- c(6, 30)
  yl <- c(0.79, 1.02)
  if (nchar(fig4)) {
    png(fig4, 1500, 1200, res=180)
    xl <- c(7, 24)
    yl <- c(0.795, 1.01)
  }
  par(mai=c(0.06, 0.96, 0.06, 0.06), family="sans")
  plot(0, 0, xlim=xl, ylim=yl, xaxs="i", yaxs="i", xaxt="n", yaxt="n",
       xlab="", ylab="Red List Index", bty="n", cex.axis=1.2, cex.lab=1.8)
  axis(2, seq(0.7, 1, 0.1), F, T, lwd=1.5, lend=1)
  axis(2, seq(0.7, 1, 0.05), T, F, cex.axis=1.2, lwd=1.5, lend=1)
  axis(2, seq(0.7, 1, 0.01), F, T, tcl=-0.25, lwd=1.5, lend=1)
  rect(11, 0, 14, 1, col=grey(0.96))
  for (i in 5:1) {
    rect(11, 1 - sum(drli.[i:1]), 14, 1, lwd=1.2, col=grey(0.96 - i * 0.12))
  }
  rect(11, 0, 14, 1, lwd=2.4, col=NA)
  x1 <- 7.4; x2 <- 10.6
  rli <- RLI(RL$Categ18)
  lines(c(x1, x2), rep(rli, 2), lwd=2.4)
  lines(c(x1, x2), rep(1, 2), lwd=2.4)
  polygon(x1 + c(-0.2,  0.1,  0.1), rli + c(0, 0.002, -0.002), col="black")
  polygon(x2 + c( 0.2, -0.1, -0.1), rli + c(0, 0.002, -0.002), col="black")
  polygon(x1 + c(-0.2,  0.1,  0.1),  1  + c(0, 0.002, -0.002), col="black")
  polygon(x2 + c( 0.2, -0.1, -0.1),  1  + c(0, 0.002, -0.002), col="black")
  text(mean(c(x1, x2)), 1, "reference value", pos=3, cex=1.2)
  text(mean(c(x1, x2)), 1, "1.0000", pos=1, cex=1.2)
  text(mean(c(x1, x2)), rli, "RLI 2018", pos=3, cex=1.2)
  text(mean(c(x1, x2)), rli, "0.8069", pos=1, cex=1.2)
  text(15, 1 - 0.5 * sum(        0) - 0.5 * sum(drli.[1:1]),
       expression(paste("Land-use change (", bold("8.3"), " ecosystems lost)")), 
       pos=4, cex=1.2)
  text(15, 1 - 0.5 * sum(drli.[1:1]) - 0.5 * sum(drli.[1:2]),
       expression(paste("Climate change (", bold("3.2"), " ecosystems lost)")), 
       pos=4, cex=1.2)
  text(15, 1 - 0.5 * sum(drli.[1:2]) - 0.5 * sum(drli.[1:3]),
       expression(paste("Other threats (", bold("2.3"), " ecosystems lost)")), 
       pos=4, cex=1.2)
  text(15, 1 - 0.5 * sum(drli.[1:3]) - 0.5 * sum(drli.[1:4]),
       expression(paste("Human disturbance (", bold("1.8"), " ecosystems lost)")), 
       pos=4, cex=1.2)
  text(15, 1 - 0.5 * sum(drli.[1:4]) - 0.5 * sum(drli.[1:5]),
       expression(paste("Pollution (", bold("1.4"), " ecosystems lost)")), 
       pos=4, cex=1.2)
  text(15, 1 - 0.5 * sum(drli.[1:5]) - 0.5 * 0.205,
       expression(paste(bold("212"), " ecosystems ", italic("not"), " lost")), 
       pos=4, cex=1.2)
  lines(14:15, rep(1 - 0.5 * sum(         0) - 0.5 * sum(drli.[1:1]), 2), lwd=1.2)
  lines(14:15, rep(1 - 0.5 * sum(drli.[1:1]) - 0.5 * sum(drli.[1:2]), 2), lwd=1.2)
  lines(14:15, rep(1 - 0.5 * sum(drli.[1:2]) - 0.5 * sum(drli.[1:3]), 2), lwd=1.2)
  lines(14:15, rep(1 - 0.5 * sum(drli.[1:3]) - 0.5 * sum(drli.[1:4]), 2), lwd=1.2)
  lines(14:15, rep(1 - 0.5 * sum(drli.[1:4]) - 0.5 * sum(drli.[1:5]), 2), lwd=1.2)
  lines(14:15, rep(1 - 0.5 * sum(drli.[1:5]) - 0.5 * sum(     0.205), 2), lwd=1.2)
  if (nchar(fig4)) {
    dev.off()
  }
}
```



## Disaggregation

The ecosystems assessed for the Norwegian Red List are at different levels 
of the underlying [EcoSyst](https://doi.org/10.1111/geb.13164) framework. 
The following analyses show the results obtained when major ecosystem types 
are disaggregated into their subordinated minor ecosystem types:
```{r}
{
  RL. <- disaggrMajorTypes(RL, minor = "MnTypes", type = "TypeCode", 
                           id = "FileID", categ = "Categ18.18")
  RL. <-    calcLoss(RL.)
  RL. <-  addThreats(RL.)
  cat("Summary table:\n")
  tab <- summariseRL(RL., exclude = NULL)
  drli. <- dRLI(RL.)
  cat("\ndRLI after disaggregation into minor types:\n")
  print(drli.$dRLI)
  cat("\nELS50 after disaggregation into minor types:\n")
  print(drli.$ELS50)
}
```



### Appendix S10

The following script recreates Appendix S10 (Figure S4).

Simplify the data by collapsing minor threats:
```{r}
ELS.  <- drli.$ELS50[,1]
drli. <- drli.$dRLI[, 1]
ELS. ["otherthr"] <- sum(ELS. [c("alienspe", "huntgath", "natcatas",
                                 "nativesp", "otherthr", "unknownf")])
drli.["otherthr"] <- sum(drli.[c("alienspe", "huntgath", "natcatas",
                                 "nativesp", "otherthr", "unknownf")])
ELS.  <-     ELS. [-which(names(ELS.)  %in% c("alienspe", "huntgath", "natcatas",
                                              "nativesp", "unknownf"))]
drli. <-     drli.[-which(names(drli.) %in% c("alienspe", "huntgath", "natcatas",
                                              "nativesp", "unknownf"))]
ELS.  <- ELS. [c(3, 1, 4, 2, 5)]
drli. <- drli.[c(3, 1, 4, 2, 5)]
```

Plot a graph for δRLI:
```{r}
{
  xl <- c(6.5, 26.5)
  yl <- c(0.792, 1.02)
  if (nchar(figS4)) {
    png("c:\\art\\threats\\figS4.png", 1500, 1200, res=180)
    xl <- c(7, 24)
    yl <- c(0.795, 1.01)
  }
  par(mai=c(0.06, 0.96, 0.06, 0.06), family="sans")
  plot(0, 0, xlim=xl, ylim=yl, 
       xaxs="i", yaxs="i", xaxt="n", yaxt="n",
       xlab="", ylab="Red List Index", bty="n", cex.axis=1.2, cex.lab=1.8)
  axis(2, seq(0.7, 1, 0.10), F, T, lwd=1.5, lend=1)
  axis(2, seq(0.7, 1, 0.05), T, F, cex.axis=1.2, lwd=1.5, lend=1)
  axis(2, seq(0.7, 1, 0.01), F, T, tcl=-0.25, lwd=1.5, lend=1)
  rect(11, 0, 14, 1, col=grey(0.96))
  for (i in 5:1) {
    rect(11, 1 - sum(drli.[i:1]), 14, 1, lwd=1.2, col=grey(0.96 - i * 0.12))
  }
  rect(11, 0, 14, 1, lwd=2.4, col=NA)
  x1 <- 7.4; x2 <- 10.6
  rli <- RLI(RL.$Categ18.18)
  lines(c(x1, x2), rep(rli, 2), lwd=2.4)
  lines(c(x1, x2), rep(1, 2), lwd=2.4)
  polygon(x1 + c(-0.2,  0.1,  0.1), rli + c(0, 0.002, -0.002), col="black")
  polygon(x2 + c( 0.2, -0.1, -0.1), rli + c(0, 0.002, -0.002), col="black")
  polygon(x1 + c(-0.2,  0.1,  0.1),  1  + c(0, 0.002, -0.002), col="black")
  polygon(x2 + c( 0.2, -0.1, -0.1),  1  + c(0, 0.002, -0.002), col="black")
  text(mean(c(x1, x2)), 1, "reference value", pos=3, cex=1.2)
  text(mean(c(x1, x2)), 1, "1.0000", pos=1, cex=1.2)
  text(mean(c(x1, x2)), rli, "RLI 2018", pos=3, cex=1.2)
  text(mean(c(x1, x2)), rli, "0.8812", pos=1, cex=1.2)
  text(15, 1 - 0.5 * sum(         0) - 0.5 * sum(drli.[1:1]),
       expression(paste("Land-use change (", bold("17"), " ecosystems lost)")),
       pos=4, cex=1.2)
  text(15, 1 - 0.5 * sum(drli.[1:1]) - 0.5 * sum(drli.[1:2]),
       expression(paste("Climate change (", bold("8"), " ecosystems lost)")), 
       pos=4, cex=1.2)
  text(15, 1 - 0.5 * sum(drli.[1:2]) - 0.5 * sum(drli.[1:3]),
       expression(paste("Other threats (", bold("4"), " ecosystems lost)")), 
       pos=4, cex=1.2)
  text(15, 1 - 0.5 * sum(drli.[1:3]) - 0.5 * sum(drli.[1:4]),
       expression(paste("Human disturbance (", bold("2"), " ecosystems lost)")), 
       pos=4, cex=1.2)
  text(15, 1 - 0.5 * sum(drli.[1:4]) - 0.5 * sum(drli.[1:5]),
       expression(paste("Pollution (", bold("3"), " ecosystems lost)")), 
       pos=4, cex=1.2)
  text(15, 1 - 0.5 * sum(drli.[1:5]) - 0.5 * 0.205,
       expression(paste(bold("774"), " ecosystems ", italic("not"), " lost")), 
       pos=4, cex=1.2)
  lines(14:15, rep(1 - 0.5 * sum(         0) - 0.5 * sum(drli.[1:1]), 2), lwd=1.2)
  lines(14:15, rep(1 - 0.5 * sum(drli.[1:1]) - 0.5 * sum(drli.[1:2]), 2), lwd=1.2)
  lines(14:15, rep(1 - 0.5 * sum(drli.[1:2]) - 0.5 * sum(drli.[1:3]), 2), lwd=1.2)
  lines(14:15, rep(1 - 0.5 * sum(drli.[1:3]) - 0.5 * sum(drli.[1:4]), 2), lwd=1.2)
  lines(14:15, rep(1 - 0.5 * sum(drli.[1:4]) - 0.5 * sum(drli.[1:5]), 2), lwd=1.2)
  lines(14:15, rep(1 - 0.5 * sum(drli.[1:5]) - 0.5 * sum(     0.205), 2), lwd=1.2)
  if (nchar(figS4)) {
    dev.off()
  }
}
```