---
title: "The direct drivers of recent global anthropogenic biodiversity loss"
subtitle: "As submitted to Science Advances"
author: "Andy Purvis & Nicolas Titeux"
date: "20 September 2022"
output:
  
  html_document: 
    df_print: paged
    code_folding: hide
    toc: true
    toc_depth: 2
    toc_float: true
    number_sections: true
  pdf_document:
    df_print: kable
    toc: true
    toc_depth: 2
    number_sections: true
---

```{r setup, include=TRUE}
knitr::opts_chunk$set(echo = FALSE)
set.seed(11) # Facilitating reproduction of results
```

# Ranking of direct drivers using David's score
This code produces the analyses reported in Jaureguiberry, Titeux et al., The direct drivers of recent global anthropogenic biodiversity loss, as published in $\it{Science Advances}$. The data underpinning the analysis are all in the `data` folder.

The code was written in a series of functions and blocks, not all of which were retained in the analyses in the paper (hence the gaps in the numbering). Block 07 sets the analysis options, which can then be changed programmatically to tweak details of the analysis. The initial set-up runs all realms (terrestrial, freshwater and marine) and all regions (the IPBES regions), counting ties as 0.5 of a win for each driver. Each head-to-head comparison receives weights based on its geographic scope and how seldom the indicator is reported in the dataset (which is recalculated within each bootstrap replicate). Additionally, assessments that compare more drivers yield more head-to-head comparisons between driver pairs, so automatically also receive more weight.
  
There are five direct drivers recognised. Some of the data include a sixth driver category, Other, but head-to-head comparisons between Other and any of the recognised direct drivers are omitted. The David score matrices among the five recognised drivers is the same whether or not Other drivers are included (though the normalised David scores change).
  
Bootstrap analyses use, by default, 10,000 trials, which is sufficient to give adequate precision for p-values near 0.05.

As we did not find a method in the literature for testing formally whether two or more sets of normalised David's scores are significantly different from each other, we developed a randomisation test. Our strategy was to stack the assessments from each subset, together with a factor indicating which subset each row came from. The values in this factor are then randomised (shuffled) `nreps` times, to provide a null distribution of a suitable test statistic. As the test statistic, we have used the residual standard error ($\sigma$) of a linear model predicting `normDS` from `driver`, in a dataframe that pools the normalised David scores from the results of each subset; this choice is suitable however many subsets are being compared.
  
The test is one-tailed: under the alternate hypothesis, the observed $\sigma$ will be larger than the null distribution (because each subset will give much the same normalised David score for any given driver), so all of the 5% zone of rejection is in one tail.
  
Post-hoc tests were conducted likewise, between each pair of subsets, using Benjamini & Yakutieli's (2001) method to control the false discovery rate by adjusting p-values.

## Figure 1A: Overall driver dominance hierarchy, pooling all realms and regions

Overall, land/sea use change is the dominant driver of biodiversity change, with direct exploitation second, pollution third, climate change fourth and invasive alien species fifth. The steepness of the normalised David's scores far exceeds those expected in the absence of any real hierarchy, and many but far from all driver pairs differ significantly in their impact (controlling for multiple testing). 

```{r, run-overall, include=TRUE, message=FALSE, cache=FALSE}

source("01_FUNCTION-generate-bootstrap-replicate.r")
source("02_FUNCTION-analyse-long-data.r")
source("03a_FUNCTION-new-analyse-DS-mat.r")
source("04_FUNCTION-compare-subsets.r")
source("05_load-packages.r")
source("06_read-csv-files.r")
source("07_set-analysis-options-overall.r")
source("08_combine-and-select-data.r")
source("09_rank-and-prepare-long-data.r")
source("10_call-analysis-and-bootstrap.r")
```

