# Supplementary Code | Source code for the functions to simulate species extinctions in the R language for statistical computing
# Explanatory comments (#)
# Author: Isabel Donoso
# R version 3.6.1 (2019-07-05)
# RCode related to:
# Title: Downsizing of animal communities triggers stronger functional than structural decay in seed-dispersal networks
# Authors: Donoso, I., Sorensen, M.C., Blendinger, P.G., Kissling, W.D., Neuschulz, E.L., Mueller, T. and Schleuning, M.
#----------------------------------------------------------------------------------------------------------------------------------------#
# Function to calculate % bird sps lost for each network depending on number of species(nrows). It is implemented in the functions below
#----------------------------------------------------------------------------------------------------------------------------------------#
def_pct <- function (x) {
pct <- list()
for (i in 1:(nrow(x)-1)){
pct[[i]] <- ((c(1:(nrow(x)-1)))[i]/(nrow(x)-1))*100 # it divides (Rich/total rich) * 100
}
return(pct)
}
#-----------------------------------------------#
#--- (A) SIZE-STRUCTURED EXTINCTION SCENARIO ---#
#-----------------------------------------------#
# def_Scn function to estimate structural and functional changes for all the networks under a size-structured extinction scenario
#
# INPUT DATA:
# (I) x is a list with a dataframes per web, where each row represents a dispersal event. Columns are:
# (1) plant species ($plant);
# (2) bird species ($bird);
# (3) species body mass in kilograms ($EBodyMass_kg) and
# (4) seed dispersal distance for that specific dispersal event ($dispdist).
#
# (II) webs is a list containing the matrices corresponding to each network, where the columns were already sorted according to birds’ body
# mass: from the largest to smallest bird species.
#
# OUTPUT DATA:
# Def_losses is an array with two dimensions: 1) Bird richness; 2) eight columns with all the information needed to estimate structural and
# functional changes. Columns are:
# (1) the percentage of bird species lost ($Def_pct);
# (2) the value of seed dispersal distance corresponding to the 0.95quantile of the total community wide seed dispersal kernel ($95%);
# (3) functional changes after scaling ($95%) to tge value of the original full assemblage ($Change95);
# (4) the number of plant species lost ($NPlant_Loss);
# (5) the percentage of plant species lost ($Pct_Plant_Loss);
# (6) the number of interactions that are lost as a consequence of losing each specific bird species from the community ($Int_Loss);
# (7) the percentage of interaction loss ($Pct_Int_Loss);
# (8) the cumulative percentage of interactions that are lost ($Cum_inter_lost).
def_Scn <- function (x, webs) {
# Vector of birds’ names determining the extinction sequence (for seed dispersal distance, SDD)
bird_sort <- x [order(x$EBodymass_kg,decreasing=T),]
def_seq <- unique(bird_sort$bird)
# Vector of birds’ names determining the extinction sequence (for Nplants and NInter lost)
matrix <- webs
def_col <- colnames(matrix)
# Create the array to save all the information
Def_losses <- array(NA,dim = c(length(def_seq)+1,8))
dimnames(Def_losses)[[1]] <- paste("Bird_rich",length(def_seq):0,sep="")
dimnames(Def_losses)[[2]] <- c("Def_pct","95%", "Change95", "NPlant_Loss","Pct_Plant_Loss", "Int_Loss", "Pct_Int_Loss","Cum_inter_lost")
Def_losses[length(def_seq)+1,"95%"] <- 0 # last row would be 0 (no birds, no distance)
Def_losses[length(def_col)+1, "NPlant_Loss"] <- nrow(webs) # when bird richness == 0, all plant species should have been removed
Def_losses[1, "Int_Loss"] <- 0 # when all birds sps are still there, the number of interactions lost should be 0
Def_losses[1, "Pct_Int_Loss"] <- 0 # when all bird sps are still there, the % inter loss should be 0
for (k in 1:length(def_seq)) {
# for SDD
Def_losses[k,"95%"] <- quantile(bird_sort$dispdist, probs = 0.95, na.rm=T)
bird_sort <- bird_sort[-which(bird_sort$bird == def_seq[k]), ]
# for N Plants; N_Interactions & Pct_Int_Loss
tmp <- rowSums(matrix)
Def_losses[k,"NPlant_Loss"] <- length (which (tmp==0))
removed_col <- matrix [,which (colnames(matrix) == def_col[k]),drop=FALSE]
Def_losses[k+1,"Int_Loss"] <- sum(removed_col)
Def_losses[k+1,"Pct_Int_Loss"] <- (sum(removed_col)*100) / sum(webs)
matrix <- matrix [,-which (colnames(matrix) == def_col[k]),drop=FALSE]
}
# Estimate (from the output), the % of Def birds, the % of plant Lost, the Cummulative % of Int Loss and the of functional changes (to have a complete final table)
tmp <- unlist(def_pct(Def_losses))
Def_losses[,"Def_pct"] <- c(0,tmp)
Def_losses[,"Pct_Plant_Loss"] <- Def_losses[,"NPlant_Loss"] / max (Def_losses[,"NPlant_Loss"]) * 100
Def_losses[,"Cum_inter_lost"] <- cumsum(Def_losses[,"Pct_Int_Loss"])
Def_losses[,"Change95"] <- (Def_losses[1,"95%"] - Def_losses[,"95%"]) / Def_losses[1,"95%"] * 100
return(Def_losses)
}
# ------------------------------------------------------------------ #
# Apply the size-structured defaunation scenario to all the networks #
# ------------------------------------------------------------------ #
# list_SDD is the list containing the dataframes with all dispersal events (as rows) for each of the networks.
# allnet_sort is the list containing the matrices corresponding to the networks, where the columns (bird species) are sorted according to birds’ body mass: from the largest to smallest bird)
# DEF_ALL_f will be the list containing the final table (i.e. Def_losses, derived from the def_Scn function) with the structural and functional changes values corresponding to each value of % bird species lost (Def_pct) for each empirical network.
DEF_ALL_f <- mapply(def_Scn, x = list_SDD, webs = allnet_sort)
#--------------------------------------#
#--- (B) RANDOM EXTINCTION SCENARIO ---#
#--------------------------------------#
# The random_Scn function was created similar to the def_Scn. However in this case the output is the array ran_losses with a third extra
# dimension for the number of replicates (i.e. results from each of the 1000 repetitions)
random_Scn <- function (x, webs, Nrep=1000) {
birds<-levels(x$bird)
ran_losses<-array(NA,dim=c(length(birds)+1,8,Nrep))
dimnames(ran_losses)[[1]]<-paste("Bird_rich",length(birds):0,sep="")
dimnames(ran_losses)[[2]]<-c("Def_pct","95%","Change95", "NPlant_Loss", "Pct_Plant_Loss", "Int_Loss", "Pct_Int_Loss", "Cum_inter_lost")
dimnames(ran_losses)[[3]]<-paste("Rep",1:Nrep,sep="")
ran_losses[length(birds)+1,"95%",] <- 0 # last row would be 0 (no birds, no distance)
ran_losses[length(birds)+1, "NPlant_Loss", ] <- nrow(webs) # when bird richness == 0, all plant species should have been removed
ran_losses[1, "Int_Loss", ] <- 0 # when all sps are there, there inter loss should be 0
ran_losses[1, "Pct_Int_Loss", ] <- 0 # when all sps are there, there % inter loss should be 0
for (n in 1:Nrep) {
seq.ran <- sample(birds)
data_ran <- x[order(match(x$bird, seq.ran)),]
matrix <- webs # for numer of plants and interactions
for (k in 1:length (seq.ran)) {
# for SDD
ran_losses[k,"95%",n] <- quantile(data_ran$dispdist, probs = 0.95, na.rm=T)
data_ran <- data_ran[-which(data_ran$bird == seq.ran[k]), ]
# for the number of plants and interactions
tmp <- rowSums(matrix)
ran_losses[k,"NPlant_Loss",n] <- length (which (tmp==0)) # cumulative number of plants that are lost
removed_col <- matrix [,which (colnames(matrix) == seq.ran[k]),drop=FALSE]
ran_losses[k+1,"Int_Loss",n] <- sum(removed_col)
ran_losses[k+1,"Pct_Int_Loss",n] <- (sum(removed_col)*100) / sum(webs)
matrix <- matrix[,-which(colnames(matrix) == seq.ran[k]), drop=FALSE]
}
tmp <- unlist(def_pct(ran_losses))
ran_losses[,"Def_pct",n] <- c(0,tmp)
ran_losses[,"Pct_Plant_Loss",n] <- ran_losses[,"NPlant_Loss",n] / max (ran_losses[,"NPlant_Loss",n]) * 100
ran_losses[,"Cum_inter_lost",n] <- cumsum(ran_losses[,"Pct_Int_Loss",n])
ran_losses[,"Change95",n] <- ((ran_losses[1,"95%",n] - ran_losses[,"95%",n]) / ran_losses[1,"95%",n]) * 100
}
return(ran_losses)
}
# -------------------------------------------------------- #
# Apply the random extinction scenario to all the networks #
# -------------------------------------------------------- #
# Because the random_Scn function computes 1000 repetitions per network, we used the parallel package to reduce computational time.
require(parallel)
detectCores()
nWebs <- length(list_SDD) # equal to the number of networks (n=8)
# Create a cluster with nWebs
cl <- makeCluster(nWebs)
# Export the input data and the functions to the cluster
clusterExport(cl, c("list_SDD", "allnet_sort", "random_Scn", "def_pct"))
# Estimate for all the networks #
RAN_ALL <- parLapply(cl, 1:nWebs, function(x) (random_Scn(x = list_SDD[[x]], webs=allnet_sort[[x]], Nrep = 1000)))
stopCluster(cl)
# For each network, compute the mean across repetitions
RAN_ALL_MEAN_rep <- list()
for (i in 1:8){
RAN_ALL_MEAN_rep[[i]] <- as.data.frame(apply(RAN_ALL[[i]], 1:2, mean))
}