# 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)) 
  }