```{r, summarise-overall, include=TRUE, cache=FALSE}
source("11_save-objects-and-plots.r")
cat("David's score results for overall analysis:\n")
print(observed.david.score)
observed.DS.overall <- observed.david.score
cat("\n")
cat(paste("Mean and 95% confidence intervals of David's score from", nreps, "replicates:\n"))
sb.overall <- summarise_bootstraps(david.score.bootstrap)
cat("Steepness of hierarchy and test of randomness:\n")
steepness(observed.david.score$david_matrix, nrand=10000)
cat("\n")

# Calculate all pairwise p-values, reporting both face-value and adjusted p-values
cat(paste("Significance testing of David's scores, based on", nreps, "replicates.\n"))
cat("p-values adjusted using the Benjamini & Yakutieli method:\n")
nd <- ifelse(exclude.other.drivers=="Yes", 5, 6)
p.names <- NULL
p.boot <- NULL
for (i in c(1:(nd-1))){
  for (j in c((i+1):nd)){
    p.names <- append(p.names, paste(names(david.score.bootstrap)[i], names(david.score.bootstrap)[j], sep="_v_"))
    diffs <- david.score.bootstrap[[i]]-david.score.bootstrap[[j]]
    two.tailed.p <- 2 * min(sum(diffs >= 0), sum(diffs <= 0))/length(diffs)
    p.boot <- append(p.boot, two.tailed.p)
  }
}
p.adj <- p.adjust(p.boot, method="BY")
ps <- data.frame(comparison=p.names, uncorrected.p=p.boot, adjusted.p=p.adj)
print(ps)
ps.overall <- ps

write.results.csv(observed.DS.overall, sb.overall, ps.overall, "Overall_analysis")
```


## Figure 2: Comparing driver dominance hierarchies among regions

```{r analyse-and-compare-regions, include=TRUE, cache=FALSE}

source("07_set-analysis-options-overall.r")  #Reload default
region.selection <- "Yes" 
regions <- c("Americas", "Africa", "Europe and Central Asia", "Asia-Pacific")
for (selected.region in regions){
  print(paste("Selected region:", selected.region))
  source("08_combine-and-select-data.r")
  source("09_rank-and-prepare-long-data.r")
  source("10_call-analysis-and-bootstrap.r")
  source("11_save-objects-and-plots.r")
  print(observed.david.score) # Will also put out the attributes, so I can easily check the size of each data set
  cat("Steepness of hierarchy and test of randomness:\n")
  print(steepness(observed.david.score$david_matrix, nrand=10000))

  sb <- summarise_bootstraps(david.score.bootstrap)
  p.names <- NULL
  p.boot <- NULL
  for (i in c(1:(nd-1))){
    for (j in c((i+1):nd)){
      p.names <- append(p.names, paste(names(david.score.bootstrap)[i], names(david.score.bootstrap)[j], sep="_v_"))
      diffs <- david.score.bootstrap[[i]]-david.score.bootstrap[[j]]
      two.tailed.p <- 2 * min(sum(diffs > 0), sum(diffs < 0))/length(diffs)
      p.boot <- append(p.boot, two.tailed.p)
    }
  }
  p.adj <- p.adjust(p.boot, method="BY")
  ps <- data.frame(comparison=p.names, uncorrected.p=p.boot, adjusted.p=p.adj)
  print(ps)
  write.results.csv(observed.david.score, sb, ps, selected.region)
  attplot.wide$source <- selected.region
  if (selected.region == "Americas"){
    regions.data <- attplot.wide
  }else{
    regions.data <- rbind(regions.data, attplot.wide)
  }
}

# Omnibus test
regions.data$source <- as.factor(regions.data$source)
regiondiffs <- compare_subsets(regions.data, what="among IPBES regions")
cat("Testing for differences in normalised David's scores among IPBES regions\n")
print(regiondiffs[-4])

# Planned comparison of the two most data-rich regions
top.2 <- names(rev(sort(table(regions.data$source))))[1:2]
top2.regions.data <- droplevels(subset(regions.data, source %in% top.2))
top2regiondiffs <- compare_subsets(top2.regions.data, what="between Asia-Pacific and ECA")
cat("Testing for differences in normalised David's scores beetween two most data-rich IPBES regions\n")
print(top2regiondiffs[-4])

```


## Figure 3: Comparing driver dominance hierarchies among EBV classes

First each EBV class is analysed in turn. Then consistency among EBV classes is tested, first among all six classes, then dropping Genetic composition as it has very few data

### Dominance hierarchy of each EBV class in turn

