diff --git a/.Rhistory b/.Rhistory index 845dfa5..f322df2 100644 --- a/.Rhistory +++ b/.Rhistory @@ -1,7 +1,3 @@ -pkgdown::build_articles() -knitr::opts_chunk$set( -collapse = TRUE, -comment = "#>", warning = FALSE ) knitr::opts_chunk$set(fig.width=6, fig.height=4, message = FALSE) @@ -509,4 +505,8 @@ library(NetworkExtinction) library(NetworkExtinction) library(NetworkExtinction) devtools::install_github("r-lib/usethis") -library(NetworkExtinction) +plot(1:1e3, 1-pexp(1:1e3, rate = 1)) +plot(1:1e3, 1-pexp(1:1e3, rate = 100)) +plot(1:1e3, 1-pexp(1:1e3, rate = 0.0001)) +plot(1:1e3, 1-pexp(1:1e3, rate = 0.001)) +plot(1:1e3, 1-pexp(1:1e3, rate = 0.01)) diff --git a/.Rproj.user/shared/notebooks/paths b/.Rproj.user/shared/notebooks/paths index 7b7da71..92dfa09 100644 --- a/.Rproj.user/shared/notebooks/paths +++ b/.Rproj.user/shared/notebooks/paths @@ -1,2 +1,7 @@ -C:/Users/Erik Kusch/Desktop/[R Package] Network Extinction/R/Extintions.R="A28E7F0E" -D:/Documents/Work/[R Package] NetworkExtinction/R/Extintions.R="1D1BB781" +C:/Users/erike/Documents/Work/[Project] QuantiDoc Fish Health Score/PittmanIndex/PreShiny_Deploy.R="75B2FA30" +C:/Users/erike/Documents/Work/[Project] QuantiDoc Fish Health Score/PittmanIndex/Preamble.R="FB520A55" +C:/Users/erike/Documents/Work/[R Package] NetworkExtinction/NAMESPACE="2E8D7EF2" +C:/Users/erike/Documents/Work/[R Package] NetworkExtinction/R/Data.R="F2702D27" +C:/Users/erike/Documents/Work/[R Package] NetworkExtinction/R/NewExtinctionOrder.R="471A23E9" +C:/Users/erike/Documents/Work/[R Package] NetworkExtinction/R/NewSimulateExtinctions.R="BECBAC64" +C:/Users/erike/Documents/Work/[R Package] NetworkExtinction/vignettes/How_to_use_the_NetworkExtinction_Package.Rmd="4301D9B7" diff --git a/NAMESPACE b/NAMESPACE index d42d0f9..1f3a4ab 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -14,7 +14,6 @@ importFrom(doParallel,registerDoParallel) importFrom(dplyr,arrange) importFrom(dplyr,bind_rows) importFrom(dplyr,case_when) -importFrom(dplyr,desc) importFrom(dplyr,filter) importFrom(dplyr,full_join) importFrom(dplyr,group_by) @@ -50,7 +49,6 @@ importFrom(network,as.matrix.network.adjacency) importFrom(network,as.matrix.network.edgelist) importFrom(network,as.network) importFrom(network,delete.vertices) -importFrom(network,network.density) importFrom(network,network.edgecount) importFrom(network,network.size) importFrom(parallel,makeCluster) diff --git a/R/Data.R b/R/Data.R index ee4e726..f1ef6d6 100644 --- a/R/Data.R +++ b/R/Data.R @@ -6,6 +6,14 @@ #' "net" +#' A toymodel distance matrix +#' +#' A distance matrix used for demonstration of rewiring capabilities +#' +#' @format a distance matrix +#' +"dist" + #' The foodweb of the intertidal zone in central chile #' #' A trophic network with 107 species present in the intertidal zone of central Chile. @@ -26,8 +34,8 @@ #' A sparsely connected foodweb #' -#' A trophic network with 30 species and 47 trophic interactions. -#' This foodweb has a connectance of 0.03 +#' A network with 30 species and 47 interactions. +#' This network has a connectance of 0.03 #' @format a network diff --git a/R/Extintions.R b/R/Extintions.R index e0a9a58..b31ce84 100644 --- a/R/Extintions.R +++ b/R/Extintions.R @@ -1,37 +1,32 @@ -#' Extinctions analysis for trophic networks +#' Extinctions analysis for ecological networks #' #' The SimulateExtinctions function, can be used to test how the order of species -#' extinctions might affect the stability of the network by comparing The extinction history +#' extinctions, species-dependency on existing interaction strength, and rewiring potential might affect the stability of the network by comparing The extinction history #' and checking for secondary extinctions. #' -#' @param Method a character with the options Mostconnected, Oredered, or Random -#' @param Network a network representation as a an adyacency matrix, edgelist, +#' @param Network a network representation as a an adjacency matrix, edgelist, #' or a network object -#' @param Order this should be NULL, unless using the Ordered method, in that case -#' it should be a vector with the order of extinctions by ID +#' @param Method a character with the options Mostconnected and Ordered +#' @param Order a numeric vector indexing order of primary extinctions. For Method = Mostconnected Order must be NULL. If Order is not NULL, Method is internally forced to be Ordered. +#' @param NetworkType a character with the options Trophic and Mutualistic - is used to calculate secondary extinctions. #' @param clust.method a character with the options cluster_edge_betweenness, cluster_spinglass, #' cluster_label_prop or cluster_infomap, defaults to cluster_infomap -#' #' @param IS either numeric or a named vector of numerics. Identifies the threshold of relative interaction strength which species require before being considered secondarily extinct (i.e. IS = 0.7 leads to removal of all nodes which lose 30% of their interaction strength in the Network argument). If a named vector, names must correspond to vertex names in Network argument. +#' @param IS either numeric or a named vector of numerics. Identifies the threshold of relative interaction strength which species require before being considered secondarily extinct (i.e. IS = 0.7 leads to removal of all nodes which lose 30% of their interaction strength in the Network argument). If a named vector, names must correspond to vertex names in Network argument. +#' @param Rewiring either a function or a named vector of functions. Signifies how rewiring probabilities are calculated from the RewiringDist argument. If FALSE, no rewiring is carried out. +#' @param RewiringDist a numeric matrix of NxN dimension (N... number of nodes in Network). Contains, for example, phylogenetic or functional trait distances between nodes in Network which are used by the Rewiring argument to calculate rewiring probabilities. +#' @param RewiringProb a numeric which identifies the threshold at which to assume rewiring potential is met. #' @param verbose Logical. Whether to report on function progress or not. #' @return exports list containing a data frame with the characteristics of the network after every extinction and a network object containing the final network. The resulting data frame contains 11 columns that incorporate the topological index, the secondary extinctions, predation release, and total extinctions of the network in each primary extinction. +#' @details When method is Mostconnected, the function takes the network and calculates which node is the most connected of the network, using total degree. Then remove the most connected node, and calculates the the topological indexes of the network and the number of secondary extinctions. #' -#' @details When method is Mostconnected, it takes a network and it calculates which node is the most connected -#' of the network, using total degree. Then remove the most connected node, -#' and calculates the the topological indexes of the network and the number of -#' secundary extintions (how many species have indegree 0, without considered -#' primary producers). After that, remove the nodes that were secondarily extinct -#' in the previous step and recalculate which is the new most connected -#' node and so on, until the number of links in the network is zero. +#' When method is Ordered, it takes a network, and extinguishes nodes using a custom order, then it calculates the secondary extinctions and plots the accumulated secondary extinctions. #' -#' When method is Ordered, it takes a network, and extinguishes nodes using a custom order, -#' then it calculates the secondary extinctions and plots the accumulated -#' secondary extinctions. +#' When NetworkType = Trophic, secondary extinctions only occur for any predator, but not producers. If NetworkType = Mutualistic, secondary extinctions occur for all species in the network. #' #' When clust.method = cluster_edge_betweenness computes the network modularity using cluster_edge_betweenness methods from igraph to detect communities #' When clust.method = cluster_spinglass computes the network modularity using cluster_spinglass methods from igraph to detect communities, here the number of spins are equal to the network size #' When clust.method = cluster_label_prop computes the network modularity using cluster_label_prop methods from igraph to detect communities #' When clust.method = cluster_infomap computes the network modularity using cluster_infomap methods from igraph to detect communities, here the number of nb.trials are equal to the network size -#' #' @examples #' # Mostconnected example #' data("net") @@ -47,16 +42,39 @@ #' data("net") #' SimulateExtinctions(Network = net, Order = c(2,8,9), #' Method = "Ordered", clust.method = "cluster_infomap") - +#' +#' #Network-Dependency Example +#' data("net") +#' SimulateExtinctions(Network = net, Order = c(2,8), IS = 0.3, +#' Method = "Ordered", clust.method = "cluster_infomap") +#' +#' #Rewiring +#' data("net") +#' data(dist) +#' SimulateExtinctions(Network = net, Order = c(2,8), IS = 0.3, +#' Rewiring = function(x){1-pexp(x, rate = 1/0.5)}, # assuming an exponential decline in rewiring potential as values in RewiringDist increase +#' RewiringDist = dist, # distance matrix +#' RewiringProb = 0.2, # low threshold for rewiring potential +#' Method = "Ordered", clust.method = "cluster_infomap") +#' +#' #Rewiring, assuming dist contains probabilities +#' #' data("net") +#' data(dist) +#' SimulateExtinctions(Network = net, Order = c(2,8), IS = 0.3, +#' Rewiring = function(x){x}, # no changes to the RewiringDist object means +#' RewiringDist = dist, RewiringProb = 0.2, +#' Method = "Ordered", clust.method = "cluster_infomap") +#' #' @author Derek Corcoran #' @author M. Isidora Ávila-Thieme #' @author Erik Kusch #' @export - -SimulateExtinctions <- function(Network, Method, - Order = NULL, clust.method = "cluster_infomap", - IS = 0, verbose = TRUE){ +SimulateExtinctions <- function(Network, Method, Order = NULL, + NetworkType = "Trophic", clust.method = "cluster_infomap", + IS = 0, + Rewiring = FALSE, RewiringDist, RewiringProb = 0.5, + verbose = TRUE){ Network <- .DataInit(x = Network) if(!is.null(Order)){Method <- "Ordered"} @@ -69,328 +87,41 @@ SimulateExtinctions <- function(Network, Method, Conected <- data.frame(ID = 1:network.size(Network), Grado = degree(edgelist, c("total"))) Conected <- arrange(Conected, desc(Grado)) DF <- ExtinctionOrder(Network = Network, Order = Conected$ID, clust.method = clust.method, - IS = IS, verbose = verbose) + IS = IS, Rewiring = Rewiring, RewiringDist = RewiringDist, + verbose = verbose, RewiringProb = RewiringProb, NetworkType = NetworkType) } if(Method == "Ordered"){ DF <- ExtinctionOrder(Network = Network, Order = Order, clust.method = clust.method, - IS = IS, verbose = verbose) + IS = IS, Rewiring = Rewiring, RewiringDist = RewiringDist, + verbose = verbose, RewiringProb = RewiringProb, NetworkType = NetworkType) } return(DF) } -#' #' Extinctions analysis from most connected to less connected nodes in the network -#' #' -#' #' It takes a network and it calculates wich node is the most connected -#' #' of the network, using total degree. Then remove the most connected node, -#' #' and calculates the the topological indexes of the network and the number of -#' #' secundary extintions (how many species have indegree 0, without considered -#' #' primary producers). After that, remove the nodes that were secondarily extinct -#' #' in the previous step and recalculate which is the new most connected -#' #' node and so on, until the number of links in the network is zero. -#' #' -#' # -#' #' @param Network a trophic network of class network -#' #' @return exports data frame with the characteristics of the network after every -#' #' extintion. The resulting data frame contains 11 columns that incorporate the -#' #' topological index, the secondary extinctions, predation release, and total extinctions of the network -#' #' in each primary extinction. -#' #' @importFrom network as.matrix.network.edgelist -#' #' @importFrom network delete.vertices -#' #' @importFrom network network.density -#' #' @importFrom network network.edgecount -#' #' @importFrom network network.size -#' #' @importFrom sna degree -#' #' @importFrom stats complete.cases -#' #' @importFrom dplyr arrange -#' #' @importFrom dplyr desc -#' #' @author Derek Corcoran -#' #' @author M. Isidora Ávila-Thieme -#' #' @seealso [NetworkExtinction::ExtinctionOrder()] -#' #' @export -#' -#' -#' Mostconnected <- function(Network, IS = 0, verbose = TRUE){ -#' Grado <- NULL -#' Network <- .DataInit(x = Network) -#' edgelist <- as.matrix.network.edgelist(Network,matrix.type="edgelist") #Prey - Predator -#' Conected <- data.frame(ID = 1:network.size(Network), Grado = degree(edgelist, c("total"))) -#' Conected <- arrange(Conected, desc(Grado)) -#' Conected1<- c(Conected$ID) -#' indegreebasenet <- degree(Network, cmode = "indegree") -#' indegreebasenetzeros <- sum(degree(Network, cmode = "indegree") == 0) -#' indegreetopnetzeros <- sum(degree(Network, cmode = "outdegree") == 0) -#' Producers <- (1:length(degree(Network, cmode = "indegree")))[degree(Network, cmode = "indegree") == 0] -#' TopPredators <- (1:length(degree(Network, cmode = "outdegree")))[degree(Network, cmode = "outdegree") == 0] -#' DF <- data.frame(Spp = rep(NA, network.size(Network)), S = rep(NA, network.size(Network)), L = rep(NA, network.size(Network)), C = rep(NA, network.size(Network)), Link_density = rep(NA, network.size(Network)),SecExt = rep(NA,network.size(Network)), Pred_release = rep(NA,network.size(Network)), Iso_nodes =rep (NA,network.size(Network))) -#' -#' Secundaryext <- c() -#' Predationrel <- c() -#' accExt <- c() -#' totalExt <- c() -#' FinalExt <- list() -#' Conected3 <- c() -#' -#' ####LOOP#### -#' for (i in 1:network.size(Network)){ -#' #esta lista tiene el mismo orden que conected 1, hay que -#' #volver a hacer la red y calcular el grado -#' if (length(accExt)==0){ -#' Temp <- Network -#' DF$Spp[i] <- Conected1[i] -#' delete.vertices(Temp, c(DF$Spp[1:i])) -#' } -#' if (length(accExt)>0){ -#' Temp <- Network -#' Temp <- delete.vertices(Temp, c(accExt)) -#' edgelist <- as.matrix.network.edgelist(Temp,matrix.type="edgelist") -#' Conected2 <- data.frame(ID = 1:network.size(Temp), Grado = degree(edgelist, c("total"))) -#' Conected2 <- arrange(Conected2, desc(Grado)) -#' for(j in sort(accExt)){ -#' Conected2$ID <- ifelse(Conected2$ID < j, Conected2$ID, Conected2$ID + 1) -#' } -#' -#' DF$Spp[i] <- Conected2$ID[1] -#' Temp <- Network -#' -#' delete.vertices(Temp, unique(c(c(DF$Spp[1:i]),accExt))) -#' } -#' -#' DF$S[i] <- network.size(Temp) -#' DF$L[i] <- network.edgecount(Temp) -#' DF$C[i] <- network.density(Temp) -#' DF$Link_density [i] <- DF$L[i]/DF$S[i] -#' SecundaryextTemp <- (1:length(degree(Temp, cmode = "indegree")))[degree(Temp, cmode = "indegree") == 0] -#' for(j in sort(unique(c(c(DF$Spp[1:i]),accExt)))){ -#' SecundaryextTemp <- ifelse(SecundaryextTemp < j, SecundaryextTemp, SecundaryextTemp + 1) -#' } -#' Secundaryext <- SecundaryextTemp -#' Secundaryext <- Secundaryext[!(Secundaryext %in% Producers)] -#' DF$SecExt[i]<- length(Secundaryext) -#' -#' PredationrelTemp <- (1:length(degree(Temp, cmode = "outdegree")))[degree(Temp, cmode = "outdegree") == 0] -#' for(j in sort(unique(c(c(DF$Spp[1:i]),accExt)))){ -#' PredationrelTemp <- ifelse(PredationrelTemp < j, PredationrelTemp, PredationrelTemp + 1) -#' } -#' Predationrel <- PredationrelTemp -#' Predationrel <- Predationrel[!(Predationrel %in% TopPredators)] -#' DF$Pred_release[i]<- length(Predationrel) -#' -#' DF$Iso_nodes[i] <- sum(degree(Temp) == 0) -#' print(i) -#' FinalExt[[i]] <-(Secundaryext) -#' accExt <- append(accExt, DF$Spp[1:i]) -#' accExt <- unique(append(accExt,Secundaryext)) -#' -#' if (DF$L[i] == 0) break -#' } -#' DF <- DF[complete.cases(DF),] -#' DF$AccSecExt<- cumsum(DF$SecExt) -#' DF$NumExt <- 1:nrow(DF) -#' DF$TotalExt <- DF$AccSecExt + DF$NumExt -#' class(DF) <- c("data.frame", "Mostconnected") -#' .Deprecated("SimulateExtinctions") -#' return(DF) -#' } - -#' #' Extinctions analysis from most connected to less connected nodes in the network -#' #' -#' #' It takes a network and it calculates wich node is the most connected -#' #' of the network, using total degree. Then remove the most connected node, -#' #' and calculates the the topological indexes of the network and the number of -#' #' secundary extintions (how many species have indegree 0, without considered -#' #' primary producers). After that, remove the nodes that were secondarily extinct -#' #' in the previous step and recalculate which is the new most connected -#' #' node and so on, until the number of links in the network is zero. -#' #' -#' # -#' #' @param Network a trophic network of class network -#' #' @param clust.method a character with the options cluster_edge_betweenness, cluster_spinglass, cluster_label_prop or cluster_infomap -#' #' @param IS either numeric or a named vector of numerics. Identifies the threshold of relative interaction strength which species require before being considered secondarily extinct (i.e. IS = 0.7 leads to removal of all nodes which lose 30% of their interaction strength in the Network argument). If a named vector, names must correspond to vertex names in Network argument. -#' #' @param verbose Logical. Whether to report on function progress or not. -#' #' @return exports data frame with the characteristics of the network after every -#' #' extinction. The resulting data frame contains 11 columns that incorporate the -#' #' topological index, the secondary extinctions, predation release, and total extinctions of the network -#' #' in each primary extinction. -#' #' @importFrom network as.matrix.network.edgelist -#' #' @importFrom network delete.vertices -#' #' @importFrom network network.density -#' #' @importFrom network network.edgecount -#' #' @importFrom network network.size -#' #' @importFrom sna degree -#' #' @importFrom stats complete.cases -#' #' @importFrom dplyr arrange -#' #' @importFrom dplyr desc -#' #' @importFrom dplyr relocate -#' #' @importFrom network as.matrix.network.adjacency -#' #' @importFrom igraph graph_from_adjacency_matrix -#' #' @importFrom igraph cluster_edge_betweenness -#' #' @importFrom igraph cluster_spinglass -#' #' @importFrom igraph cluster_label_prop -#' #' @importFrom igraph cluster_infomap -#' #' @importFrom igraph modularity -#' #' @author Derek Corcoran -#' #' @author M. Isidora Ávila-Thieme -#' #' @seealso [NetworkExtinction::ExtinctionOrder()] -#' -#' .Mostconnected <- function(Network, clust.method = "cluster_infomap", -#' IS = IS, verbose = verbose){ -#' ## setting up objects for function run -#' Link_density <- Modularity <- Grado <- NULL -#' Network <- Network -#' edgelist <- as.matrix.network.edgelist(Network,matrix.type="edgelist") #Prey - Predator -#' Conected <- data.frame(ID = 1:network.size(Network), Grado = degree(edgelist, c("total"))) -#' Conected <- arrange(Conected, desc(Grado)) -#' Conected1<- c(Conected$ID) -#' if(length(IS )== 1){ -#' IS <- rep(IS, network.size(Network)) -#' names(IS) <- get.vertex.attribute(Network, "vertex.names") -#' } -#' -#' ## Base net calculations -#' ### identify base interaction strengths per node -#' if(sum(IS) != 0){ -#' if(sum(get.edge.attribute(Network, "weight"), na.rm = TRUE) == 0){ -#' stop("Either your network does not contain any edges with weights or your network does not have the edge attribute `weight` required for calculation of extinctions based on relative interaction strength loss.") -#' } -#' net <- as.matrix.network.adjacency(Network, attrname = "weight") -#' netgraph <- suppressMessages(graph_from_adjacency_matrix(net, weighted = TRUE)) -#' strengthbasenet <- igraph::strength(netgraph) -#' } -#' -#' ### identification of producers and top predators -#' indegreebasenet <- degree(Network, cmode = "indegree") -#' indegreebasenetzeros <- sum(degree(Network, cmode = "indegree") == 0) -#' indegreetopnetzeros <- sum(degree(Network, cmode = "outdegree") == 0) -#' Producers <- (1:length(degree(Network, cmode = "indegree")))[degree(Network, cmode = "indegree") == 0] -#' TopPredators <- (1:length(degree(Network, cmode = "outdegree")))[degree(Network, cmode = "outdegree") == 0] -#' -#' ## output object -#' DF <- data.frame(Spp = rep(NA, network.size(Network)), S = rep(NA, network.size(Network)), L = rep(NA, network.size(Network)), C = rep(NA, network.size(Network)), Link_density = rep(NA, network.size(Network)),SecExt = rep(NA,network.size(Network)), Pred_release = rep(NA,network.size(Network)), Iso_nodes =rep (NA,network.size(Network))) -#' Secundaryext <- c() -#' Predationrel <- c() -#' accExt <- c() -#' totalExt <- c() -#' FinalExt <- list() -#' Conected3 <- c() -#' -#' ## Sequential extinction simulation -#' if(verbose){ProgBar <- txtProgressBar(max = network.size(Network), style = 3)} -#' -#' ### creating temporary network representations and deleting vertices if they have been set to go extinct -#' for (i in 1:network.size(Network)){ -#' if (length(accExt)==0){ -#' Temp <- Network -#' DF$Spp[i] <- Conected1[i] -#' delete.vertices(Temp, c(DF$Spp[1:i])) -#' } -#' if (length(accExt)>0){ -#' Temp <- Network -#' Temp <- delete.vertices(Temp, c(accExt)) -#' edgelist <- as.matrix.network.edgelist(Temp,matrix.type="edgelist") -#' Conected2 <- data.frame(ID = 1:network.size(Temp), Grado = degree(edgelist, c("total"))) -#' Conected2 <- arrange(Conected2, desc(Grado)) -#' for(j in sort(accExt)){ -#' Conected2$ID <- ifelse(Conected2$ID < j, Conected2$ID, Conected2$ID + 1) -#' } -#' -#' DF$Spp[i] <- Conected2$ID[1] -#' Temp <- Network -#' -#' delete.vertices(Temp, unique(c(c(DF$Spp[1:i]),accExt))) -#' } -#' -#' ### network metrics to output object -#' DF$S[i] <- network.size(Temp) -#' DF$L[i] <- network.edgecount(Temp) -#' DF$C[i] <- network.density(Temp) -#' DF$Link_density [i] <- DF$L[i]/DF$S[i] -#' -#' if(i > 1 ){ -#' if(DF$L[i-1] == 0){ -#' if(verbose){setTxtProgressBar(ProgBar, length(Order))} -#' break -#' } -#' } -#' -#' -#' ### calculating modularity -#' Networkclass = class(Temp) -#' if (Networkclass[1] == "matrix"){ -#' netgraph = graph_from_adjacency_matrix(Temp, mode = "directed", weighted = TRUE) -#' } -#' -#' if (Networkclass[1] == "network"){ -#' net = as.matrix.network.adjacency(Temp) -#' netgraph = graph_from_adjacency_matrix(net, mode = "directed", weighted = TRUE) -#' } -#' -#' if (clust.method == "cluster_edge_betweenness"){ -#' Membership = suppressWarnings(cluster_edge_betweenness(netgraph, weights = TRUE, directed = FALSE, edge.betweenness = TRUE, -#' merges = TRUE, bridges = TRUE, modularity = TRUE, membership = TRUE)) -#' } else if (clust.method == "cluster_spinglass"){ -#' spins = network.size(Temp) -#' Membership = suppressWarnings(cluster_spinglass(netgraph, spins=spins)) #spins could be the Richness -#' }else if (clust.method == "cluster_label_prop"){ -#' Membership = suppressWarnings(cluster_label_prop(netgraph, weights = TRUE, initial = NULL, -#' fixed = NULL)) -#' }else if (clust.method == "cluster_infomap"){ -#' nb.trials = network.size(Temp) -#' Membership = suppressWarnings(cluster_infomap(as.undirected(netgraph), -#' e.weights = E(netgraph)$weight, -#' v.weights = NULL, -#' nb.trials = nb.trials, -#' modularity = TRUE)) -#' -#' } else stop('Select a valid method for clustering. ?SimulateExtinction') -#' DF$Modularity[i] <- Membership$modularity -#' -#' SecundaryextTemp <- (1:length(degree(Temp, cmode = "indegree")))[degree(Temp, cmode = "indegree") == 0] -#' for(j in sort(unique(c(c(DF$Spp[1:i]),accExt)))){ -#' SecundaryextTemp <- ifelse(SecundaryextTemp < j, SecundaryextTemp, SecundaryextTemp + 1) -#' } -#' Secundaryext <- SecundaryextTemp -#' Secundaryext <- Secundaryext[!(Secundaryext %in% Producers)] -#' DF$SecExt[i]<- length(Secundaryext) -#' -#' PredationrelTemp <- (1:length(degree(Temp, cmode = "outdegree")))[degree(Temp, cmode = "outdegree") == 0] -#' for(j in sort(unique(c(c(DF$Spp[1:i]),accExt)))){ -#' PredationrelTemp <- ifelse(PredationrelTemp < j, PredationrelTemp, PredationrelTemp + 1) -#' } -#' Predationrel <- PredationrelTemp -#' Predationrel <- Predationrel[!(Predationrel %in% TopPredators)] -#' DF$Pred_release[i]<- length(Predationrel) -#' -#' DF$Iso_nodes[i] <- sum(degree(Temp) == 0) -#' print(i) -#' FinalExt[[i]] <-(Secundaryext) -#' accExt <- append(accExt, DF$Spp[1:i]) -#' accExt <- unique(append(accExt,Secundaryext)) -#' -#' if (DF$L[i] == 0) break -#' } -#' DF <- DF[complete.cases(DF),] -#' DF$AccSecExt<- cumsum(DF$SecExt) -#' DF$NumExt <- 1:nrow(DF) -#' DF$TotalExt <- DF$AccSecExt + DF$NumExt -#' DF <- relocate(DF, Modularity, .after = Link_density) -#' class(DF) <- c("data.frame", "SimulateExt") -#' return(DF) -#' } - #' Extinctions analysis from custom order #' -#' It takes a network, and extinguishes nodes using a custom order, -#' then it calculates the secondary extinctions and plots the accumulated -#' secondary extinctions. +#' This function takes a network and eliminates nodes using a custom order. Subsequently, secondary extinctions are tallied up. Secondary extinction severity can be targeted by manipulating the node-dependency on network edges (IS) and node-rewiring potential upon loss of links (Rewiring). #' -#' @param Network a network of class network -#' @param Order Vector with the order of extinctions by ID +#' @param Network a network representation as a an adjacency matrix, edgelist, +#' or a network object +#' @param Method a character with the options Mostconnected and Ordered +#' @param Order a numeric vector indexing order of primary extinctions. For Method = Mostconnected Order must be NULL. If Order is not NULL, Method is internally forced to be Ordered. +#' @param NetworkType a character with the options Trophic and Mutualistic - is used to calculate secondary extinctions. +#' @param clust.method a character with the options cluster_edge_betweenness, cluster_spinglass, +#' cluster_label_prop or cluster_infomap, defaults to cluster_infomap #' @param IS either numeric or a named vector of numerics. Identifies the threshold of relative interaction strength which species require before being considered secondarily extinct (i.e. IS = 0.7 leads to removal of all nodes which lose 30% of their interaction strength in the Network argument). If a named vector, names must correspond to vertex names in Network argument. +#' @param Rewiring either a function or a named vector of functions. Signifies how rewiring probabilities are calculated from the RewiringDist argument. If FALSE, no rewiring is carried out. +#' @param RewiringDist a numeric matrix of NxN dimension (N... number of nodes in Network). Contains, for example, phylogenetic or functional trait distances between nodes in Network which are used by the Rewiring argument to calculate rewiring probabilities. +#' @param RewiringProb a numeric which identifies the threshold at which to assume rewiring potential is met. #' @param verbose Logical. Whether to report on function progress or not. -#' @param clust.method a character with the options cluster_edge_betweenness, cluster_spinglass, -#' cluster_label_prop or cluster_infomap #' @return exports list containing a data frame with the characteristics of the network after every extinction and a network object containing the final network. The resulting data frame contains 11 columns that incorporate the topological index, the secondary extinctions, predation release, and total extinctions of the network in each primary extinction. +#' @details When NetworkType = Trophic, secondary extinctions only occur for any predator, but not producers. If NetworkType = Mutualistic, secondary extinctions occur for all species in the network. +#' +#' When clust.method = cluster_edge_betweenness computes the network modularity using cluster_edge_betweenness methods from igraph to detect communities +#' When clust.method = cluster_spinglass computes the network modularity using cluster_spinglass methods from igraph to detect communities, here the number of spins are equal to the network size +#' When clust.method = cluster_label_prop computes the network modularity using cluster_label_prop methods from igraph to detect communities +#' When clust.method = cluster_infomap computes the network modularity using cluster_infomap methods from igraph to detect communities, here the number of nb.trials are equal to the network size #' #' @importFrom dplyr relocate #' @importFrom ggplot2 aes_string @@ -416,40 +147,66 @@ SimulateExtinctions <- function(Network, Method, #' @author M. Isidora Ávila-Thieme #' @author Erik Kusch #' @export - -ExtinctionOrder <- function(Network, Order, IS = 0, verbose = TRUE, - clust.method = "cluster_infomap"){ - ## Setting up Objects for function run +ExtinctionOrder <- function(Network, Order, NetworkType = "Trophic", clust.method = "cluster_infomap", + IS = 0, + Rewiring = FALSE, RewiringDist, RewiringProb = 0.5, + verbose = TRUE +){ + # Setting up Objects for function run ++++++++++ ++++++++++ ++++++++++ ++++++++++ Link_density <- Modularity <- Grado <- NULL Network <- .DataInit(x = Network) edgelist <- as.matrix.network.edgelist(Network,matrix.type="edgelist") #Prey - Predator Conected <- data.frame(ID = 1:network.size(Network), Grado = degree(edgelist, c("total"))) - Conected1 <- Order - if(length(IS )== 1){ - IS <- rep(IS, network.size(Network)) - names(IS) <- get.vertex.attribute(Network, "vertex.names") + Conected1 <- Order + + ## Interaction Strength Loss Preparations ++++++++++ ++++++++++ + if(length(IS )== 1){ # when the same dependency is to be applied for all species + IS <- rep(IS, network.size(Network)) # repeat IS argument for all species + names(IS) <- network::get.vertex.attribute(Network, "vertex.names") # assign species names to IS arguments + }else{ + if(sum(!(network::get.vertex.attribute(Network, "vertex.names") %in% names(IS))) != 0){stop("Please ensure that the names of the nodes in your Network are matched by names of the elements in your IS argument vector.")} + } + ## Rewiring Preparations ++++++++++ ++++++++++ + if(!isFALSE(Rewiring)){ # if rewiring has been specified + if(!exists("RewiringDist")){stop("To execute rewiring simulations, you need to specify the RewiringDist argument as well as the Rewiring argument.")} + diag(RewiringDist)<- NA # set diagonal to NA as one can't rewire to oneself + if(length(Rewiring) == 1){ # when the same Rewiring happen for all species + fun <- deparse1(Rewiring) # turn function into string + Rewiring <- rep(fun, network.size(Network)) # repeat function for each species + names(Rewiring) <- network::get.vertex.attribute(Network, "vertex.names") # assign names of species in network to rewiring functions + } + if(sum(!(network::get.vertex.attribute(Network, "vertex.names") %in% names(Rewiring))) != 0){stop("Please ensure that the names of the nodes in your Network are matched by names of the elements in your Rewiring argument vector.")} } - ## Base net calculations - ### identify base interaction strengths per node + # Base net calculations ++++++++++ ++++++++++ ++++++++++ ++++++++++ + ## base interaction strengths per node ++++++++++ ++++++++++ + options(warn=-1) # turn warning off + Weight_mat <- net <- as.matrix.network.adjacency(Network, attrname = "weight") + options(warn=0) # turn warnings on if(sum(IS) != 0){ - if(sum(get.edge.attribute(Network, "weight"), na.rm = TRUE) == 0){ - stop("Either your network does not contain any edges with weights or your network does not have the edge attribute `weight` required for calculation of extinctions based on relative interaction strength loss.") - } - net <- as.matrix.network.adjacency(Network, attrname = "weight") + # if(sum(network::get.edge.attribute(Network, "weight"), na.rm = TRUE) == 0){ + # stop("Either your network does not contain any edges with weights or your network does not have the edge attribute `weight` required for calculation of extinctions based on relative interaction strength loss.") + # } netgraph <- suppressMessages(graph_from_adjacency_matrix(net, weighted = TRUE)) strengthbasenet <- igraph::strength(netgraph) } - ### identification of producers and top predators + ## identification of producers and top predators ++++++++++ ++++++++++ indegreebasenet <- degree(Network, cmode = "indegree") indegreebasenetzeros <- sum(degree(Network, cmode = "indegree") == 0) indegreetopnetzeros <- sum(degree(Network, cmode = "outdegree") == 0) - Producers <- (1:length(degree(Network, cmode = "indegree")))[degree(Network, cmode = "indegree") == 0] - TopPredators <- (1:length(degree(Network, cmode = "outdegree")))[degree(Network, cmode = "outdegree") == 0] - - ### output object - DF <- data.frame(Spp = rep(NA, length(Order)), S = rep(NA, length(Order)), L = rep(NA, length(Order)), C = rep(NA, length(Order)), Link_density = rep(NA, length(Order)),SecExt = rep(NA,length(Order)), Pred_release = rep(NA,length(Order)), Iso_nodes =rep (NA,length(Order))) + Producers <- network::get.vertex.attribute(Network, "vertex.names")[degree(Network, cmode = "indegree") == 0] + TopPredators <- network::get.vertex.attribute(Network, "vertex.names")[degree(Network, cmode = "outdegree") == 0] + + ## output object ++++++++++ ++++++++++ + DF <- data.frame(Spp = rep(NA, length(Order)), + S = rep(NA, length(Order)), + L = rep(NA, length(Order)), + C = rep(NA, length(Order)), + Link_density = rep(NA, length(Order)), + SecExt = rep(NA,length(Order)), + Pred_release = rep(NA,length(Order)), + Iso_nodes = rep (NA,length(Order))) Secundaryext <- c() Predationrel <- c() accExt <- c() @@ -457,41 +214,33 @@ ExtinctionOrder <- function(Network, Order, IS = 0, verbose = TRUE, FinalExt <- list() Conected3 <- c() - ## Sequential extinction simulation + # Sequential extinction simulation ++++++++++ ++++++++++ ++++++++++ ++++++++++ if(verbose){ProgBar <- txtProgressBar(max = length(Order), style = 3)} for (i in 1:length(Order)){ # print(i) - - ### creating temporary network representations and deleting vertices if they have been set to go extinct - if (length(accExt)==0){ + ### creating temporary network + deleting vertices if they have been set to go extinct ++++++++++ ++++++++++ + if(length(accExt)==0){ # on first iteration Temp <- Network DF$Spp[i] <- Conected1[i] - delete.vertices(Temp, c(DF$Spp[1:i])) + network::delete.vertices(Temp, c(DF$Spp[1:i])) } - if (length(accExt)>0){ + if (length(accExt)>0){ # on any subsequent iteration Temp <- Network - Temp <- delete.vertices(Temp, c(accExt)) + Temp <- network::delete.vertices(Temp, c(accExt)) edgelist <- as.matrix.network.edgelist(Temp,matrix.type="edgelist") - # Conected2 <- data.frame(ID =1:network.size(Temp), Grado = degree(edgelist, c("total"))) - # for(j in sort(accExt)){ - # Conected2$ID <- ifelse(Conected2$ID < j, Conected2$ID, Conected2$ID + 1) - # } - DF$Spp[i] <- Conected1[i] Temp <- Network - - delete.vertices(Temp, unique(c(c(DF$Spp[1:i]),accExt))) - + network::delete.vertices(Temp, unique(c(c(DF$Spp[1:i]),accExt))) } - # print(Temp) - ### network metrics to output object + ## network metrics to output object ++++++++++ ++++++++++ DF$S[i] <- network.size(Temp) DF$L[i] <- network.edgecount(Temp) DF$C[i] <- network.density(Temp) DF$Link_density[i] <- DF$L[i]/DF$S[i] - if(i > 1 ){ + ## premature complete annihilation message ++++++++++ ++++++++++ + if(i > 1){ if(DF$L[i-1] == 0){ if(verbose){setTxtProgressBar(ProgBar, length(Order))} warning(paste("All species in network went extinct through secondary extinction before all primary extinctions were simulated. This happened at extinction step", i-1, "out of", length(Order))) @@ -499,18 +248,15 @@ ExtinctionOrder <- function(Network, Order, IS = 0, verbose = TRUE, } } - - ### calculating modularity + ## calculating modularity ++++++++++ ++++++++++ Networkclass = class(Temp) if (Networkclass[1] == "matrix"){ netgraph = graph_from_adjacency_matrix(Temp, mode = "directed", weighted = TRUE) } - if (Networkclass[1] == "network"){ net = as.matrix.network.adjacency(Temp) netgraph = suppressMessages(graph_from_adjacency_matrix(net, mode = "directed", weighted = TRUE)) } - if (clust.method == "cluster_edge_betweenness"){ Membership = suppressWarnings(cluster_edge_betweenness(netgraph, weights = TRUE, directed = TRUE, edge.betweenness = TRUE, merges = TRUE, bridges = TRUE, modularity = TRUE, membership = TRUE)) @@ -531,55 +277,110 @@ ExtinctionOrder <- function(Network, Order, IS = 0, verbose = TRUE, } else if (clust.method == "none"){ Membership = NA }else stop('Select a valid method for clustering. ?SimulateExtinction') - # if(is.na(Membership)){ DF$Modularity[i] <- Membership$modularity - # }else{ - # DF$Modularity[i] <- suppressWarnings(modularity(Membership)) - # } - ### identifying secondary extinctions - #### Producers - SecundaryextTemp <- (1:length(degree(Temp, cmode = "indegree")))[degree(Temp, cmode = "indegree") == 0] - for(j in sort(unique(c(c(DF$Spp[1:i]),accExt)))){ - SecundaryextTemp <- ifelse(SecundaryextTemp < j, SecundaryextTemp, SecundaryextTemp + 1) - } - Secundaryext <- SecundaryextTemp - Secundaryext <- Secundaryext[!(Secundaryext %in% Producers)] - DF$SecExt[i] <- length(Secundaryext) - - #### Predators - PredationrelTemp <- (1:length(degree(Temp, cmode = "outdegree")))[degree(Temp, cmode = "outdegree") == 0] - for(j in sort(unique(c(c(DF$Spp[1:i]),accExt)))){ - PredationrelTemp <- ifelse(PredationrelTemp < j, PredationrelTemp, PredationrelTemp + 1) + ## rewiring ++++++++++ ++++++++++ + accExt <- unique(append(accExt, DF$Spp[1:i])) + if(!isFALSE(Rewiring)){ + ### identify rewiring potential ++++++++++ + Rewiring_df <- data.frame(Direction = NA, + Species = NA, + NewPartner = NA, + LostPartner = NA, + IS = NA) + Rewiring_df <- na.omit(Rewiring_df) + #### loop over all deleted vertices and the connections lost because of their exclusion + for(Iter_PrimaryExt in 1:length(accExt)){ + # Iter_PrimaryExt = 1 + LostPartner <- network::get.vertex.attribute(Network, "vertex.names")[accExt[Iter_PrimaryExt]] # name of primary extinction species + LostISCol <- Weight_mat[, LostPartner] # lost interaction strength with nodes now slated for secondary extinction + LostISRow <- Weight_mat[LostPartner, ] + Lost_df <- data.frame(LostIS = c(LostISCol, LostISRow), + Direction = rep(c(1,2), c(length(LostISCol), length(LostISRow))), + names = c(names(LostISCol), names(LostISRow)) + ) + Lost_df <- Lost_df[Lost_df$LostIS != 0, ] + if(nrow(Lost_df)!=0){ + for(Iter_LostIS in 1:nrow(Lost_df)){ ## looping over all species that were linked to the current primary extinction + # Iter_LostIS = 1 + LostPartnerSim <- eval(str2lang(Rewiring[which(names(Rewiring) == Lost_df$names[Iter_LostIS])]))(RewiringDist[,LostPartner]) # probability of rewiring too each node in network given rewiring function and species similraity + RewiringCandidates <- LostPartnerSim[LostPartnerSim > RewiringProb & names(LostPartnerSim) %in% network::get.vertex.attribute(Temp, "vertex.names")] # rewiring probability for nodes still in temporary network and having a higher rewiring probability than 0.5 + RewiredPartner <- names(which.max(RewiringCandidates)) # most likely rewiring partner + if(!is.null(RewiredPartner)){ # if a rewired partner has been found + Rewiring_df <- rbind(Rewiring_df, + data.frame(Direction = Lost_df$Direction[Iter_LostIS], + Species = Lost_df$names[Iter_LostIS], + NewPartner = RewiredPartner, + LostPartner = LostPartner, + IS = Lost_df$LostIS[Iter_LostIS]) + ) + } + } + } + } + + ### shift rewired interaction strengths ++++++++++ + if(nrow(Rewiring_df) != 0){ + #### shift interaction weights in Weight_mat + for(Iter_Rewiring in 1:nrow(Rewiring_df)){ + # Iter_Rewiring = 1 + ## assigning shifted interaction strength + ColSpec <- Rewiring_df[Iter_Rewiring,4-Rewiring_df[Iter_Rewiring,"Direction"]] + RowSpec <- Rewiring_df[Iter_Rewiring,1+Rewiring_df[Iter_Rewiring,"Direction"]] + Weight_mat[RowSpec, ColSpec] <- Weight_mat[RowSpec, ColSpec] + Rewiring_df[Iter_Rewiring,"IS"] + ## deleting shiften interaction strength + + ColLost <- ifelse(Rewiring_df[Iter_Rewiring, "Direction"] == 1, + Rewiring_df[Iter_Rewiring, "LostPartner"], + Rewiring_df[Iter_Rewiring, "Species"]) + RowLost <- ifelse(Rewiring_df[Iter_Rewiring, "Direction"] == 1, + Rewiring_df[Iter_Rewiring, "Species"], + Rewiring_df[Iter_Rewiring, "LostPartner"]) + Weight_mat[RowLost, ColLost] <- 0 + } + #### establishing rewired network and deleting primary extinction nodes + Network <- as.network(Weight_mat, matrix.type = "adjacency", ignore.eval=FALSE, names.eval='weight') + Temp <- Network + network::delete.vertices(Temp, unique(c(c(DF$Spp[1:i]),accExt))) + } } - Predationrel <- PredationrelTemp - Predationrel <- Predationrel[!(Predationrel %in% TopPredators)] - DF$Pred_release[i]<- length(Predationrel) - DF$Iso_nodes[i] <- sum(degree(Temp) == 0) - #### Relative Interaction Strength loss + ## identify secondary extinctions ++++++++++ ++++++++++ + ### Relative Interaction Strength loss ++++++++++ if(sum(IS) == 0){ - Secundaryext <- get.vertex.attribute(Temp, "vertex.names")[which(degree(Temp) == 0)] - Secundaryext <- match(Secundaryext, get.vertex.attribute(Network, "vertex.names")) + SecundaryextNames <- network::get.vertex.attribute(Temp, "vertex.names")[which(degree(Temp) == 0)] + Secundaryext <- match(SecundaryextNames, network::get.vertex.attribute(Network, "vertex.names")) }else{ AbsIS <- igraph::strength(suppressMessages(graph_from_adjacency_matrix( as.matrix.network.adjacency(Temp, attrname = "weight"), weighted = TRUE) )) - RelISloss <- AbsIS / strengthbasenet[names(strengthbasenet) %in% get.vertex.attribute(Temp, "vertex.names")] - Secundaryext <- which(AbsIS == 0 | RelISloss < IS[match(names(RelISloss), names(IS))]) - Secundaryext <- match(names(Secundaryext), get.vertex.attribute(Network, "vertex.names")) + RelISloss <- AbsIS / strengthbasenet[names(strengthbasenet) %in% network::get.vertex.attribute(Temp, "vertex.names")] + SecundaryextNames <- which(AbsIS == 0 | (1-RelISloss) < IS[match(names(RelISloss), names(IS))]) + Secundaryext <- match(names(SecundaryextNames), network::get.vertex.attribute(Network, "vertex.names")) + } + + ### for trophic networks ++++++++++ + + if(NetworkType == "Trophic"){ + DF$SecExt[i] <- length(SecundaryextNames[!(names(SecundaryextNames) %in% as.character(Producers))]) + DF$Pred_release[i] <- length(SecundaryextNames[!(names(SecundaryextNames) %in% as.character(TopPredators))]) } - DF$SecExt[i] <- length(Secundaryext) + ### for mutualistic networks ++++++++++ + if(NetworkType == "Mutualistic"){ + DF$SecExt[i] <- length(Secundaryext) + } + DF$Iso_nodes[i] <- sum(degree(Temp) == 0) - ### Return of objects - FinalExt[[i]] <- (Secundaryext) + ## Return of objects ++++++++++ ++++++++++ + FinalExt[[i]] <- Secundaryext accExt <- append(accExt, DF$Spp[1:i]) accExt <- unique(append(accExt,Secundaryext)) if(verbose){setTxtProgressBar(ProgBar, i)} - } - DF <- DF[complete.cases(DF),] + + # return of final data objects ++++++++++ ++++++++++ + DF <- DF[!is.na(DF$Spp),] DF$AccSecExt <- cumsum(DF$SecExt) DF$NumExt <- 1:nrow(DF) DF$TotalExt <- DF$AccSecExt + DF$NumExt @@ -591,21 +392,34 @@ ExtinctionOrder <- function(Network, Order, IS = 0, verbose = TRUE, #' Random extinction #' -#' Generates a null model by generating random extinction histories and calculating -#' the mean and standard deviation of the accumulated secondary extinctions developed -#' by making n random extinction histories +#' Generates a null model by generating random extinction histories and calculating the mean and standard deviation of the accumulated secondary extinctions developed by making n random extinction histories. #' -#' @param Network a trophic network of class network -#' @param nsim number of simulations -#' @param IS either numeric or a named vector of numerics. Identifies the threshold of relative interaction strength which species require before being considered secondarily extinct (i.e. IS = 0.7 leads to removal of all nodes which lose 30% of their interaction strength in the Network argument). If a named vector, names must correspond to vertex names in Network argument. -#' @param verbose Logical. Whether to report on function progress or not. -#' @param parallel if TRUE, it will use parallel procesing, if FALSE (default) it will run -#' sequentially -#' @param ncores number of cores to use if using parallel procesing +#' @param Network a network representation as a an adjacency matrix, edgelist, +#' or a network object +#' @param nsim numeric, number of simulations #' @param Record logical, if TRUE, records every simulation and you can read the #' raw results in the object FullSims -#' @param plot logical if true, will add a graph to the results +#' @param plot logical if TRUE, will add a graph to the results +#' @param SimNum numeric, how many nodes to register for primary extinction. By default sets all of them. +#' @param NetworkType a character with the options Trophic and Mutualistic - is used to calculate secondary extinctions. +#' @param clust.method a character with the options cluster_edge_betweenness, cluster_spinglass, +#' cluster_label_prop or cluster_infomap, defaults to cluster_infomap +#' @param parallel if TRUE, it will use parallel procesing, if FALSE (default) it will run +#' sequentially +#' @param ncores numeric, number of cores to use if using parallel procesing +#' @param IS either numeric or a named vector of numerics. Identifies the threshold of relative interaction strength which species require before being considered secondarily extinct (i.e. IS = 0.7 leads to removal of all nodes which lose 30% of their interaction strength in the Network argument). If a named vector, names must correspond to vertex names in Network argument. +#' @param Rewiring either a function or a named vector of functions. Signifies how rewiring probabilities are calculated from the RewiringDist argument. If FALSE, no rewiring is carried out. +#' @param RewiringDist a numeric matrix of NxN dimension (N... number of nodes in Network). Contains, for example, phylogenetic or functional trait distances between nodes in Network which are used by the Rewiring argument to calculate rewiring probabilities. +#' @param RewiringProb a numeric which identifies the threshold at which to assume rewiring potential is met. +#' @param verbose Logical. Whether to report on function progress or not. #' @return exports list containing a data frame with the characteristics of the network after every extinction, a network object containing the final network, and a graph with the mean and 95% interval. The resulting data frame contains 11 columns that incorporate the topological index, the secondary extinctions, predation release, and total extinctions of the network in each primary extinction. +#' @details When NetworkType = Trophic, secondary extinctions only occur for any predator, but not producers. If NetworkType = Mutualistic, secondary extinctions occur for all species in the network. +#' +#' When clust.method = cluster_edge_betweenness computes the network modularity using cluster_edge_betweenness methods from igraph to detect communities +#' When clust.method = cluster_spinglass computes the network modularity using cluster_spinglass methods from igraph to detect communities, here the number of spins are equal to the network size +#' When clust.method = cluster_label_prop computes the network modularity using cluster_label_prop methods from igraph to detect communities +#' When clust.method = cluster_infomap computes the network modularity using cluster_infomap methods from igraph to detect communities, here the number of nb.trials are equal to the network size +#' #' @examples #' #first example #' data("More_Connected") @@ -643,11 +457,14 @@ ExtinctionOrder <- function(Network, Order, IS = 0, verbose = TRUE, #' @author M. Isidora Ávila-Thieme #' @author Erik Kusch #' @export - -RandomExtinctions <- function(Network, nsim = 10, parallel = FALSE, ncores, +RandomExtinctions <- function(Network, nsim = 10, Record = FALSE, plot = FALSE, SimNum = NULL, - IS = 0, verbose = TRUE){ + NetworkType = "Trophic", clust.method = "cluster_infomap", + parallel = FALSE, ncores, + IS = 0, + Rewiring = FALSE, RewiringDist, RewiringProb = 0.5, + verbose = TRUE){ ## setting up objects NumExt <- sd <- AccSecExt <- AccSecExt_95CI <- AccSecExt_mean <- Lower <- NULL network <- .DataInit(x = Network) @@ -661,7 +478,10 @@ RandomExtinctions <- function(Network, nsim = 10, parallel = FALSE, ncores, cl <- makeCluster(ncores) registerDoParallel(cl) sims <- foreach(i=1:nsim, .packages = "NetworkExtinction")%dopar%{ - sims <- try(ExtinctionOrder(Network = network, Order = sample(1:network.size(network), size = SimNum), IS = IS, verbose = FALSE), silent = TRUE) + sims <- try(ExtinctionOrder(Network = network, Order = sample(1:network.size(network), size = SimNum), + IS = IS, NetworkType = NetworkType, + Rewiring = Rewiring, RewiringDist = RewiringDist, + verbose = FALSE, RewiringProb = RewiringProb), silent = TRUE) try({sims$simulation <- i}, silent = TRUE) sims } @@ -669,7 +489,10 @@ RandomExtinctions <- function(Network, nsim = 10, parallel = FALSE, ncores, }else{ sims <- list() for(i in 1:nsim){ - sims[[i]] <- try(ExtinctionOrder(Network = network, Order = sample(1:network.size(network), size = SimNum), IS = IS, verbose = FALSE), silent = TRUE) + sims[[i]] <- try(ExtinctionOrder(Network = network, Order = sample(1:network.size(network), size = SimNum), + IS = IS, NetworkType = NetworkType, + Rewiring = Rewiring, RewiringDist = RewiringDist, + verbose = FALSE, RewiringProb = RewiringProb), silent = TRUE) try({sims[[i]]$simulation <- i}, silent = TRUE) if(verbose){setTxtProgressBar(ProgBar, i)} } @@ -686,12 +509,12 @@ RandomExtinctions <- function(Network, nsim = 10, parallel = FALSE, ncores, FullSims <- sims } - sims <- sims %>% group_by(NumExt) %>% summarise(AccSecExt_95CI = 1.96*sd(AccSecExt), AccSecExt_mean = mean(AccSecExt)) %>% mutate(Upper = AccSecExt_mean + AccSecExt_95CI, Lower = AccSecExt_mean - AccSecExt_95CI, Lower = ifelse(Lower < 0, 0, Lower)) + sims <- sims[!is.na(sims$SecExt), ] %>% group_by(NumExt) %>% summarise(AccSecExt_95CI = 1.96*sd(AccSecExt), AccSecExt_mean = mean(AccSecExt)) %>% mutate(Upper = AccSecExt_mean + AccSecExt_95CI, Lower = AccSecExt_mean - AccSecExt_95CI, Lower = ifelse(Lower < 0, 0, Lower)) ## plot output if(plot == TRUE){ - g <- ggplot(sims, aes_string(x = "NumExt", y = "AccSecExt_mean")) + geom_ribbon(aes_string(ymin = "Lower", ymax = "Upper"), fill = muted("red")) + geom_line() + ylab("Acc. Secondary extinctions") + xlab("Primary extinctions") + theme_bw() - g + g <- ggplot(sims, aes_string(x = "NumExt", y = "AccSecExt_mean")) + geom_ribbon(aes_string(ymin = "Lower", ymax = "Upper"), fill = scales::muted("red")) + geom_line() + ylab("Acc. Secondary extinctions") + xlab("Primary extinctions") + theme_bw() + print(g) } ## object output @@ -719,12 +542,10 @@ RandomExtinctions <- function(Network, nsim = 10, parallel = FALSE, ncores, #' with the observed extinction history. #' #' @examples -#' data("net") -#' History <- SimulateExtinctions(Network = net, Method = "Mostconnected") -#' -#' NullHyp <- RandomExtinctions(Network = net, nsim = 100, plot = TRUE) -#' -#' CompareExtinctions(Nullmodel = NullHyp[[1]], Hypothesis = History) +#' data("Less_Connected") +#' History <- SimulateExtinctions(Network = Less_Connected, Method = "Mostconnected", NetworkType = "Mutualistic") +#' NullHyp <- RandomExtinctions(Network = Less_Connected, nsim = 100, NetworkType = "Mutualistic") +#' CompareExtinctions(Nullmodel = NullHyp, Hypothesis = History) #' @importFrom broom tidy #' @importFrom ggplot2 aes #' @importFrom ggplot2 geom_line @@ -744,7 +565,7 @@ CompareExtinctions <- function(Nullmodel, Hypothesis){ g <- Nullmodel$graph + geom_line(aes(color = "blue")) g <- g + geom_point(data = Hypothesis, aes(y = AccSecExt), color = "black") + geom_line(data = Hypothesis, aes(y = AccSecExt, color = "black")) + scale_color_manual(name = "Comparison",values =c("black", "blue"), label = c("Observed","Null hypothesis")) } else { - g <- ggplot(Nullmodel$sims, aes(x = NumExt, y = AccSecExt_mean)) + geom_ribbon(aes_string(ymin = "Lower", ymax = "Upper"), fill = muted("red")) + geom_line() + ylab("Acc. Secondary extinctions") + xlab("Primary extinctions") + theme_bw() + g <- ggplot(Nullmodel$sims, aes(x = NumExt, y = AccSecExt_mean)) + geom_ribbon(aes_string(ymin = "Lower", ymax = "Upper"), fill = scales::muted("red")) + geom_line() + ylab("Acc. Secondary extinctions") + xlab("Primary extinctions") + theme_bw() g <- g + geom_point(data = Hypothesis$sims, aes(y = AccSecExt), color = "black") + geom_line(data = Hypothesis$sims, aes(y = AccSecExt, color = "black")) + scale_color_manual(name = "Comparison",values =c("black", "blue"), label = c("Observed","Null hypothesis")) g } diff --git a/R/NewExtinctionOrder.R b/R/NewExtinctionOrder.R new file mode 100644 index 0000000..15f6f66 --- /dev/null +++ b/R/NewExtinctionOrder.R @@ -0,0 +1,264 @@ +ExtinctionOrder <- function(Network, Order, IS = 0, + Rewiring = NULL, RewiringDist = NULL, + verbose = TRUE, clust.method = "cluster_infomap"){ + ## Setting up Objects for function run + Link_density <- Modularity <- Grado <- NULL + Network <- .DataInit(x = Network) + edgelist <- as.matrix.network.edgelist(Network,matrix.type="edgelist") #Prey - Predator + Conected <- data.frame(ID = 1:network.size(Network), Grado = degree(edgelist, c("total"))) + Conected1 <- Order + if(length(IS )== 1){ + IS <- rep(IS, network.size(Network)) + names(IS) <- get.vertex.attribute(Network, "vertex.names") + } + + ## Base net calculations + ### identify base interaction strengths per node + if(sum(IS) != 0){ + if(sum(get.edge.attribute(Network, "weight"), na.rm = TRUE) == 0){ + stop("Either your network does not contain any edges with weights or your network does not have the edge attribute `weight` required for calculation of extinctions based on relative interaction strength loss.") + } + net <- as.matrix.network.adjacency(Network, attrname = "weight") + netgraph <- suppressMessages(graph_from_adjacency_matrix(net, weighted = TRUE)) + strengthbasenet <- igraph::strength(netgraph) + } + Weight_mat <- as.matrix.network.adjacency(Network, attrname = "weight") + + ### identification of producers and top predators + indegreebasenet <- degree(Network, cmode = "indegree") + indegreebasenetzeros <- sum(degree(Network, cmode = "indegree") == 0) + indegreetopnetzeros <- sum(degree(Network, cmode = "outdegree") == 0) + Producers <- (1:length(degree(Network, cmode = "indegree")))[degree(Network, cmode = "indegree") == 0] + TopPredators <- (1:length(degree(Network, cmode = "outdegree")))[degree(Network, cmode = "outdegree") == 0] + + ### output object + DF <- data.frame(Spp = rep(NA, length(Order)), + S = rep(NA, length(Order)), + L = rep(NA, length(Order)), + C = rep(NA, length(Order)), + Link_density = rep(NA, length(Order)), + SecExt = rep(NA,length(Order)), + Pred_release = rep(NA,length(Order)), + Iso_nodes =rep (NA,length(Order))) + Secundaryext <- c() + Predationrel <- c() + accExt <- c() + totalExt <- c() + FinalExt <- list() + Conected3 <- c() + + ## Rewiring + if(!isFALSE(Rewiring)){ + ### default decay parameter + diag(RewiringDist)<- NA + if(is.null(Rewiring)){ + dist_10 <- quantile(RewiringDist, na.rm = TRUE, 0.1) # 10% distance quantile + decay <- 1/as.numeric(dist_10) # assuming mean rewiring capability lies at dist_10 + Rewiring <- function(x){1-pexp(x, rate = decay)} + } + + if(length(Rewiring )== 1){ + fun <- deparse1(Rewiring) + Rewiring <- rep(fun, network.size(Network)) + names(Rewiring) <- get.vertex.attribute(Network, "vertex.names") + } + # plot_seq <- seq(from = 0, to = max(RewiringDist[RewiringDist != 0], na.rm = TRUE), length = 1e3) + # plot(plot_seq, + # eval(str2lang(Rewiring[1]))(x) + # ) + } + + + ## Sequential extinction simulation + if(verbose){ProgBar <- txtProgressBar(max = length(Order), style = 3)} + for (i in 1:length(Order)){ + # print(i) + + ### creating temporary network representations and deleting vertices if they have been set to go extinct + if (length(accExt)==0){ + Temp <- Network + DF$Spp[i] <- Conected1[i] + delete.vertices(Temp, c(DF$Spp[1:i])) + } + if (length(accExt)>0){ + Temp <- Network + Temp <- delete.vertices(Temp, c(accExt)) + edgelist <- as.matrix.network.edgelist(Temp,matrix.type="edgelist") + # Conected2 <- data.frame(ID =1:network.size(Temp), Grado = degree(edgelist, c("total"))) + # for(j in sort(accExt)){ + # Conected2$ID <- ifelse(Conected2$ID < j, Conected2$ID, Conected2$ID + 1) + # } + + DF$Spp[i] <- Conected1[i] + Temp <- Network + + delete.vertices(Temp, unique(c(c(DF$Spp[1:i]),accExt))) + + } + # print(Temp) + + ### network metrics to output object + DF$S[i] <- network.size(Temp) + DF$L[i] <- network.edgecount(Temp) + DF$C[i] <- network.density(Temp) + DF$Link_density[i] <- DF$L[i]/DF$S[i] + + if(i > 1 ){ + if(DF$L[i-1] == 0){ + if(verbose){setTxtProgressBar(ProgBar, length(Order))} + warning(paste("All species in network went extinct through secondary extinction before all primary extinctions were simulated. This happened at extinction step", i-1, "out of", length(Order))) + break + } + } + + + ### calculating modularity + Networkclass = class(Temp) + if (Networkclass[1] == "matrix"){ + netgraph = graph_from_adjacency_matrix(Temp, mode = "directed", weighted = TRUE) + } + + if (Networkclass[1] == "network"){ + net = as.matrix.network.adjacency(Temp) + netgraph = suppressMessages(graph_from_adjacency_matrix(net, mode = "directed", weighted = TRUE)) + } + + if (clust.method == "cluster_edge_betweenness"){ + Membership = suppressWarnings(cluster_edge_betweenness(netgraph, weights = TRUE, directed = TRUE, edge.betweenness = TRUE, + merges = TRUE, bridges = TRUE, modularity = TRUE, membership = TRUE)) + } else if (clust.method == "cluster_spinglass"){ + spins = 107#network.size(Temp) + Membership = suppressWarnings(cluster_spinglass(netgraph, spins=spins)) #spins could be the Richness + }else if (clust.method == "cluster_label_prop"){ + Membership = suppressWarnings(cluster_label_prop(netgraph, weights = TRUE, initial = NULL, + fixed = NULL)) + }else if (clust.method == "cluster_infomap"){ + nb.trials = 107#network.size(Temp) + Membership = suppressWarnings(cluster_infomap(as.undirected(netgraph), + e.weights = E(netgraph)$weight, + v.weights = NULL, + nb.trials = nb.trials, + modularity = TRUE)) + + } else if (clust.method == "none"){ + Membership = NA + }else stop('Select a valid method for clustering. ?SimulateExtinction') + # if(is.na(Membership)){ + DF$Modularity[i] <- Membership$modularity + # }else{ + # DF$Modularity[i] <- suppressWarnings(modularity(Membership)) + # } + + ### identifying secondary extinctions + #### Producers + SecundaryextTemp <- (1:length(degree(Temp, cmode = "indegree")))[degree(Temp, cmode = "indegree") == 0] + for(j in sort(unique(c(c(DF$Spp[1:i]),accExt)))){ + SecundaryextTemp <- ifelse(SecundaryextTemp < j, SecundaryextTemp, SecundaryextTemp + 1) + } + Secundaryext <- SecundaryextTemp + Secundaryext <- Secundaryext[!(Secundaryext %in% Producers)] + DF$SecExt[i] <- length(Secundaryext) + + #### Predators + PredationrelTemp <- (1:length(degree(Temp, cmode = "outdegree")))[degree(Temp, cmode = "outdegree") == 0] + for(j in sort(unique(c(c(DF$Spp[1:i]),accExt)))){ + PredationrelTemp <- ifelse(PredationrelTemp < j, PredationrelTemp, PredationrelTemp + 1) + } + Predationrel <- PredationrelTemp + Predationrel <- Predationrel[!(Predationrel %in% TopPredators)] + DF$Pred_release[i]<- length(Predationrel) + DF$Iso_nodes[i] <- sum(degree(Temp) == 0) + + ### rewiring + accExt <- append(accExt, DF$Spp[1:i]) + if(!isFALSE(Rewiring)){ + Rewiring_df <- data.frame(Direction = NA, + Species = NA, + NewPartner = NA, + LostPartner = NA, + IS = NA) + Rewiring_df <- na.omit(Rewiring_df) + #### loop over all deleted vertices and the connections lost because of their exclusion + for(Iter_PrimaryExt in 1:length(accExt)){ + # Iter_PrimaryExt = 1 + LostPartner <- get.vertex.attribute(Network, "vertex.names")[accExt[Iter_PrimaryExt]] # name of primary extinction species + LostIS <- Weight_mat[, LostPartner] # lost interaction strength with nodes now slated for secondary extinction + Direction <- 1 # identify column-driven loss + if(sum(abs(LostIS))==0){# if the primary species is not an animal, LostIS will be filled with 0s, so we need to look for LosTIS in other orientiation in Weight_mat + LostIS <- Weight_mat[LostPartner, ] + Direction <- 2 # identify row-driven loss + } + + + for(Iter_LostIS in 1:length(LostIS)){ ## looping over all species that were linked to the current primary extinction + # Iter_LostIS = 1 + LostPartnerSim <- eval(str2lang(Rewiring[Iter_LostIS]))(dist_mat[,LostPartner]) # probability of rewiring too each node in network given rewiring function and species similraity + RewiringCandidates <- LostPartnerSim[LostPartnerSim > 0.3 & names(LostPartnerSim) %in% get.vertex.attribute(Temp, "vertex.names")] # rewiring probability for nodes still in temporary network and having a higher rewiring probability than 0.3 + RewiredPartner <- names(which.max(RewiringCandidates)) # most likely rewiring partner + if(!is.null(RewiredPartner)){ # if a rewired partner has been found + Rewiring_df <- rbind(Rewiring_df, + data.frame(Direction = Direction, + Species = names(LostIS[Iter_LostIS]), + NewPartner = RewiredPartner, + LostPartner = LostPartner, + IS = LostIS[Iter_LostIS]) + ) + } + } + } + + + #### shift interaction weights in Weight_mat + for(Iter_Rewiring in 1:nrow(Rewiring_df)){ + # Iter_Rewiring = 1 + ## assigning shifted interaction strength + ColSpec <- Rewiring_df[Iter_Rewiring,4-Rewiring_df[Iter_Rewiring,"Direction"]] + RowSpec <- Rewiring_df[Iter_Rewiring,1+Rewiring_df[Iter_Rewiring,"Direction"]] + Weight_mat[RowSpec, ColSpec] <- Weight_mat[RowSpec, ColSpec] + Rewiring_df[Iter_Rewiring,"IS"] + ## deleting shiften interaction strength + + ColLost <- ifelse(Rewiring_df[Iter_Rewiring, "Direction"] == 1, + Rewiring_df[Iter_Rewiring, "LostPartner"], + Rewiring_df[Iter_Rewiring, "Species"]) + RowLost <- ifelse(Rewiring_df[Iter_Rewiring, "Direction"] == 1, + Rewiring_df[Iter_Rewiring, "Species"], + Rewiring_df[Iter_Rewiring, "LostPartner"]) + Weight_mat[RowLost, ColLost] <- 0 + } + #### establishing rewired network and deleting primary extinction nodes + Network <- as.network(Weight_mat, matrix.type = "adjacency", ignore.eval=FALSE, names.eval='weight') + Temp <- Network + delete.vertices(Temp, unique(c(c(DF$Spp[1:i]),accExt))) + } + + #### Relative Interaction Strength loss + if(sum(IS) == 0){ + Secundaryext <- get.vertex.attribute(Temp, "vertex.names")[which(degree(Temp) == 0)] + Secundaryext <- match(Secundaryext, get.vertex.attribute(Network, "vertex.names")) + }else{ + AbsIS <- igraph::strength(suppressMessages(graph_from_adjacency_matrix( + as.matrix.network.adjacency(Temp, attrname = "weight"), + weighted = TRUE) + )) + RelISloss <- AbsIS / strengthbasenet[names(strengthbasenet) %in% get.vertex.attribute(Temp, "vertex.names")] + Secundaryext <- which(AbsIS == 0 | RelISloss < IS[match(names(RelISloss), names(IS))]) + Secundaryext <- match(names(Secundaryext), get.vertex.attribute(Network, "vertex.names")) + } + DF$SecExt[i] <- length(Secundaryext) + + ### Return of objects + FinalExt[[i]] <- (Secundaryext) + accExt <- append(accExt, DF$Spp[1:i]) + accExt <- unique(append(accExt,Secundaryext)) + if(verbose){setTxtProgressBar(ProgBar, i)} + + } + DF <- DF[complete.cases(DF),] + DF$AccSecExt <- cumsum(DF$SecExt) + DF$NumExt <- 1:nrow(DF) + DF$TotalExt <- DF$AccSecExt + DF$NumExt + DF <- relocate(DF, Modularity, .after = Link_density) + class(DF) <- c("data.frame", "SimulateExt") + return(list(sims = DF, + Network = Temp)) +} diff --git a/R/NewSimulateExtinctions.R b/R/NewSimulateExtinctions.R new file mode 100644 index 0000000..e7395d8 --- /dev/null +++ b/R/NewSimulateExtinctions.R @@ -0,0 +1,27 @@ +SimulateExtinctions <- function(Network, Method, Order = NULL, + clust.method = "cluster_infomap", + IS = 0, + Rewiring = FALSE, RewiringDist = NULL, + verbose = TRUE){ + Network <- .DataInit(x = Network) + + if(!is.null(Order)){Method <- "Ordered"} + + '%ni%'<- Negate('%in%') + if(Method %ni% c("Mostconnected", "Ordered")) stop('Choose the right method. See ?SimulateExtinction.') + + if(Method == "Mostconnected"){ + edgelist <- as.matrix.network.edgelist(Network,matrix.type="edgelist") #Prey - Predator + Conected <- data.frame(ID = 1:network.size(Network), Grado = degree(edgelist, c("total"))) + Conected <- arrange(Conected, desc(Grado)) + DF <- ExtinctionOrder(Network = Network, Order = Conected$ID, clust.method = clust.method, + IS = IS, Rewiring = Rewiring, RewiringDist = RewiringDist, verbose = verbose) + } + if(Method == "Ordered"){ + DF <- ExtinctionOrder(Network = Network, Order = Order, clust.method = clust.method, + IS = IS, Rewiring = Rewiring, RewiringDist = RewiringDist, verbose = verbose) + } + + return(DF) +} + diff --git a/data/Less_Connected.rda b/data/Less_Connected.rda index 6b30a01..e1a528e 100644 Binary files a/data/Less_Connected.rda and b/data/Less_Connected.rda differ diff --git a/data/More_Connected.rda b/data/More_Connected.rda index ed9f905..152f34b 100644 Binary files a/data/More_Connected.rda and b/data/More_Connected.rda differ diff --git a/data/dist.rda b/data/dist.rda new file mode 100644 index 0000000..0c14f60 Binary files /dev/null and b/data/dist.rda differ diff --git a/man/CompareExtinctions.Rd b/man/CompareExtinctions.Rd index 48d9eea..b1dd98d 100644 --- a/man/CompareExtinctions.Rd +++ b/man/CompareExtinctions.Rd @@ -22,11 +22,9 @@ with a null hypothesis generated by the RandomExtinctions function it is importa RandomExtinctions is in plot = T. } \examples{ -data("net") -History <- SimulateExtinctions(Network = net, Method = "Mostconnected") - -NullHyp <- RandomExtinctions(Network = net, nsim = 100, plot = TRUE) - +data("Less_Connected") +History <- SimulateExtinctions(Network = Less_Connected, Method = "Mostconnected", NetworkType = "Mutualistic") +NullHyp <- RandomExtinctions(Network = Less_Connected, nsim = 100, NetworkType = "Mutualistic") CompareExtinctions(Nullmodel = NullHyp, Hypothesis = History) } \author{ diff --git a/man/ExtinctionOrder.Rd b/man/ExtinctionOrder.Rd index 5cdb5ed..68573c6 100644 --- a/man/ExtinctionOrder.Rd +++ b/man/ExtinctionOrder.Rd @@ -4,24 +4,57 @@ \alias{ExtinctionOrder} \title{Extinctions analysis from custom order} \usage{ -ExtinctionOrder(Network, Order) +ExtinctionOrder( + Network, + Order, + IS = 0, + Rewiring = NULL, + RewiringDist = NULL, + verbose = TRUE, + clust.method = "cluster_infomap" +) } \arguments{ -\item{Network}{a network of class network} +\item{Network}{a network representation as a an adjacency matrix, edgelist, +or a network object} -\item{Order}{Vector with the order of extinctions by ID} +\item{Order}{a numeric vector indexing order of primary extinctions. For Method = Mostconnected Order must be NULL. If Order is not NULL, Method is internally forced to be Ordered.} + +\item{IS}{either numeric or a named vector of numerics. Identifies the threshold of relative interaction strength which species require before being considered secondarily extinct (i.e. IS = 0.7 leads to removal of all nodes which lose 30% of their interaction strength in the Network argument). If a named vector, names must correspond to vertex names in Network argument.} + +\item{Rewiring}{either a function or a named vector of functions. Signifies how rewiring probabilities are calculated from the RewiringDist argument. If FALSE, no rewiring is carried out.} + +\item{RewiringDist}{a numeric matrix of NxN dimension (N... number of nodes in Network). Contains, for example, phylogenetic or functional trait distances between nodes in Network which are used by the Rewiring argument to calculate rewiring probabilities.} + +\item{verbose}{Logical. Whether to report on function progress or not.} + +\item{clust.method}{a character with the options cluster_edge_betweenness, cluster_spinglass, +cluster_label_prop or cluster_infomap, defaults to cluster_infomap} + +\item{Method}{a character with the options Mostconnected and Ordered} + +\item{NetworkType}{a character with the options Trophic and Mutualistic - is used to calculate secondary extinctions.} + +\item{RewiringProb}{a numeric which identifies the threshold at which to assume rewiring potential is met.} } \value{ -exports data frame with the characteristics of the network after every -extintion, and a graph with the mean and 95% interval +exports list containing a data frame with the characteristics of the network after every extinction and a network object containing the final network. The resulting data frame contains 11 columns that incorporate the topological index, the secondary extinctions, predation release, and total extinctions of the network in each primary extinction. } \description{ -It takes a network, and extinguishes nodes using a custom order, -then it calculates the secondary extinctions and plots the accumulated -secondary extinctions. +This function takes a network and eliminates nodes using a custom order. Subsequently, secondary extinctions are tallied up. Secondary extinction severity can be targeted by manipulating the node-dependency on network edges (IS) and node-rewiring potential upon loss of links (Rewiring). +} +\details{ +When NetworkType = Trophic, secondary extinctions only occur for any predator, but not producers. If NetworkType = Mutualistic, secondary extinctions occur for all species in the network. + +When clust.method = cluster_edge_betweenness computes the network modularity using cluster_edge_betweenness methods from igraph to detect communities +When clust.method = cluster_spinglass computes the network modularity using cluster_spinglass methods from igraph to detect communities, here the number of spins are equal to the network size +When clust.method = cluster_label_prop computes the network modularity using cluster_label_prop methods from igraph to detect communities +When clust.method = cluster_infomap computes the network modularity using cluster_infomap methods from igraph to detect communities, here the number of nb.trials are equal to the network size } \author{ Derek Corcoran M. Isidora Ávila-Thieme + +Erik Kusch } diff --git a/man/Less_Connected.Rd b/man/Less_Connected.Rd index d315905..28b4757 100644 --- a/man/Less_Connected.Rd +++ b/man/Less_Connected.Rd @@ -11,8 +11,8 @@ a network Less_Connected } \description{ -A trophic network with 30 species and 47 trophic interactions. -This foodweb has a connectance of 0.03 +A network with 30 species and 47 interactions. +This network has a connectance of 0.03 } \seealso{ \code{\link{More_Connected}} diff --git a/man/Mostconnected.Rd b/man/Mostconnected.Rd deleted file mode 100644 index 9129f6f..0000000 --- a/man/Mostconnected.Rd +++ /dev/null @@ -1,34 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/Extintions.R -\name{Mostconnected} -\alias{Mostconnected} -\title{Extinctions analysis from most connected to less connected nodes in the network} -\usage{ -Mostconnected(Network) -} -\arguments{ -\item{Network}{a trophic network of class network} -} -\value{ -exports data frame with the characteristics of the network after every -extintion. The resulting data frame contains 11 columns that incorporate the -topological index, the secondary extinctions, predation release, and total extinctions of the network -in each primary extinction. -} -\description{ -It takes a network and it calculates wich node is the most connected -of the network, using total degree. Then remove the most connected node, -and calculates the the topological indexes of the network and the number of -secundary extintions (how many species have indegree 0, without considered -primary producers). After that, remove the nodes that were secondarily extinct -in the previous step and recalculate which is the new most connected -node and so on, until the number of links in the network is zero. -} -\seealso{ -[NetworkExtinction::ExtinctionOrder()] -} -\author{ -Derek Corcoran - -M. Isidora Ávila-Thieme -} diff --git a/man/RandomExtinctions.Rd b/man/RandomExtinctions.Rd index 3e3f743..ae622df 100644 --- a/man/RandomExtinctions.Rd +++ b/man/RandomExtinctions.Rd @@ -7,35 +7,66 @@ RandomExtinctions( Network, nsim = 10, + Record = FALSE, + plot = FALSE, + SimNum = NULL, + NetworkType = "Trophic", + clust.method = "cluster_infomap", parallel = FALSE, ncores, - Record = F, - plot = F + IS = 0, + Rewiring = FALSE, + RewiringDist, + RewiringProb = 0.5, + verbose = TRUE ) } \arguments{ -\item{Network}{a trophic network of class network} +\item{Network}{a network representation as a an adjacency matrix, edgelist, +or a network object} -\item{nsim}{number of simulations} +\item{nsim}{numeric, number of simulations} + +\item{Record}{logical, if TRUE, records every simulation and you can read the +raw results in the object FullSims} + +\item{plot}{logical if TRUE, will add a graph to the results} + +\item{SimNum}{numeric, how many nodes to register for primary extinction. By default sets all of them.} + +\item{NetworkType}{a character with the options Trophic and Mutualistic - is used to calculate secondary extinctions.} + +\item{clust.method}{a character with the options cluster_edge_betweenness, cluster_spinglass, +cluster_label_prop or cluster_infomap, defaults to cluster_infomap} \item{parallel}{if TRUE, it will use parallel procesing, if FALSE (default) it will run sequentially} -\item{ncores}{number of cores to use if using parallel procesing} +\item{ncores}{numeric, number of cores to use if using parallel procesing} -\item{Record}{logical, if TRUE, records every simulation and you can read the -raw results in the object FullSims} +\item{IS}{either numeric or a named vector of numerics. Identifies the threshold of relative interaction strength which species require before being considered secondarily extinct (i.e. IS = 0.7 leads to removal of all nodes which lose 30% of their interaction strength in the Network argument). If a named vector, names must correspond to vertex names in Network argument.} + +\item{Rewiring}{either a function or a named vector of functions. Signifies how rewiring probabilities are calculated from the RewiringDist argument. If FALSE, no rewiring is carried out.} + +\item{RewiringDist}{a numeric matrix of NxN dimension (N... number of nodes in Network). Contains, for example, phylogenetic or functional trait distances between nodes in Network which are used by the Rewiring argument to calculate rewiring probabilities.} + +\item{RewiringProb}{a numeric which identifies the threshold at which to assume rewiring potential is met.} -\item{plot}{logical if true, will add a graph to the results} +\item{verbose}{Logical. Whether to report on function progress or not.} } \value{ -exports data frame with the characteristics of the network after every -extintion, and a graph with the mean and 95% interval +exports list containing a data frame with the characteristics of the network after every extinction, a network object containing the final network, and a graph with the mean and 95% interval. The resulting data frame contains 11 columns that incorporate the topological index, the secondary extinctions, predation release, and total extinctions of the network in each primary extinction. } \description{ -Generates a null model by generating random extinction histories and calculating -the mean and standard deviation of the accumulated secondary extinctions developed -by making n random extinction histories +Generates a null model by generating random extinction histories and calculating the mean and standard deviation of the accumulated secondary extinctions developed by making n random extinction histories. +} +\details{ +When NetworkType = Trophic, secondary extinctions only occur for any predator, but not producers. If NetworkType = Mutualistic, secondary extinctions occur for all species in the network. + +When clust.method = cluster_edge_betweenness computes the network modularity using cluster_edge_betweenness methods from igraph to detect communities +When clust.method = cluster_spinglass computes the network modularity using cluster_spinglass methods from igraph to detect communities, here the number of spins are equal to the network size +When clust.method = cluster_label_prop computes the network modularity using cluster_label_prop methods from igraph to detect communities +When clust.method = cluster_infomap computes the network modularity using cluster_infomap methods from igraph to detect communities, here the number of nb.trials are equal to the network size } \examples{ #first example @@ -55,4 +86,6 @@ RandomExtinctions(Network = More_Connected, nsim = 20, parallel = TRUE, ncores = Derek Corcoran M. Isidora Ávila-Thieme + +Erik Kusch } diff --git a/man/SimulateExtinctions.Rd b/man/SimulateExtinctions.Rd index 0c21812..751934a 100644 --- a/man/SimulateExtinctions.Rd +++ b/man/SimulateExtinctions.Rd @@ -2,55 +2,61 @@ % Please edit documentation in R/Extintions.R \name{SimulateExtinctions} \alias{SimulateExtinctions} -\title{Extinctions analysis for trophic networks} +\title{Extinctions analysis for ecological networks} \usage{ SimulateExtinctions( Network, Method, Order = NULL, - clust.method = "cluster_infomap" + clust.method = "cluster_infomap", + IS = 0, + Rewiring = FALSE, + RewiringDist = NULL, + verbose = TRUE ) } \arguments{ -\item{Network}{a network representation as a an adyacency matrix, edgelist, +\item{Network}{a network representation as a an adjacency matrix, edgelist, or a network object} -\item{Method}{a character with the options Mostconnected, Oredered, or Random} +\item{Method}{a character with the options Mostconnected and Ordered} -\item{Order}{this should be NULL, unless using the Ordered method, in that case -it should be a vector with the order of extinctions by ID} +\item{Order}{a numeric vector indexing order of primary extinctions. For Method = Mostconnected Order must be NULL. If Order is not NULL, Method is internally forced to be Ordered.} \item{clust.method}{a character with the options cluster_edge_betweenness, cluster_spinglass, cluster_label_prop or cluster_infomap, defaults to cluster_infomap} + +\item{IS}{either numeric or a named vector of numerics. Identifies the threshold of relative interaction strength which species require before being considered secondarily extinct (i.e. IS = 0.7 leads to removal of all nodes which lose 30% of their interaction strength in the Network argument). If a named vector, names must correspond to vertex names in Network argument.} + +\item{Rewiring}{either a function or a named vector of functions. Signifies how rewiring probabilities are calculated from the RewiringDist argument. If FALSE, no rewiring is carried out.} + +\item{RewiringDist}{a numeric matrix of NxN dimension (N... number of nodes in Network). Contains, for example, phylogenetic or functional trait distances between nodes in Network which are used by the Rewiring argument to calculate rewiring probabilities.} + +\item{verbose}{Logical. Whether to report on function progress or not.} + +\item{NetworkType}{a character with the options Trophic and Mutualistic - is used to calculate secondary extinctions.} + +\item{RewiringProb}{a numeric which identifies the threshold at which to assume rewiring potential is met.} } \value{ -exports data frame with the characteristics of the network after every -extintion. The resulting data frame contains 11 columns that incorporate the -topological index, the secondary extinctions, predation release, and total extinctions of the network -in each primary extinction. +exports list containing a data frame with the characteristics of the network after every extinction and a network object containing the final network. The resulting data frame contains 11 columns that incorporate the topological index, the secondary extinctions, predation release, and total extinctions of the network in each primary extinction. } \description{ The SimulateExtinctions function, can be used to test how the order of species -extinctions might affect the stability of the network by comparing The extintion history +extinctions, species-dependency on existing interaction strength, and rewiring potential might affect the stability of the network by comparing The extinction history and checking for secondary extinctions. } \details{ -When method is Mostconnected, it takes a network and it calculates wich node is the most connected -of the network, using total degree. Then remove the most connected node, -and calculates the the topological indexes of the network and the number of -secundary extintions (how many species have indegree 0, without considered -primary producers). After that, remove the nodes that were secondarily extinct -in the previous step and recalculate which is the new most connected -node and so on, until the number of links in the network is zero. - -When method is Ordered, it takes a network, and extinguishes nodes using a custom order, -then it calculates the secondary extinctions and plots the accumulated -secondary extinctions. +When method is Mostconnected, the function takes the network and calculates which node is the most connected of the network, using total degree. Then remove the most connected node, and calculates the the topological indexes of the network and the number of secondary extinctions. + +When method is Ordered, it takes a network, and extinguishes nodes using a custom order, then it calculates the secondary extinctions and plots the accumulated secondary extinctions. + +When NetworkType = Trophic, secondary extinctions only occur for any predator, but not producers. If NetworkType = Mutualistic, secondary extinctions occur for all species in the network. When clust.method = cluster_edge_betweenness computes the network modularity using cluster_edge_betweenness methods from igraph to detect communities -When clust.method = cluster_spinglass computes the network modularity using cluster_spinglass methods from igraph to detect communities, here the number of spins are equal to the nerwork size +When clust.method = cluster_spinglass computes the network modularity using cluster_spinglass methods from igraph to detect communities, here the number of spins are equal to the network size When clust.method = cluster_label_prop computes the network modularity using cluster_label_prop methods from igraph to detect communities -When clust.method = cluster_infomap computes the network modularity using cluster_infomap methods from igraph to detect communities, here the number of nb.trials are equal to the nerwork size +When clust.method = cluster_infomap computes the network modularity using cluster_infomap methods from igraph to detect communities, here the number of nb.trials are equal to the network size } \examples{ # Mostconnected example @@ -67,9 +73,34 @@ Method = "Ordered" , clust.method = "cluster_infomap") data("net") SimulateExtinctions(Network = net, Order = c(2,8,9), Method = "Ordered", clust.method = "cluster_infomap") + +#Network-Dependency Example +data("net") +SimulateExtinctions(Network = net, Order = c(2,8), IS = 0.3, +Method = "Ordered", clust.method = "cluster_infomap") + + #Rewiring +data("net") +data(dist) +SimulateExtinctions(Network = net, Order = c(2,8), IS = 0.3, +Rewiring = function(x){1-pexp(x, rate = 1/0.5)}, # assuming an exponential decline in rewiring potential as values in RewiringDist increase +RewiringDist = dist, # distance matrix +RewiringProb = 0.2, # low threshold for rewiring potential +Method = "Ordered", clust.method = "cluster_infomap") + +#Rewiring, assuming dist contains probabilities +#' data("net") +data(dist) +SimulateExtinctions(Network = net, Order = c(2,8), IS = 0.3, +Rewiring = function(x){x}, # no changes to the RewiringDist object means +RewiringDist = dist, RewiringProb = 0.2, +Method = "Ordered", clust.method = "cluster_infomap") + } \author{ Derek Corcoran M. Isidora Ávila-Thieme + +Erik Kusch } diff --git a/man/dist.Rd b/man/dist.Rd new file mode 100644 index 0000000..6e8a827 --- /dev/null +++ b/man/dist.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Data.R +\docType{data} +\name{dist} +\alias{dist} +\title{A toymodel distance matrix} +\format{ +a distance matrix +} +\usage{ +dist +} +\description{ +A distance matrix used for demonstration of rewiring capabilities +} +\keyword{datasets} diff --git a/man/dot-ExtinctionOrder.Rd b/man/dot-ExtinctionOrder.Rd deleted file mode 100644 index c4f2f5b..0000000 --- a/man/dot-ExtinctionOrder.Rd +++ /dev/null @@ -1,30 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/Extintions.R -\name{.ExtinctionOrder} -\alias{.ExtinctionOrder} -\title{Extinctions analysis from custom order} -\usage{ -.ExtinctionOrder(Network, Order, clust.method = "cluster_infomap") -} -\arguments{ -\item{Network}{a network of class network} - -\item{Order}{Vector with the order of extinctions by ID} - -\item{clust.method}{a character with the options cluster_edge_betweenness, cluster_spinglass, -cluster_label_prop or cluster_infomap} -} -\value{ -exports data frame with the characteristics of the network after every -extintion, and a graph with the mean and 95% interval -} -\description{ -It takes a network, and extinguishes nodes using a custom order, -then it calculates the secondary extinctions and plots the accumulated -secondary extinctions. -} -\author{ -Derek Corcoran - -M. Isidora Ávila-Thieme -} diff --git a/man/dot-Mostconnected.Rd b/man/dot-Mostconnected.Rd deleted file mode 100644 index 6462d09..0000000 --- a/man/dot-Mostconnected.Rd +++ /dev/null @@ -1,36 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/Extintions.R -\name{.Mostconnected} -\alias{.Mostconnected} -\title{Extinctions analysis from most connected to less connected nodes in the network} -\usage{ -.Mostconnected(Network, clust.method = "cluster_infomap") -} -\arguments{ -\item{Network}{a trophic network of class network} - -\item{clust.method}{a character with the options cluster_edge_betweenness, cluster_spinglass, cluster_label_prop or cluster_infomap} -} -\value{ -exports data frame with the characteristics of the network after every -extintion. The resulting data frame contains 11 columns that incorporate the -topological index, the secondary extinctions, predation release, and total extinctions of the network -in each primary extinction. -} -\description{ -It takes a network and it calculates wich node is the most connected -of the network, using total degree. Then remove the most connected node, -and calculates the the topological indexes of the network and the number of -secundary extintions (how many species have indegree 0, without considered -primary producers). After that, remove the nodes that were secondarily extinct -in the previous step and recalculate which is the new most connected -node and so on, until the number of links in the network is zero. -} -\seealso{ -[NetworkExtinction::ExtinctionOrder()] -} -\author{ -Derek Corcoran - -M. Isidora Ávila-Thieme -} diff --git a/vignettes/How_to_use_the_NetworkExtinction_Package.Rmd b/vignettes/How_to_use_the_NetworkExtinction_Package.Rmd index cfb92c7..60b3c2b 100644 --- a/vignettes/How_to_use_the_NetworkExtinction_Package.Rmd +++ b/vignettes/How_to_use_the_NetworkExtinction_Package.Rmd @@ -189,7 +189,7 @@ Test <- RandomExtinctions(Network= net, nsim= 100, plot = T) Test$graph ``` - + ### Comparison of Null hypothesis with other extinction histories The *RandomExtinctons* function generates a null hypothesis for us to compare it with either an extinction history generated by the *ExtinctionOrder* function or the *Mostconnected* function. In order to compare the expected extinctions developed by our null hypothesis with the observed extinction history, we developed the *CompareExtinctions* function. The way to use this last function is to first create the extinction history and the null hypothesis, and then the *CompareExtinctins* function to compare both extinction histories. @@ -203,7 +203,7 @@ NullHyp <- RandomExtinctions(Network = net, nsim = 100, plot = T) Comparison <- CompareExtinctions(Nullmodel = NullHyp, Hypothesis = History) ``` -The result will be a graph (Figue 6) with a dashed line showing the observed extinction history and a solid line showing the expected value of secondary extinctions randomly generated. +The result will be a graph (Figue 6) with a dashed line showing the observed extinction history and a solid line showing the expected value of secondary extinctions randomly generated. ```{r, echo=FALSE, fig.cap= "Figure 6. The resulting graph of the CompareExtinctions function, where the dashed line shows the observed extinction history, and a solid line shows the expected value of secondary extinctions originated at random"} Comparison @@ -227,7 +227,7 @@ ExtinctionPlot(History = history, Variable = "Link_density") ## Degree distribution function -The *DegreeDistribution* function calculates the cumulative distribution of the number of links that each species in the food network has [@estrada2007food]. Then, the observed distribution is fitted to the exponential, and power law models. +The *DegreeDistribution* function calculates the cumulative distribution of the number of links that each species in the food network has [@estrada2007food]. Then, the observed distribution is fitted to the exponential, and power law models. The results of this function are shown in figure 9 and table 4. The graph shows the observed degree distribution in a log log scale fitting the three models mentioned above, for this example we use an example dataset of Chilean litoral rocky shores [@kefi2015network]. The table shows the fitted model information ordered by descending AIC, that is, the model in the first row is the most probable distribution, followed by the second an finally the third distribution in this case (Table 3), the Exponential distribution would be the best model, followed by the Power law model.