species.Rmd

---
title: "Threats to species"
subtitle: "Quantification of threats on the Norwegian Red Lists of species"
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 Lists 
for species 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 species 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/species.csv"
file <- "species.csv"
```

**(2) Handling of DD species.** 
Decides whether Data Deficient species 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) 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", "C" and "D" as well as "Ev2", "Ev3".)
```{r}
weightingRLI <- "equal-steps"
weightingELS <- "E"
```

**(5) 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 "Categ21" or "Categ2021", 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 threats. 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 species 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 species.
* Generation time, measured in _years_, must be a numerical variable.

**(6) Abbreviations used for threats.**
What are the abbreviations used for unknown threats? 
They 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.)
```{r}
unknownThreat   <- "unknownf"
unknownTiming   <- "unknownt"
unknownScope    <- "unknownp"
unknownSeverity <- "unknownd"
```

**(7) Abbreviation used for real change.**
What is the abbreviations used for real population changes? 
This 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 item (5). 
(Defaults to the abbreviation used in the dataset analysed in the paper.
May need to be adjusted for other datasets.)
```{r}
realChange <- "realpopu"
```

**(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 population 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 Lists analysed here, severity was used to describe the decline
of the _entire population_. The default value (FALSE) assumes the latter situation, 
in which severity should _not_ be multiplied with scope.
(See [here](scopesev.md) for some more detail.)

**(10) Number of simulations.** 
NB: the default takes several hours!
For exploration purposes, `nsim <- 1000` will suffice.
For pure illustration, `nsim <- 100` is enough.
```{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}
fig1  <- ""
fig2  <- ""
fig3  <- ""
figS1 <- ""
figS2 <- ""
figS3 <- ""
```



## 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 species that have been evaluated:
```{r}
RLcateg <- data.frame(
  name  = c(  "LC",   "NT",   "VU",   "EN",   "CR",   "RE",   "EW",   "EX"),
  LC    = c(  TRUE,  FALSE,  FALSE,  FALSE,  FALSE,  FALSE,  FALSE,  FALSE),
  EX    = c( FALSE,  FALSE,  FALSE,  FALSE,  FALSE,   TRUE,   TRUE,   TRUE),
  wt    = c(     0,      1,      2,      3,      4,      5,      5,      5),
  lowP  = c(  0.00,   0.05,   0.10,   0.20,   0.50,   1.00,   1.00,   1.00),
  uppP  = c(  0.00,   0.10,   0.20,   0.50,   1.00,   1.00,   1.00,   1.00),
  lowT  = c(   100,    100,    100,     20,     10,     10,     10,     10),
  uppT  = c(   100,    100,     20,     10,     10,     10,     10,     10),
  lowG  = c(     0,      0,      0,      5,      3,      1,      1,      1),
  uppG  = c(     0,      0,      5,      3,      1,      1,      1,      1),
  lowA1 = c(  0.00,   0.25,   0.50,   0.70,   0.90,   1.00,   1.00,   1.00),
  uppA1 = c(  0.00,   0.50,   0.70,   0.90,   1.00,   1.00,   1.00,   1.00),
  lowA2 = c(  0.00,   0.15,   0.30,   0.50,   0.80,   1.00,   1.00,   1.00),
  uppA2 = c(  0.00,   0.30,   0.50,   0.80,   1.00,   1.00,   1.00,   1.00),
  lowB1 = c( 40000,  40000,  20000,   5000,    100,      0,      0,      0),
  uppB1 = c( 40000,  20000,   5000,    100,      0,      0,      0,      0),
  lowB2 = c(  4000,   4000,   2000,    500,     10,      0,      0,      0),
  uppB2 = c(  4000,   2000,    500,     10,      0,      0,      0,      0),
  lowC  = c( 20000,  20000,  10000,   2500,    250,      0,      0,      0),
  uppC  = c( 20000,  10000,   2500,    250,      0,      0,      0,      0),
  lowD  = c(  2000,   2000,   1000,    250,     50,      0,      0,      0),
  uppD  = c(  2000,   1000,    250,     50,      0,      0,      0,      0),
  distr = c("unif", "unif", "unif", "unif", "decr", "unif", "unif", "unif"),
  beta  = c(    NA,     NA,     NA,     NA,     NA,     NA,     NA,     NA),
  stringsAsFactors = FALSE
)
```
The values in the dataframe are based on [IUCN (2012a)](https://portals.iucn.org/library/node/10315), 
[IUCN (2012b)](https://portals.iucn.org/library/node/10336), 
[IUCN (2022)](https://www.iucnredlist.org/resources/redlistguidelines) and the 
Norwegian guidance document ([Artsdatabanken 2020](https://artsdatabanken.no/Files/41216/)). 
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 extinction (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" provide the lower and upper threshold values for extinction time in _generations_ according to IUCN Red List Criterion E.
* The columns "lowA1" and "uppA1" provide the lower and upper threshold values for population reduction according to IUCN Red List Criterion A1.
* The columns "lowA2" and "uppA2" provide the lower and upper threshold values for population reduction according to IUCN Red List Criterion A2.
* 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 population size according to IUCN Red List Criterion C.
* The columns "lowD" and "uppD" provide the lower and upper threshold values for population size according to IUCN Red List Criterion D.
* 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 species and for species 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 is followed by this symbol, it is assumed to have been 
_downlisted_ by _one_ Red List Category.
```{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 frames 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("neglprop", "minority", "majority", "wholepop", "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", "wholepop", "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 population affected 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 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 population size over 10 years or
3 generations (whichever is largest) 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 Species, 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 species of Least Concern
extinct <- RLcateg$name[RLcateg$EX]  # abbreviation(s) for extinct species
LC.EX   <- RLcateg$name              # Red List Categories of evaluated species
RedListCat <- c(LC.EX, DD, notEval)  # all Red List Categories
```



## Read and check the data

Read the dataset "Norwegian Red List for species":
```{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)
}
```

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}
# Take a backup of the downlisted RL Categories in 2012
RL$Cat21orig <- downlist(RL)$Categ21
# Then reverse downlisting
RL <- uplist(RL) 
```
This step is only needed if the dataset contains species that have been downlisted 
due to rescue effects in other countries. Undoing this downlisting is motivated 
by the wish to quantify threats that can be addressed by management authorities 
in a given country.

If you want to analyse the data with downlisting retained, use instead the command 
`RL <- downlist(RL)`.

**(2) Back-cast** knowledge from the most recent Red List to earlier ones:
```{r}
RL <- backCast(RL)
```
This step is only needed if the dataset contains Red List Categories from 
different Red List Assessments. It corrects earlier assessments for the 
best available (i.e. most recent) information.

**(3) Calculate extinction probabilites** for all species:
```{r}
RL <- calcLoss(RL)
```

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



## Summarise the data

Summarise the Red Lists:
```{r}
tab <- summariseRL(RL)
```

In the specific case of the three Norwegian Red Lists analysed, the dataset only 
contains the species included in the current Red List (2021). For different 
reasons (such as taxonomic change), the earlier Red Lists contained species (names) 
that are not included in the most recent one. Therefore, the above summary table 
is not entirely correct for the earlier Red Lists. This has to be corrected 
manually by adding data for the Red Lists 2010 and 2015 (prior to back-casting).
The sources for these data are:

* [Artsdatabanken (2010)](http://www.artsportalen.artsdatabanken.no/)
* [Artsdatabanken (2015)](https://www.artsdatabanken.no/Rodlista2015) 

```{r}
# Create a new table for the results:
Table3 <- matrix(as.numeric(NA), 9, length(RedListCat) + 3, dimnames=list(
  "RL" %+% c("2010" %+% downlistSymbol, "2010", "2010(15)", "2010(21)",
             "2015" %+% downlistSymbol, "2015", "2015(21)", 
             "2021" %+% downlistSymbol, "2021"),
  c("N", RedListCat, "RLI", "Cum.ELS50")
))
alphabetic <- sort(RedListCat)

# Manually add the figures for the earlier Red Lists:
Table3[1, match(alphabetic, colnames(Table3))] <-
  c(284, 809, 890, 16762, 2580, 6528, 1310, 127, 1265)
Table3[2, match(alphabetic, colnames(Table3))] <-
  c(290, 809, 908, 16745,   NA,   NA, 1302, 127, 1266)
Table3[5, match(alphabetic, colnames(Table3))] <-
  c(247, 755, 901, 17594, 3018, 6095, 1302, 119, 1294)
Table3[6, match(alphabetic, colnames(Table3))] <- 
  c(252, 755, 916, 17579,   NA,   NA, 1297, 119, 1294)
tb <- table(RL$Cat21orig)
Table3[8, match(names(tb), colnames(Table3))] <- tb
Loss21o <- LoS(RL$Cat21orig, RL$GenTime)

# Insert the previous summary into this table:
Table3[c(4, 7, 9), colnames(tab)] <- tab[c(3, 5, 6), ]
Table3[3, ] <- Table3[4, ] - Table3[7, ] + Table3[6, ]
Table3[, "N"] <- apply(Table3[, RedListCat], 1, sum, na.rm=T)
Table3[, "RLI"] <- 1 - 
  apply(t(Table3[, LC.EX]) * RLW(LC.EX), 2, sum, na.rm=T) / 
  apply(  Table3[, LC.EX], 1, sum, na.rm=T) / max(RLW(LC.EX))
Table3[c(4, 7, 9), colnames(tab)] <- tab[c(3, 5, 6), ]
Table3 <- Table3[, !is.na(apply(Table3 > 0, 2, any))]
RLI21 <- RLI(RL$Categ21.21, RL$GenTime)
RLI15 <- RLI(RL$Categ15.21, RL$GenTime)
RLI10 <- RLI(RL$Categ10.21, RL$GenTime)

# Calculate means per Red List Category
# (needed to approximate species loss for data 
#  that are not based on the 2021 Red List):
mn10 <- mn15 <- rep(0, length(RedListCat))
names(mn10) <- names(mn15) <- RedListCat
for (i in RedListCat) {
  mn10[i] <- mean(RL$Loss10.21[which(RL$Categ10.21 == i)])
  mn15[i] <- mean(RL$Loss15.21[which(RL$Categ15.21 == i)])
}
Table3[1, "Cum.ELS50"] <- sum(mn10 * Table3[1, RedListCat], na.rm=T)
Table3[2, "Cum.ELS50"] <- sum(mn10 * Table3[2, RedListCat], na.rm=T)
Table3[3, "Cum.ELS50"] <- sum(mn15 * Table3[3, RedListCat], na.rm=T)
Table3[5, "Cum.ELS50"] <- sum(mn15 * Table3[5, RedListCat], na.rm=T)
Table3[6, "Cum.ELS50"] <- sum(mn15 * Table3[6, RedListCat], na.rm=T)
Loss21o[which(isDD(RL$Cat21orig))] <- RL$Loss21[which(isDD(RL$Cat21orig))]
Table3[8, "Cum.ELS50"] <- sum(Loss21o, na.rm=T)
rm(alphabetic, tb, tab, mn10, mn15, Loss21o)
```

Print the corrected table (which underlies Table 3 of the paper):
```{r}
print(Table3)
```



## Analysis of threat factors

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

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

Confidence intervals on RLI:
```{r, eval=FALSE, echo=TRUE}
print(confidenceRLI(RL, nsim, "Categ21"))
```
```{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(conf)
```

Confidence intervals on ΔRLI, δRLI and ELS~50~:
```{r, eval=FALSE, echo=TRUE}
results <- simulateDRLI(RL, nsim)
```
```{r, eval=TRUE, echo=FALSE}
# With the default `nsim`, the above line would take a looong time... 
# Instead, cached results are loaded:
load("cache.rdata")
for (i in 1:length(introductions)) {
  cat("\n\nConfidence intervals for " %+% introductions[i] %+% ":\n")
  print(simul[[i]])
}
```



## Figures

### Figure 1

The following script recreates Figure 1.

Simplify the table by collapsing minor threats:
```{r}
DRLI. <- DRLI
DRLI.[, 1]  <- 0
DRLI.[10, ] <- apply(DRLI.[c(1, 2, 5, 7, 9, 10, 11, 13), ], 2, sum)
DRLI. <- DRLI.[c(6, 10, 3, 12, 8, 4), ]
DRLI. <- DRLI.[, c(1:3, 3)]
DRLI. <- rbind(0, DRLI.)
DRLI. <- rbind(DRLI., 0)
```

Plot a graph for ΔRLI:
```{r}
{
  xl <- c(2009.5, 2028)
  yl <- c(0.9198, 0.92023)
  if (nchar(fig1)) {
    png(fig1, 1500, 1200, res = 180)
    xl <- c(2009.5, 2027.5)
    yl <- c(0.9198, 0.92022)
  }
  par(mai = c(0.96, 0.96, 0.12, 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(1, c(2009, 2021.5),              F, T, tcl = 0,        lwd = 1.5, lend = 1)
  axis(1, 2009:2021,                    F, T,                 lwd = 1.5, lend = 1)
  axis(1, c(2010, 2015, 2021),          T, T, cex.axis = 1.2, lwd = 1.5, lend = 1)
  axis(2, seq(0.9198, 0.9202, 0.00010), T, T, cex.axis = 1.2, lwd = 1.5, lend = 1)
  axis(2, seq(0.9198, 0.9202, 0.00001), F, T, tcl = -0.25,    lwd = 1.5, lend = 1)
  mtext("Year of Red List assessment", 1, 3, F, 2015.5, cex = 1.8)
  x <- c(2010, 2015, 2021, 2023)
  DRLI.[, 4] <- DRLI.[, 3] + c(0, 0, 0, 0, 0, 0, -0.00000887, 0)
  lines(x[4:1], rep(RLI10, 4), lty = "12",     lwd = 9.6, col = grey(0.84))
  lines(x[1:4], c(RLI10, RLI15, RLI21, RLI21), lwd = 9.6, col = grey(0.84))
  for (i in 2:7) {
    lines(x, RLI10 + DRLI.[i, ], lty = i - 1, lwd = 2.4)
    points(x[2:3], RLI10 + DRLI.[i, 2:3], pch = c(1, 4, 21, 24, 22, 25, 23)[i],
      cex = 1.8, bg = "black", lwd = 2.4, ljoin = 1)
    text(2023, RLI10 + DRLI.[i, 4],
      c("", "land-use change", "other/unknown", "climate change", "pollution",
        "native species", "disturbance")[i],
      pos = 4, cex = 1.2)
  }
  text(2023, RLI21, "RLI", pos = 4, cex = 1.2)
  text(2023, RLI10, "no change", pos = 4, cex = 1.2)
  text(2015.5, 0.92021, expression(bold(Delta*RLI)), cex = 1.8)
  if (nchar(fig1)) {
    dev.off()
  }
}
```



### Figure 2

The following script recreates Figure 2.

Simplify the table by collapsing minor threats:
```{r}
drli. <- drli$dRLI
ELS. <- drli$ELS50
for (i in 1:3) {
  drli.["otherthr", i] <- sum(drli.[c("otherthr", "unknownf", "alienspe",
    "huntgath", "outsiden", "natcatas", "bycatchc", "nothreat"), i])
  ELS. ["otherthr", i] <- sum(ELS. [c("otherthr", "unknownf", "alienspe",
    "huntgath", "outsiden", "natcatas", "bycatchc", "nothreat"), i])
}
drli. <- drli.[-which(rownames(drli.) %in%
  c("unknownf","alienspe","huntgath","outsiden","natcatas","bycatchc","nothreat")),]
ELS.  <- ELS. [-which(rownames(ELS.)  %in%
  c("unknownf","alienspe","huntgath","outsiden","natcatas","bycatchc","nothreat")),]
ELS.  <- ELS. [order(drli.[, 3], decreasing=T),]
ELS.  <- ELS. [, c(1:3, 3)]
drli. <- drli.[order(drli.[, 3], decreasing=T),]
drli. <- drli.[, 3:1]
drli. <- rbind(0, drli.)
drli. <- rbind(0, drli.)
```

Plot a graph for δRLI:
```{r}
{
  xl <- c(2009.5, 2028)
  yl <- c(0.91, 1.008)
  if (nchar(fig2)) {
    png(fig2, 1500, 1200, res=180)
    xl <- c(2009.5, 2027.5)
    yl <- c(0.91, 1.002)
  }
  par(mai=c(0.96, 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(1, c(2009, 2021.5),     F, T, tcl=0,        lwd=1.5, lend=1)
  axis(1, 2009:2021,           F, T,               lwd=1.5, lend=1)
  axis(1, c(2010, 2015, 2021), T, T, cex.axis=1.2, lwd=1.5, lend=1)
  axis(2, seq(0.8, 1, 0.01),   F, T,               lwd=1.5, lend=1)
  axis(2, seq(0.8, 1, 0.02),   T, T, cex.axis=1.2, lwd=1.5, lend=1)
  mtext("Year of Red List assessment", 1, 3, F, 2015.5, cex=1.8)
  x <- c(2021, 2015, 2010)
  for (i in 3:nrow(drli.)) {
    polygon(c(x, rev(x)),
            1 - c(apply(drli.[1:(i-1),], 2, sum), rev(apply(drli.[1:i,], 2, sum))),
            border=NA, col=grey(1.32 - i * 0.12))
  }
  y0 <- rep(RLI21, 3)
  y2 <- 1 - (sum(drli.[1:8,1]) + sum(drli.[1:7,1])) / 2
  #y2 <- RLI21 + 0.005
  for (i in 7:3) {
    y1 <- 1 - apply(drli.[1:i,], 2, sum)
    lines(x, y1, lwd=2.4)
    lines(c(2021, 2023), c(mean(c(y0[1], y1[1])), y2), lwd=1.8)
    text(2023, y2, 
      c("", "land-use change", "other/unknown", "climate change",
        "pollution", "native species", "disturbance")[i],
      pos=4, cex=1.2)
    y0 <- y1
    y2 <- y2 + 0.005
  }
  lines(c(2021, 2023), rep(y2, 2), lwd=1.8)
  text(2023, y2, "land-use change", pos=4, cex=1.2)
  lines(x, rep(1, 3), lwd=4.8)
  lines(c(2021, 2023), rep(1, 2), lwd=1.8)
  text(2023, 1, "reference value", font=2, pos=4, cex=1.2)
  lines(x, c(RLI21, RLI15, RLI10), lwd=4.8)
  lines(c(2021, 2023), c(RLI21, y2 - 0.03), lwd=1.8)
  text(2023, y2 - 0.03, "RLI", font=2, pos=4, cex=1.2)
  text(2015.5, 1.005, expression(bold(delta*RLI)), cex = 1.8)
  if (nchar(fig2)) {
    dev.off()
  }
}
```



### Figure 3

The following script recreates Figure 3.

Plot a graph for ELS~50~:
```{r}
{
  xl <- c(2009.5, 2028)
  yl <- c(lg(8), 3.1)
  if (nchar(fig3)) {
    png(fig3, 1500, 1200, res = 180)
    xl <- c(2009.5, 2027.5)
    yl <- c(lg(8), 3)
  }
  par(mai = c(0.96, 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 = "Expected loss of species",
       bty = "n", cex.axis = 1.2, cex.lab = 1.8)
  axis(1, c(2009, 2021.5),          F, T, tcl = 0,        lwd = 1.5, lend = 1)
  axis(1, 2009:2021,                F, T,                 lwd = 1.5, lend = 1)
  axis(1, c(2010, 2015, 2021),      T, T, cex.axis = 1.2, lwd = 1.5, lend = 1)
  axis(2, 0:3, c(1, 10, 100, 1000), T,    cex.axis = 1.2, lwd = 1.5, lend = 1)
  axis(2, lg(c(2:9, seq(20, 90, 10), seq(200, 900, 100))),
                                    F, T, tcl = -0.25,    lwd = 1.5, lend = 1)
  mtext("Year of Red List assessment", 1, 3, F, 2015.5, cex = 1.8)
  x <- c(2010, 2015, 2021, 2023)
  ELS.[3, 4] <- 10^(1/3 * lg(ELS.[4, 4]) + 2/3 * lg(ELS.[6, 4]))
  ELS.[5, 4] <- 10^(2/3 * lg(ELS.[4, 4]) + 1/3 * lg(ELS.[6, 4]))
  for (i in 1:6) {
    lines (x,      lg(ELS.[i, ]), lty = i, lwd = 2.4)
    points(x[1:3], lg(ELS.[i, 1:3]), pch = c(4, 21, 24, 22, 25, 23)[i],
      cex = 1.8, bg = "black", lwd = 2.4, ljoin = 1)
    text(2023, lg(ELS.[i, 4]), 
      c("land-use change", "other/unknown", "climate change",
        "pollution", "native species", "disturbance")[i],
      pos = 4, cex = 1.2)
  }
  text(2015.5, 3.01, expression(bold(ELS[50])), cex = 1.8)
  if (nchar(fig3)) {
    dev.off()
  }
}
```



## Analysis with DD species excluded

If the above analyses have been run with Data Deficient species included, 
they can be re-run with Data Deficient species excluded:
```{r}
{
  includeDD <- FALSE
  RL. <- calcLoss(RL)
  RL. <- addThreats(RL.)
  cat("\nDeltaRLI excluding DD species:\n")
  DRLI. <- DeltaRLI(RL.)
  print(DRLI.)
  drli. <- dRLI(RL.)
  cat("\ndRLI excluding DD species:\n")
  print(drli.$dRLI)
  cat("\nELS50 excluding DD species:\n")
  print(drli.$ELS50)
  includeDD <- TRUE
}
```



## Analysis with unknown threats inferred

If the above analyses have been run with unknown threats as a separate category, 
they can be re-run with unknown threats distributed over known threats:
```{r}
{
  inferThreats <- TRUE
  RL. <- calcLoss(RL)
  RL. <- addThreats(RL.)
  cat("\nDeltaRLI with unknown threats inferred:\n")
  DRLI. <- DeltaRLI(RL.)
  print(DRLI.)
  drli. <- dRLI(RL.)
  cat("\ndRLI with unknown threats inferred:\n")
  print(drli.$dRLI)
  cat("\nELS50 with unknown threats inferred:\n")
  print(drli.$ELS50)
  inferThreats <- FALSE
}
```



### Appendix S7
The following script recreates Appendix S7 (Figure S1).

Simplify the table by collapsing minor threats:
```{r}
DRLI.[,  1] <- 0
DRLI.[10, ] <- apply(DRLI.[c(1, 2, 5, 7, 9, 10, 11, 13), ], 2, sum)
DRLI. <- DRLI.[c(6, 10, 3, 12, 8, 4), ]
DRLI. <- DRLI.[, c(1:3, 3)]
DRLI. <- rbind(0, DRLI.)
DRLI. <- rbind(DRLI., 0)
```

Plot a graph for ΔRLI:
```{r}
{
  xl <- c(2009.5, 2028)
  yl <- c(0.9198, 0.92034)
  if (nchar(figS1)) {
    png(figS1, 1500, 1200, res = 180)
    xl <- c(2009.5, 2027.5)
    yl <- c(0.9198, 0.92033)
  }
  par(mai = c(0.96, 0.96, 0.12, 0.06), family = "sans") # ylim=c(0.9199, 0.92015)
  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(1, c(2009, 2021.5),              F, T, tcl = 0,        lwd = 1.5, lend = 1)
  axis(1, 2009:2021,                    F, T,                 lwd = 1.5, lend = 1)
  axis(1, c(2010, 2015, 2021),          T, T, cex.axis = 1.2, lwd = 1.5, lend = 1)
  axis(2, seq(0.9198, 0.9203, 0.00010), T, T, cex.axis = 1.2, lwd = 1.5, lend = 1)
  axis(2, seq(0.9198, 0.9203, 0.00001), F, T, tcl = -0.25,    lwd = 1.5, lend = 1)
  mtext("Year of Red List assessment", 1, 3, F, 2015.5, cex = 1.8)
  x <- c(2010, 2015, 2021, 2023)
  DRLI.[,4] <- DRLI.[,3] + c(0, 0, 0, 0, 0, 0, -0.00000876, 0)
  lines(x[3:1], rep(RLI10, 3), lty = "12", lwd = 9.6, col = grey(0.84))
  lines(x[1:3], c(RLI10, RLI15, RLI21),    lwd = 9.6, col = grey(0.84))
  for (i in 2:7) {
    lines(x, RLI10 + DRLI.[i, ], lty = i - 1, lwd = 2.4)
    points(x[2:3], RLI10 + DRLI.[i, 2:3], pch=c(1, 4, 21, 24, 22, 25, 23)[i],
      cex = 1.8, bg = "black", lwd = 2.4, ljoin = 1)
    text(2023, RLI10 + DRLI.[i, 4],
      c("", "land-use change", "other/unknown", "climate change", "pollution",
        "native species", "disturbance")[i],
      pos = 4, cex = 1.2)
  }
  text(2015.5, 0.92032, expression(bold(Delta*RLI)), cex = 1.8)
  if (nchar(figS1)) {
    dev.off()
  }
}
```



### Appendix S8
The following script recreates Appendix S8 (Figure S2).

Simplify the table by collapsing minor threats:
```{r}
ELS. <- drli.$ELS50
drli. <- drli.$dRLI
for (i in 1:3) {
  drli.["otherthr", i] <- sum(drli.[c("otherthr", "unknownf", "alienspe", 
    "huntgath", "outsiden", "natcatas", "bycatchc", "nothreat"), i])
  ELS. ["otherthr", i] <- sum(ELS. [c("otherthr", "unknownf", "alienspe",
    "huntgath", "outsiden", "natcatas", "bycatchc", "nothreat"), i])
}
drli. <- drli.[-which(rownames(drli.) %in%
  c("unknownf", "alienspe", "huntgath", "outsiden", 
    "natcatas", "bycatchc", "nothreat")),]
ELS.  <- ELS. [-which(rownames(ELS.)  %in%
  c("unknownf", "alienspe", "huntgath", "outsiden",
    "natcatas", "bycatchc", "nothreat")),]
ELS. <-  ELS.[c("landusec", "otherthr", "climatec", 
                "pollutio", "nativesp", "disturba"), ]
ELS.  <- ELS. [, c(1:3, 3)]
drli. <- drli.[order(drli.[, 3], decreasing=T), ]
drli. <- drli.[, 3:1]
drli. <- rbind(0, drli.)
drli. <- rbind(0, drli.)
```

Plot a graph for δRLI:
```{r}
{
  xl <- c(2009.5, 2028)
  yl <- c(0.91, 1.008)
  if (nchar(figS2)) {
    png(figS2, 1500, 1200, res = 180)
    xl <- c(2009.5, 2027.5)
    yl <- c(0.91, 1.008)
  }
  par(mai = c(0.96, 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(1, c(2009, 2021.5),     F, T, tcl = 0,        lwd = 1.5, lend = 1)
  axis(1, 2009:2021,           F, T,                 lwd = 1.5, lend = 1)
  axis(1, c(2010, 2015, 2021), T, T, cex.axis = 1.2, lwd = 1.5, lend = 1)
  axis(2, seq(0.8, 1, 0.01),   F, T,                 lwd = 1.5, lend = 1)
  axis(2, seq(0.8, 1, 0.02),   T, T, cex.axis = 1.2, lwd = 1.5, lend = 1)
  mtext("Year of Red List assessment", 1, 3, F, 2015.5, cex = 1.8)
  x <- c(2021, 2015, 2010)
  for (i in 3:nrow(drli.)) {
    polygon(c(x, rev(x)),
            1 - c(apply(drli.[1:(i-1),], 2, sum), rev(apply(drli.[1:i,], 2, sum))),
            border = NA, col = grey(1.32 - i * 0.12))
  }
  y0 <- rep(RLI21, 3)
  y2 <- 1 - (sum(drli.[1:8, 1]) + sum(drli.[1:7, 1])) / 2
  for (i in 7:3) {
    y1 <- 1 - apply(drli.[1:i, ], 2, sum)
    lines(x, y1, lwd = 2.4)
    lines(c(2021, 2023), c(mean(c(y0[1], y1[1])), y2), lwd = 1.8)
    text(2023, y2, 
      c("", "land-use change", "other/unknown", "climate change",
        "pollution", "native species", "disturbance")[i],
      pos = 4, cex = 1.2)
    y0 <- y1
    y2 <- y2 + 0.005
  }
  lines(c(2021, 2023), rep(y2, 2), lwd = 1.8)
  text(2023, y2, "land-use change", pos = 4, cex = 1.2)
  lines(x, rep(1,3), lwd = 4.8)
  lines(x, c(RLI21, RLI15, RLI10), lwd = 4.8)
  text(2015.5, 1.005, expression(bold(delta*RLI)), cex = 1.8)
  if (nchar(figS2)) {
    dev.off()
  }
}
```



### Appendix S9
The following script recreates Appendix S9 (Figure S3).

Plot a graph for ELS~50~:
```{r}
{
  xl <- c(2009.5, 2028)
  yl <- c(lg(8), 3.25)
  if (nchar(figS3)) {
    png(figS3, 1500, 1200, res = 180)
    xl <- c(2009.5, 2027.5)
    yl <- c(lg(8), 3.1)
  }
  par(mai = c(0.96, 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="Expected loss of species",
    bty = "n", cex.axis = 1.2, cex.lab = 1.8)
  axis(1, c(2009, 2021.5),          F, T, tcl = 0,        lwd = 1.5, lend = 1)
  axis(1, 2009:2021,                F, T,                 lwd = 1.5, lend = 1)
  axis(1, c(2010, 2015, 2021),      T, T, cex.axis = 1.2, lwd = 1.5, lend = 1)
  axis(2, c(1, lg(1200)),           F, T, tcl = 0,        lwd = 1.5, lend = 1)
  axis(2, 0:3, c(1, 10, 100, 1000), T,    cex.axis = 1.2, lwd = 1.5, lend = 1)
  axis(2, lg(c(2:9, seq(20, 90, 10), seq(200, 900, 100))),
                                    F, T, tcl = -0.25,    lwd = 1.5, lend = 1)
  mtext("Year of Red List assessment", 1, 3, F, 2015.5, cex = 1.8)
  x <- c(2010, 2015, 2021, 2023)
  ELS.[2, 4] <- 10^(3/4 * lg(ELS.[4, 4]) + 1/4 * lg(ELS.[6, 4]))
  ELS.[5, 4] <- 10^(2/4 * lg(ELS.[4, 4]) + 2/4 * lg(ELS.[6, 4]))
  ELS.[3, 4] <- 10^(1/4 * lg(ELS.[4, 4]) + 3/4 * lg(ELS.[6, 4]))
  for (i in 1:6) {
    lines (x,      lg(ELS.[i, ]), lty = i, lwd = 2.4)
    points(x[1:3], lg(ELS.[i, 1:3]), pch = c(4, 21, 24, 22, 25, 23)[i],
      cex = 1.8, bg = "black", lwd = 2.4, ljoin = 1)
    text(2023, lg(ELS.[i, 4]), 
      c("land-use change", "other/unknown", "climate change",
        "pollution", "native species", "disturbance")[i],
      pos = 4, cex = 1.2)
  }
  text(2015.5, 3.12, expression(bold(ELS[50])), cex = 1.8)
  if (nchar(figS3)) {
    dev.off()
  }
}
```



## Sensitivity analysis

In a kind of sensitivity analysis it is possible to check how important the 
weighting scheme chosen is for the (ranking of the) estimates obtained. 
This is here tested for the Expected Loss of Species. The most relevant measure 
is the _relative_ importance of threats, so in addition to ELS values themselves 
we should look at the _fraction_ of the total loss of species 
attributable to the different threats. 
These fractions are directly comparable across weightings.
```{r}
for (meth in c("E", "Ev2", "Ev3", "equal-steps", "A1", "A2", "B", "C", "D")) {
  weightingELS <- meth
  RL. <- calcLoss(RL)
  RL. <- addThreats(RL.)
  drli <- sort(dRLI(RL.)$ELS50[, 3], decreasing = TRUE)
  drli <- cbind(drli, drli / sum(drli))
  colnames(drli) <- c("ELS" %+% TimeFrame, "fraction")
  drli <- rbind(drli, Cumulative = apply(drli, 2, sum))
  cat("\n\nWeighting scheme " %+% meth %+% ":\n")
  print(drli)
}
```