```{r, run_for_each_ebv_class, include=TRUE, message=FALSE, cache=FALSE}
source("07_set-analysis-options-overall.r")  #Reset to defaults, such as overall analysis
analysis.level <- "ebv" # selecting "indicators.x.ebv" produces identical results
source("08_combine-and-select-data.r")
source("09_rank-and-prepare-long-data.r")
ebvs <- levels(attplot.wide$EBV.class)
ebvdata <- attplot.wide
ebvlongdata <- attplot

for (webv in 1:length(ebvs)){
  # Starting the loop at 2 drops genetic composition, for which there are only two rows of data
  attplot.wide <- droplevels(subset(ebvdata, EBV.class == ebvs[webv]))
  attplot <- droplevels(subset(ebvlongdata, EBV.class == ebvs[webv]))
  cat(paste("David's score results for EBV class:", ebvs[webv], "\n"))
  cat(paste(ebvs[webv], ":", nrow(attplot.wide), "rows of data.\n"))
  observed.david.score <- analyse.DS.mat(attplot.wide, study.drivers, IR=indicator.reweights, summarise.evidence=TRUE, show.result=TRUE)
  
  print(observed.david.score) # Will also put out the attributes
  cat("\n")
  cat("Steepness of hierarchy and test of randomness:\n")
  print(steepness(observed.david.score$david_matrix, nrand=10000))
  cat("\n")

	if (exclude.other.drivers == "Yes"){
	  david.score.bootstrap <- data.frame(climate_change = rep(NA, nreps), direct_exploitation = rep(NA, nreps), ias = rep(NA, nreps), 
	                                      land_sea_use = rep(NA, nreps), pollution = rep(NA, nreps), n_rows = rep(NA, nreps),
	                                      n_different_studies = rep(NA, nreps))
	}else{
	  david.score.bootstrap <- data.frame(climate_change = rep(NA, nreps), direct_exploitation = rep(NA, nreps), ias = rep(NA, nreps), 
		    land_sea_use = rep(NA, nreps), pollution = rep(NA, nreps), other = rep(NA, nreps), n_rows = rep(NA, nreps), 
		    n_different_studies = rep(NA, nreps))
	}
  
  for (i in 1:nreps){
    # Get the bootstrap pseudodatasets
    this.rep <- one.bootstrap.replicate(to.bootstrap.long = attplot, to.bootstrap.wide = attplot.wide) #no need to set the studies argument

    # Do the analysis of head-to-head
    this.wide.rep.out <- analyse.DS.mat(this.rep$wide, study.drivers,
                                        IR=indicator.reweights, summarise.evidence=FALSE, show.result=FALSE)$david_scores 
    
    david.score.bootstrap$climate_change[i] <- this.wide.rep.out$normDS[this.wide.rep.out$ID=="CC"]
    david.score.bootstrap$direct_exploitation[i] <- this.wide.rep.out$normDS[this.wide.rep.out$ID=="DE"]
    david.score.bootstrap$ias[i] <- this.wide.rep.out$normDS[this.wide.rep.out$ID=="IAS"]
    david.score.bootstrap$land_sea_use[i] <- this.wide.rep.out$normDS[this.wide.rep.out$ID=="LU"]
    david.score.bootstrap$pollution[i] <- this.wide.rep.out$normDS[this.wide.rep.out$ID=="PO"]
    if (exclude.other.drivers == "No") david.score.bootstrap$other[i] <- this.wide.rep.out$normDS[this.wide.rep.out$ID=="OT"]
    david.score.bootstrap$n_rows[i] <- nrow(this.rep$wide)
    david.score.bootstrap$n_different_studies[i] <- length(unique(this.rep$wide$UI))
  }
  
  cat(paste("Mean and 95% confidence intervals of David's score from", nreps, "replicates:\n"))
  sb <- summarise_bootstraps(david.score.bootstrap)
  print(sb)
  
  # Calculate all pairwise p-values
  cat(paste("Significance testing of David's scores, based on", nreps, "replicates.\n"))
  cat("p-values adjusted using the Benjamini & Yakutieli method:\n")
  nd <- ifelse(exclude.other.drivers=="Yes", 5, 6)
  p.names <- NULL
  p.boot <- NULL
  for (k in c(1:(nd-1))){
    for (j in c((k+1):nd)){
      p.names <- append(p.names, paste(names(david.score.bootstrap)[k], names(david.score.bootstrap)[j], sep="_v_"))
      diffs <- david.score.bootstrap[[k]]-david.score.bootstrap[[j]]
      two.tailed.p <- 2 * min(sum(diffs >= 0), sum(diffs <= 0))/length(diffs)
      p.boot <- append(p.boot, two.tailed.p)
    }
  }
  p.adj <- p.adjust(p.boot, method="BY")
  ps <- data.frame(comparison=p.names, uncorrected.p=p.boot, adjusted.p=p.adj)
  names(ps)[1] <- paste(names(ps)[1], "_within_", make.names(ebvs[webv])) # To tell the tables apart
  print(ps)
  
  write.results.csv(observed.david.score, sb, ps, paste("EBV_", ebvs[webv], "_", indicator.reweights,sep=""))
}
```

### Test whether dominance hierarchies in different EBV classes differ significantly

```{r, compare_ebv_classes_overall, include=TRUE, cache=FALSE}

# Now repeat the above steps but using ebv as the analysis.level
source("07_set-analysis-options-overall.r")  #Default is the overall analysis
analysis.level <- "ebv" # selecting "indicators.x.ebv" produces identical results
source("08_combine-and-select-data.r")
source("09_rank-and-prepare-long-data.r")
ebvdata <- attplot.wide
ebvdata$source <- ebvdata$EBV.class
print("")
print(paste("analysis.level = ebv produces this many rows", nrow(ebvdata)))
table(ebvdata$source)

# First, compare dominance hierarchies across all six EBV classes
# N.B. Not doing pairwise tests because Genetic composition has so few data
ebvdiffs <- compare_subsets(ebvdata, what="among six EBV classes", pairwise=FALSE)
print(ebvdiffs[-4])

# Second, drop genetic composition and compare hierarchies across top 5 EBV classes
ebvdata <- subset(ebvdata, source != "Genetic\ncomposition")
ebvdata <- droplevels(ebvdata)
ebvdiffs <- compare_subsets(ebvdata, what="among five EBV classes")
print(ebvdiffs[-4])

# Third, planned comparison between two most data-rich EBV classes
top2ebvs <- names(rev(sort(table(ebvdata$source))))[1:2]
top2.ebvdata <- subset(ebvdata, source %in% top2ebvs)
top2.ebvdiffs <- compare_subsets(top2.ebvdata, what="between top two EBV classes")
print(top2.ebvdiffs[-4])
```

## Figure 1B: Calculate dominance hierarchy within each realm and test for significant differences

This code also tests for significant differences among EBV classes within each realm in turn.

### Terrestrial realm

Land/sea use is the dominant driver of terrestrial biodiversity change, with direct exploitation second and climate change only fifth.

```{r, terrestrial, include=TRUE, cache=FALSE}
# Analyse terrestrial data

source("07_set-analysis-options-overall.r")  #Reload default
realm.selection<-"Yes"
selected.realm <- "Terrestrial"
source("08_combine-and-select-data.r")
source("09_rank-and-prepare-long-data.r")

terr.data <- attplot.wide #For use when comparing realms later
source("10_call-analysis-and-bootstrap.r")
source("11_save-objects-and-plots.r")
sb.terrestrial <- summarise_bootstraps(david.score.bootstrap)

terr.DS <- observed.david.score$david_scores
cat("Steepness of hierarchy and test of randomness:\n")
print(steepness(observed.david.score$david_matrix, nrand=10000))
observed.DS.terrestrial <- observed.david.score

cat(paste("Significance testing of terrestrial David's scores, based on", nreps, "replicates.\n"))
cat("p-values adjusted using the Benjamini & Yakutieli method:\n")
p.names <- NULL
p.boot <- NULL
for (i in c(1:(nd-1))){
  for (j in c((i+1):nd)){
    p.names <- append(p.names, paste(names(david.score.bootstrap)[i], names(david.score.bootstrap)[j], sep="_v_"))
    diffs <- david.score.bootstrap[[i]]-david.score.bootstrap[[j]]
    two.tailed.p <- 2 * min(sum(diffs > 0), sum(diffs < 0))/length(diffs)
    p.boot <- append(p.boot, two.tailed.p)
  }
}
p.adj <- p.adjust(p.boot, method="BY")
ps <- data.frame(comparison=p.names, uncorrected.p=p.boot, adjusted.p=p.adj)
print(ps)
ps.terrestrial <- ps

write.results.csv(observed.DS.terrestrial, sb.terrestrial, ps.terrestrial, "Terrestrial_realm")

```

### Compare EBV classes within terrestrial realm

```{r, terrestrial_ebv_comparison, include=TRUE, cache=FALSE}

source("07_set-analysis-options-overall.r")  #Reload default
realm.selection<-"Yes"
selected.realm <- "Terrestrial"
analysis.level <- "indicators.x.ebv"
source("08_combine-and-select-data.r")
source("09_rank-and-prepare-long-data.r")
ebvdata <- attplot.wide
ebvdata$source <- ebvdata$EBV.class
ebvdiffs <- compare_subsets(ebvdata, what="among EBV classes: terrestrial")
print(ebvdiffs[-4])
```

### Freshwater realm

Land/sea use change is the dominant driver of freshwater biodiversity change, with pollution second and climate change fourth.
 
```{r, freshwater, include=TRUE, cache=FALSE}
source("07_set-analysis-options-overall.r")  #Reload default

realm.selection<-"Yes"
selected.realm <- "Freshwater"
source("08_combine-and-select-data.r")
source("09_rank-and-prepare-long-data.r")

fw.data <- attplot.wide #For use when comparing realms later
source("10_call-analysis-and-bootstrap.r")
source("11_save-objects-and-plots.r")
sb.freshwater <- summarise_bootstraps(david.score.bootstrap)

fw.DS <- observed.david.score$david_scores
cat("Steepness of hierarchy and test of randomness:\n")
print(steepness(observed.david.score$david_matrix, nrand=10000))
observed.DS.freshwater <- observed.david.score

cat(paste("Significance testing of freshwater David's scores, based on", nreps, "replicates.\n"))
cat("p-values adjusted using the Benjamini & Yakutieli method:\n")
p.names <- NULL
p.boot <- NULL
for (i in c(1:(nd-1))){
  for (j in c((i+1):nd)){
    p.names <- append(p.names, paste(names(david.score.bootstrap)[i], names(david.score.bootstrap)[j], sep="_v_"))
    diffs <- david.score.bootstrap[[i]]-david.score.bootstrap[[j]]
    two.tailed.p <- 2 * min(sum(diffs > 0), sum(diffs < 0))/length(diffs)
    p.boot <- append(p.boot, two.tailed.p)
  }
}
p.adj <- p.adjust(p.boot, method="BY")
ps <- data.frame(comparison=p.names, uncorrected.p=p.boot, adjusted.p=p.adj)
print(ps)
ps.freshwater <- ps
write.results.csv(observed.DS.freshwater, sb.freshwater, ps.freshwater, "Freshwater_realm")

```

### Compare EBV classes within freshwater

```{r, freshwater_ebv_comparison, include=TRUE, cache=FALSE}

source("07_set-analysis-options-overall.r")  #Reload default
realm.selection<-"Yes"
selected.realm <- "Freshwater"
analysis.level <- "indicators.x.ebv"
source("08_combine-and-select-data.r")
source("09_rank-and-prepare-long-data.r")
ebvdata <- attplot.wide
ebvdata$source <- ebvdata$EBV.class
ebvdiffs <- compare_subsets(ebvdata, what="among EBV classes: freshwater")
print(ebvdiffs[-4])
```


### Marine realm

Direct exploitation is domminant in the marine realm, with climate change second and land/sea use change only fourth.


```{r, marine, include=TRUE, cache=FALSE}

source("07_set-analysis-options-overall.r")  #Reload default
realm.selection<-"Yes"
selected.realm <- "Marine"
source("08_combine-and-select-data.r")
source("09_rank-and-prepare-long-data.r")
mar.data <- attplot.wide #For use when comparing realms later

source("10_call-analysis-and-bootstrap.r")
source("11_save-objects-and-plots.r")
sb.marine <- summarise_bootstraps(david.score.bootstrap)

mar.DS <- observed.david.score$david_scores
cat("Steepness of hierarchy and test of randomness:\n")
print(steepness(observed.david.score$david_matrix, nrand=10000))
observed.DS.marine <- observed.david.score

cat(paste("Significance testing of marine David's scores, based on", nreps, "replicates.\n"))
cat("p-values adjusted using the Benjamini & Yakutieli method:\n")
p.names <- NULL
p.boot <- NULL
for (i in c(1:(nd-1))){
  for (j in c((i+1):nd)){
    p.names <- append(p.names, paste(names(david.score.bootstrap)[i], names(david.score.bootstrap)[j], sep="_v_"))
    diffs <- david.score.bootstrap[[i]]-david.score.bootstrap[[j]]
    two.tailed.p <- 2 * min(sum(diffs > 0), sum(diffs < 0))/length(diffs)
    p.boot <- append(p.boot, two.tailed.p)
  }
}
p.adj <- p.adjust(p.boot, method="BY")
ps <- data.frame(comparison=p.names, uncorrected.p=p.boot, adjusted.p=p.adj)
print(ps)
ps.marine <- ps

write.results.csv(observed.DS.marine, sb.marine, ps.marine, "Marine_realm")

```

### Compare EBV classes within marine

```{r, marine_ebv_comparison, include=TRUE, cache=FALSE}

source("07_set-analysis-options-overall.r")  #Reload default
realm.selection<-"Yes"
selected.realm <- "Marine"
analysis.level <- "indicators.x.ebv"
source("08_combine-and-select-data.r")
source("09_rank-and-prepare-long-data.r")
ebvdata <- attplot.wide
ebvdata$source <- ebvdata$EBV.class
ebvdiffs <- compare_subsets(ebvdata, what="among EBV classes: marine")
print(ebvdiffs[-4])
```


### Testing for significantly different hierarchies among realms (terrestrial, freshwater, marine)

Realms differ significantly in the normalised David's scores for the respective drivers (randomization: p << 0.05). Post-hoc testing controls the false discovery rate by adjusting p-values given multiple non-independent tests.

```{r, compare-realms, include=TRUE, cache=FALSE}
terr.data$source <- "terrestrial"
fw.data$source <- "freshwater"
mar.data$source <- "marine"
realm.datasets <- list(terr.data, fw.data, mar.data)
realms.data <- Reduce(rbind, realm.datasets)
realms.data$source <- as.factor(realms.data$source)

realmdiffs <- compare_subsets(realms.data, what="among realms")

# This is highly significant, so make pairwise comparisons among realms

TvF.data <- subset(realms.data, source != "marine")
TvF.diffs <- compare_subsets(TvF.data, what="between terrestrial & freshwater")
TvM.data <- subset(realms.data, source != "freshwater")
TvM.diffs <- compare_subsets(TvM.data, what="between terrestrial & marine")
FvM.data <- subset(realms.data, source != "terrestrial")
FvM.diffs <- compare_subsets(FvM.data, what="between freshwater & marine")

uncorrected.p <- c(TvF.diffs$p.value, TvM.diffs$p.value, FvM.diffs$p.value)
names(uncorrected.p) <- c("TvF", "TvM", "FvM")
adjusted.p <- p.adjust(uncorrected.p, method="BY")
cat("Post-hoc tests for significant differences between pairs of realms\n")
cat("(False discovery rate controlled using Benjamini & Yakutieli's 2001 method)\n")
print(adjusted.p)
```


# Compare impacts among EBV classes within each region

```{r ebvs_within_regions, include=TRUE, cache=FALSE}

source("07_set-analysis-options-overall.r")  #Reload default
region.selection <- "Yes" 
regions <- c("Americas", "Africa", "Europe and Central Asia", "Asia-Pacific")
analysis.level <- "indicators.x.ebv"

for (selected.region in regions){
  source("08_combine-and-select-data.r")
  source("09_rank-and-prepare-long-data.r")
  ebvdata <- attplot.wide
  ebvdata$source <- ebvdata$EBV.class
  print(paste("Comparing EBVs within", selected.region))
  ebvdiffs <- compare_subsets(ebvdata, what=paste("among EBV classes:", selected.region))
  print(ebvdiffs[-4])
}

```

